Index: palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model_mod.f90	(revision 3606)
+++ palm/trunk/SOURCE/radiation_model_mod.f90	(revision 3607)
@@ -28,4 +28,7 @@
 ! -----------------
 ! $Id$
+! Output of radiation-related quantities migrated to radiation_model_mod.
+! 
+! 3589 2018-11-30 15:09:51Z suehring
 ! Remove erroneous UTF encoding
 ! 
@@ -545,5 +548,5 @@
                message_string, plant_canopy, pt_surface,                       &
                rho_surface, simulated_time, spinup_time, surface_pressure,     &
-               time_since_reference_point, urban_surface
+               time_since_reference_point, urban_surface, varnamelength
 
     USE cpulog,                                                                &
@@ -899,4 +902,12 @@
     INTEGER(iwp)                                   ::  nwalls           !< number of wall surfaces in local processor
 
+!-- indices needed for RTM netcdf output subroutines
+    INTEGER(iwp), PARAMETER                        :: nd = 5
+    CHARACTER(LEN=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
+    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     :: dirint_u = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)
+    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     :: dirint_l = (/ iup_l, isouth_l, inorth_l, iwest_l, ieast_l /)
+    INTEGER(iwp), DIMENSION(0:nd-1)                :: dirstart
+    INTEGER(iwp), DIMENSION(0:nd-1)                :: dirend
+
 !-- indices and sizes of urban and land surface models
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfl_l          !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x]
@@ -1065,4 +1076,22 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  rt2_dist
 
+!-- arrays for time averages
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfradnet_av    !< average of net radiation to local surface including radiation from reflections
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinsw_av      !< average of sw radiation falling to local surface including radiation from reflections
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlw_av      !< average of lw radiation falling to local surface including radiation from reflections
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdir_av   !< average of direct sw radiation falling to local surface
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdif_av   !< average of diffuse sw radiation from sky and model boundary falling to local surface
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwdif_av   !< average of diffuse lw radiation from sky and model boundary falling to local surface
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswref_av   !< average of sw radiation falling to surface from reflections
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwref_av   !< average of lw radiation falling to surface from reflections
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutsw_av     !< average of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutlw_av     !< average of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins_av       !< average of array of residua of sw radiation absorbed in surface after last reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinl_av       !< average of array of residua of lw radiation absorbed in surface after last reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinlw_av       !< Average of pcbinlw
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinsw_av       !< Average of pcbinsw
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdir_av    !< Average of pcbinswdir
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdif_av    !< Average of pcbinswdif
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswref_av    !< Average of pcbinswref
 
 
@@ -1208,10 +1237,5 @@
            skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,&
            zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon,       &
-           nrefsteps, nsvfl, svf,                                              &
-           svfsurf, surfinsw, surfinlw, surfins, surfinl, surfinswdir,         &
-           surfinswdif, surfoutsw, surfoutlw, surfinlwdif, rad_sw_in_dir,      &
-           rad_sw_in_diff, rad_lw_in_diff, surfouts, surfoutl, surfoutsl,      &
-           surfoutll, idir, jdir, kdir, id, iz, iy, ix,                        &
-           surf, surfl, nsurfl, pcbinswdir, pcbinswdif, pcbinsw, pcbinlw,      &
+           idir, jdir, kdir, id, iz, iy, ix,                                   &
            iup_u, inorth_u, isouth_u, ieast_u, iwest_u,                        &
            iup_l, inorth_l, isouth_l, ieast_l, iwest_l,                        &
@@ -1219,6 +1243,6 @@
            idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct,                             &
            radiation_interactions, startwall, startland, endland, endwall,     &
-           skyvf, skyvft, radiation_interactions_on, average_radiation, npcbl, &
-           pcbl
+           skyvf, skyvft, radiation_interactions_on, average_radiation
+
 
 #if defined ( __rrtmg )
@@ -1263,5 +1287,5 @@
 !> Check data output for radiation model
 !------------------------------------------------------------------------------!
-    SUBROUTINE radiation_check_data_output( var, unit, i, ilen, k )
+    SUBROUTINE radiation_check_data_output( variable, unit, i, ilen, k )
  
  
@@ -1271,95 +1295,127 @@
        IMPLICIT NONE
 
-       CHARACTER (LEN=*) ::  unit     !< 
-       CHARACTER (LEN=*) ::  var      !<
-
-       INTEGER(iwp) :: i
+       CHARACTER (LEN=*) ::  unit          !<
+       CHARACTER (LEN=*) ::  variable      !<
+
+       INTEGER(iwp) :: i, j, k, l
        INTEGER(iwp) :: ilen
-       INTEGER(iwp) :: k
-
-       SELECT CASE ( TRIM( var ) )
-
-          CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_lw_in', 'rad_lw_out', &
-                 'rad_sw_cs_hr', 'rad_sw_hr', 'rad_sw_in', 'rad_sw_out'  )
-             IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
-                message_string = '"output of "' // TRIM( var ) // '" requi' // &
-                                 'res radiation = .TRUE. and ' //              &
-                                 'radiation_scheme = "rrtmg"'
-                CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
-             ENDIF
-             unit = 'K/h'     
-
-          CASE ( 'rad_net*', 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*',      &
-                 'rrtm_asdir*', 'rad_lw_in*', 'rad_lw_out*', 'rad_sw_in*',     &
-                 'rad_sw_out*')
-             IF ( i == 0 .AND. ilen == 0 .AND. k == 0)  THEN
-                ! Workaround for masked output (calls with i=ilen=k=0)
-                unit = 'illegal'
-                RETURN
-             ENDIF
-             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
-                message_string = 'illegal value for data_output: "' //         &
-                                 TRIM( var ) // '" & only 2d-horizontal ' //   &
-                                 'cross sections are allowed for this value'
-                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
-             ENDIF
-             IF (  .NOT.  radiation  .OR.  radiation_scheme /= "rrtmg" )  THEN
-                IF ( TRIM( var ) == 'rrtm_aldif*'  .OR.                        &
-                     TRIM( var ) == 'rrtm_aldir*'  .OR.                        &
-                     TRIM( var ) == 'rrtm_asdif*'  .OR.                        &
-                     TRIM( var ) == 'rrtm_asdir*'      )                       &
-                THEN
-                   message_string = 'output of "' // TRIM( var ) // '" require'&
-                                    // 's radiation = .TRUE. and radiation_sch'&
-                                    // 'eme = "rrtmg"'
-                   CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
-                ENDIF
-             ENDIF
-
-             IF ( TRIM( var ) == 'rad_net*'      ) unit = 'W/m2'   
-             IF ( TRIM( var ) == 'rad_lw_in*'    ) unit = 'W/m2'
-             IF ( TRIM( var ) == 'rad_lw_out*'   ) unit = 'W/m2'
-             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 = '' 
-             IF ( TRIM( var ) == 'rrtm_asdif*'   ) unit = '' 
-             IF ( TRIM( var ) == 'rrtm_asdir*'   ) unit = '' 
-
-          CASE ( 'rad_mrt', 'rad_mrt_sw', 'rad_mrt_lw'  )
-
-             IF ( i == 0 .AND. ilen == 0 .AND. k == 0)  THEN
-                ! Workaround for masked output (calls with i=ilen=k=0)
-                unit = 'illegal'
-                RETURN
-             ENDIF
-
-             IF ( .NOT.  radiation ) THEN
+       CHARACTER(LEN=varnamelength) :: var          !< TRIM(variable)
+
+       var = TRIM(variable)
+
+!--    first process diractional variables
+       IF ( var(1:12) == 'rtm_rad_net_'  .OR.  var(1:13) == 'rtm_rad_insw_'  .OR.        &
+            var(1:13) == 'rtm_rad_inlw_'  .OR.  var(1:16) == 'rtm_rad_inswdir_'  .OR.    &
+            var(1:16) == 'rtm_rad_inswdif_'  .OR.  var(1:16) == 'rtm_rad_inswref_'  .OR. &
+            var(1:16) == 'rtm_rad_inlwdif_'  .OR.  var(1:16) == 'rtm_rad_inlwref_'  .OR. &
+            var(1:14) == 'rtm_rad_outsw_'  .OR.  var(1:14) == 'rtm_rad_outlw_'  .OR.     &
+            var(1:14) == 'rtm_rad_ressw_'  .OR.  var(1:14) == 'rtm_rad_reslw_'  ) THEN
+          IF ( .NOT.  radiation ) THEN
                 message_string = 'output of "' // TRIM( var ) // '" require'&
                                  // 's radiation = .TRUE.'
                 CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
-             ENDIF
-             IF ( mrt_nlevels == 0 ) THEN
+          ENDIF
+          unit = 'W/m2'
+       ELSE IF ( var(1:7) == 'rtm_svf'  .OR.  var(1:7) == 'rtm_dif'  .OR.                &
+                 var(1:9) == 'rtm_skyvf' .OR. var(1:9) == 'rtm_skyvft' )  THEN
+          IF ( .NOT.  radiation ) THEN
                 message_string = 'output of "' // TRIM( var ) // '" require'&
-                                 // 's mrt_nlevels > 0'
-                CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
-             ENDIF
-             IF ( TRIM( var ) == 'rad_mrt_sw'  .AND.  .NOT. mrt_include_sw ) THEN
-                message_string = 'output of "' // TRIM( var ) // '" require'&
-                                 // 's rad_mrt_sw = .TRUE.'
-                CALL message( 'check_parameters', 'PA0511', 1, 2, 0, 6, 0 )
-             ENDIF
-             IF ( TRIM( var ) == 'rad_mrt' ) THEN
-                unit = 'K'
-             ELSE
-                unit = 'W m-2'
-             ENDIF
-
-          CASE DEFAULT
-             unit = 'illegal'
-
-       END SELECT
-
+                                 // 's radiation = .TRUE.'
+                CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
+          ENDIF
+          unit = '1'
+       ELSE
+!--       non-directional variables
+          SELECT CASE ( TRIM( var ) )
+             CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_lw_in', 'rad_lw_out', &
+                    'rad_sw_cs_hr', 'rad_sw_hr', 'rad_sw_in', 'rad_sw_out'  )
+                IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
+                   message_string = '"output of "' // TRIM( var ) // '" requi' // &
+                                    'res radiation = .TRUE. and ' //              &
+                                    'radiation_scheme = "rrtmg"'
+                   CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
+                ENDIF
+                unit = 'K/h'
+
+             CASE ( 'rad_net*', 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*',      &
+                    'rrtm_asdir*', 'rad_lw_in*', 'rad_lw_out*', 'rad_sw_in*',     &
+                    'rad_sw_out*')
+                IF ( i == 0 .AND. ilen == 0 .AND. k == 0)  THEN
+                   ! Workaround for masked output (calls with i=ilen=k=0)
+                   unit = 'illegal'
+                   RETURN
+                ENDIF
+                IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
+                   message_string = 'illegal value for data_output: "' //         &
+                                    TRIM( var ) // '" & only 2d-horizontal ' //   &
+                                    'cross sections are allowed for this value'
+                   CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
+                ENDIF
+                IF (  .NOT.  radiation  .OR.  radiation_scheme /= "rrtmg" )  THEN
+                   IF ( TRIM( var ) == 'rrtm_aldif*'  .OR.                        &
+                        TRIM( var ) == 'rrtm_aldir*'  .OR.                        &
+                        TRIM( var ) == 'rrtm_asdif*'  .OR.                        &
+                        TRIM( var ) == 'rrtm_asdir*'      )                       &
+                   THEN
+                      message_string = 'output of "' // TRIM( var ) // '" require'&
+                                       // 's radiation = .TRUE. and radiation_sch'&
+                                       // 'eme = "rrtmg"'
+                      CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
+                   ENDIF
+                ENDIF
+
+                IF ( TRIM( var ) == 'rad_net*'      ) unit = 'W/m2'
+                IF ( TRIM( var ) == 'rad_lw_in*'    ) unit = 'W/m2'
+                IF ( TRIM( var ) == 'rad_lw_out*'   ) unit = 'W/m2'
+                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 = ''
+                IF ( TRIM( var ) == 'rrtm_asdif*'   ) unit = ''
+                IF ( TRIM( var ) == 'rrtm_asdir*'   ) unit = ''
+
+             CASE ( 'rtm_rad_pc_inlw', 'rtm_rad_pc_insw', 'rtm_rad_pc_inswdir', &
+                    'rtm_rad_pc_inswdif', 'rtm_rad_pc_inswref')
+                IF ( .NOT.  radiation ) THEN
+                   message_string = 'output of "' // TRIM( var ) // '" require'&
+                                    // 's radiation = .TRUE.'
+                   CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
+                ENDIF
+                unit = 'W'
+
+             CASE ( 'rtm_mrt', 'rtm_mrt_sw', 'rtm_mrt_lw'  )
+                IF ( i == 0 .AND. ilen == 0 .AND. k == 0)  THEN
+                   ! Workaround for masked output (calls with i=ilen=k=0)
+                   unit = 'illegal'
+                   RETURN
+                ENDIF
+
+                IF ( .NOT.  radiation ) THEN
+                   message_string = 'output of "' // TRIM( var ) // '" require'&
+                                    // 's radiation = .TRUE.'
+                   CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
+                ENDIF
+                IF ( mrt_nlevels == 0 ) THEN
+                   message_string = 'output of "' // TRIM( var ) // '" require'&
+                                    // 's mrt_nlevels > 0'
+                   CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
+                ENDIF
+                IF ( TRIM( var ) == 'rtm_mrt_sw'  .AND.  .NOT. mrt_include_sw ) THEN
+                   message_string = 'output of "' // TRIM( var ) // '" require'&
+                                    // 's rtm_mrt_sw = .TRUE.'
+                   CALL message( 'check_parameters', 'PA0511', 1, 2, 0, 6, 0 )
+                ENDIF
+                IF ( TRIM( var ) == 'rtm_mrt' ) THEN
+                   unit = 'K'
+                ELSE
+                   unit = 'W m-2'
+                ENDIF
+
+             CASE DEFAULT
+                unit = 'illegal'
+
+          END SELECT
+       ENDIF
 
     END SUBROUTINE radiation_check_data_output
@@ -5826,4 +5882,6 @@
        endwall = nsurfl
        nwalls  = endwall - startwall + 1
+       dirstart = (/ startland, startwall, startwall, startwall, startwall /)
+       dirend = (/ endland, endwall, endwall, endwall, endwall /)
 
 !--    fill gridpcbl and pcbl
@@ -8656,8 +8714,30 @@
     INTEGER(iwp) ::  l, m !< index of current surface element
 
+    INTEGER(iwp)                                       :: ids, idsint_u, idsint_l, isurf
+    CHARACTER(LEN=varnamelength)                       :: var
+
+!-- find the real name of the variable
+    ids = -1
+    l = -1
+    var = TRIM(variable)
+    DO i = 0, nd-1
+        k = len(TRIM(var))
+        j = len(TRIM(dirname(i)))
+        IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) )  THEN
+            ids = i
+            idsint_u = dirint_u(ids)
+            idsint_l = dirint_l(ids)
+            var = var(:k-j)
+            EXIT
+        ENDIF
+    ENDDO
+    IF ( ids == -1 )  THEN
+        var = TRIM(variable)
+    ENDIF
+
     IF ( mode == 'allocate' )  THEN
 
-       SELECT CASE ( TRIM( variable ) )
-
+       SELECT CASE ( TRIM( var ) )
+!--          block of large scale (e.g. RRTMG) radiation output variables
              CASE ( 'rad_net*' )
                 IF ( .NOT. ALLOCATED( rad_net_av ) )  THEN
@@ -8738,5 +8818,124 @@
                 rad_sw_hr_av = 0.0_wp
 
-             CASE ( 'rad_mrt_sw' )
+!--          block of RTM output variables
+             CASE ( 'rtm_rad_net' )
+!--              array of complete radiation balance
+                 IF ( .NOT.  ALLOCATED(surfradnet_av) )  THEN
+                     ALLOCATE( surfradnet_av(nsurfl) )
+                     surfradnet_av = 0.0_wp
+                 ENDIF
+
+             CASE ( 'rtm_rad_insw' )
+!--                 array of sw radiation falling to surface after i-th reflection
+                 IF ( .NOT.  ALLOCATED(surfinsw_av) )  THEN
+                     ALLOCATE( surfinsw_av(nsurfl) )
+                     surfinsw_av = 0.0_wp
+                 ENDIF
+
+             CASE ( 'rtm_rad_inlw' )
+!--                 array of lw radiation falling to surface after i-th reflection
+                 IF ( .NOT.  ALLOCATED(surfinlw_av) )  THEN
+                     ALLOCATE( surfinlw_av(nsurfl) )
+                     surfinlw_av = 0.0_wp
+                 ENDIF
+
+             CASE ( 'rtm_rad_inswdir' )
+!--                 array of direct sw radiation falling to surface from sun
+                 IF ( .NOT.  ALLOCATED(surfinswdir_av) )  THEN
+                     ALLOCATE( surfinswdir_av(nsurfl) )
+                     surfinswdir_av = 0.0_wp
+                 ENDIF
+
+             CASE ( 'rtm_rad_inswdif' )
+!--                 array of difusion sw radiation falling to surface from sky and borders of the domain
+                 IF ( .NOT.  ALLOCATED(surfinswdif_av) )  THEN
+                     ALLOCATE( surfinswdif_av(nsurfl) )
+                     surfinswdif_av = 0.0_wp
+                 ENDIF
+
+             CASE ( 'rtm_rad_inswref' )
+!--                 array of sw radiation falling to surface from reflections
+                 IF ( .NOT.  ALLOCATED(surfinswref_av) )  THEN
+                     ALLOCATE( surfinswref_av(nsurfl) )
+                     surfinswref_av = 0.0_wp
+                 ENDIF
+
+             CASE ( 'rtm_rad_inlwdif' )
+!--                 array of sw radiation falling to surface after i-th reflection
+                IF ( .NOT.  ALLOCATED(surfinlwdif_av) )  THEN
+                     ALLOCATE( surfinlwdif_av(nsurfl) )
+                     surfinlwdif_av = 0.0_wp
+                 ENDIF
+
+             CASE ( 'rtm_rad_inlwref' )
+!--                 array of lw radiation falling to surface from reflections
+                 IF ( .NOT.  ALLOCATED(surfinlwref_av) )  THEN
+                     ALLOCATE( surfinlwref_av(nsurfl) )
+                     surfinlwref_av = 0.0_wp
+                 ENDIF
+
+             CASE ( 'rtm_rad_outsw' )
+!--                 array of sw radiation emitted from surface after i-th reflection
+                 IF ( .NOT.  ALLOCATED(surfoutsw_av) )  THEN
+                     ALLOCATE( surfoutsw_av(nsurfl) )
+                     surfoutsw_av = 0.0_wp
+                 ENDIF
+
+             CASE ( 'rtm_rad_outlw' )
+!--                 array of lw radiation emitted from surface after i-th reflection
+                 IF ( .NOT.  ALLOCATED(surfoutlw_av) )  THEN
+                     ALLOCATE( surfoutlw_av(nsurfl) )
+                     surfoutlw_av = 0.0_wp
+                 ENDIF
+             CASE ( 'rtm_rad_ressw' )
+!--                 array of residua of sw radiation absorbed in surface after last reflection
+                 IF ( .NOT.  ALLOCATED(surfins_av) )  THEN
+                     ALLOCATE( surfins_av(nsurfl) )
+                     surfins_av = 0.0_wp
+                 ENDIF
+
+             CASE ( 'rtm_rad_reslw' )
+!--                 array of residua of lw radiation absorbed in surface after last reflection
+                 IF ( .NOT.  ALLOCATED(surfinl_av) )  THEN
+                     ALLOCATE( surfinl_av(nsurfl) )
+                     surfinl_av = 0.0_wp
+                 ENDIF
+
+             CASE ( 'rtm_rad_pc_inlw' )
+!--                 array of of lw radiation absorbed in plant canopy
+                 IF ( .NOT.  ALLOCATED(pcbinlw_av) )  THEN
+                     ALLOCATE( pcbinlw_av(1:npcbl) )
+                     pcbinlw_av = 0.0_wp
+                 ENDIF
+
+             CASE ( 'rtm_rad_pc_insw' )
+!--                 array of of sw radiation absorbed in plant canopy
+                 IF ( .NOT.  ALLOCATED(pcbinsw_av) )  THEN
+                     ALLOCATE( pcbinsw_av(1:npcbl) )
+                     pcbinsw_av = 0.0_wp
+                 ENDIF
+
+             CASE ( 'rtm_rad_pc_inswdir' )
+!--                 array of of direct sw radiation absorbed in plant canopy
+                 IF ( .NOT.  ALLOCATED(pcbinswdir_av) )  THEN
+                     ALLOCATE( pcbinswdir_av(1:npcbl) )
+                     pcbinswdir_av = 0.0_wp
+                 ENDIF
+
+             CASE ( 'rtm_rad_pc_inswdif' )
+!--                 array of of diffuse sw radiation absorbed in plant canopy
+                 IF ( .NOT.  ALLOCATED(pcbinswdif_av) )  THEN
+                     ALLOCATE( pcbinswdif_av(1:npcbl) )
+                     pcbinswdif_av = 0.0_wp
+                 ENDIF
+
+             CASE ( 'rtm_rad_pc_inswref' )
+!--                 array of of reflected sw radiation absorbed in plant canopy
+                 IF ( .NOT.  ALLOCATED(pcbinswref_av) )  THEN
+                     ALLOCATE( pcbinswref_av(1:npcbl) )
+                     pcbinswref_av = 0.0_wp
+                 ENDIF
+
+             CASE ( 'rtm_mrt_sw' )
                 IF ( .NOT. ALLOCATED( mrtinsw_av ) )  THEN
                    ALLOCATE( mrtinsw_av(nmrtbl) )
@@ -8744,5 +8943,5 @@
                 mrtinsw_av = 0.0_wp
 
-             CASE ( 'rad_mrt_lw' )
+             CASE ( 'rtm_mrt_lw' )
                 IF ( .NOT. ALLOCATED( mrtinlw_av ) )  THEN
                    ALLOCATE( mrtinlw_av(nmrtbl) )
@@ -8750,5 +8949,5 @@
                 mrtinlw_av = 0.0_wp
 
-             CASE ( 'rad_mrt' )
+             CASE ( 'rtm_mrt' )
                 IF ( .NOT. ALLOCATED( mrt_av ) )  THEN
                    ALLOCATE( mrt_av(nmrtbl) )
@@ -8763,6 +8962,6 @@
     ELSEIF ( mode == 'sum' )  THEN
 
-       SELECT CASE ( TRIM( variable ) )
-
+       SELECT CASE ( TRIM( var ) )
+!--       block of large scale (e.g. RRTMG) radiation output variables
           CASE ( 'rad_net*' )
              IF ( ALLOCATED( rad_net_av ) ) THEN
@@ -8971,4 +9170,129 @@
              ENDIF
 
+!--       block of RTM output variables
+          CASE ( 'rtm_rad_net' )
+!--           array of complete radiation balance
+              DO isurf = dirstart(ids), dirend(ids)
+                 IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                    surfradnet_av(isurf) = surfinsw(isurf) - surfoutsw(isurf) +  surfinlw(isurf) - surfoutlw(isurf)
+                 ENDIF
+              ENDDO
+
+          CASE ( 'rtm_rad_insw' )
+!--           array of sw radiation falling to surface after i-th reflection
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfinsw_av(isurf) = surfinsw_av(isurf) + surfinsw(isurf)
+                  ENDIF
+              ENDDO
+
+          CASE ( 'rtm_rad_inlw' )
+!--           array of lw radiation falling to surface after i-th reflection
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfinlw_av(isurf) = surfinlw_av(isurf) + surfinlw(isurf)
+                  ENDIF
+              ENDDO
+
+          CASE ( 'rtm_rad_inswdir' )
+!--           array of direct sw radiation falling to surface from sun
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfinswdir_av(isurf) = surfinswdir_av(isurf) + surfinswdir(isurf)
+                  ENDIF
+              ENDDO
+
+          CASE ( 'rtm_rad_inswdif' )
+!--           array of difusion sw radiation falling to surface from sky and borders of the domain
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfinswdif_av(isurf) = surfinswdif_av(isurf) + surfinswdif(isurf)
+                  ENDIF
+              ENDDO
+
+          CASE ( 'rtm_rad_inswref' )
+!--           array of sw radiation falling to surface from reflections
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfinswref_av(isurf) = surfinswref_av(isurf) + surfinsw(isurf) - &
+                                          surfinswdir(isurf) - surfinswdif(isurf)
+                  ENDIF
+              ENDDO
+
+
+          CASE ( 'rtm_rad_inlwdif' )
+!--           array of sw radiation falling to surface after i-th reflection
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfinlwdif_av(isurf) = surfinlwdif_av(isurf) + surfinlwdif(isurf)
+                  ENDIF
+              ENDDO
+!
+          CASE ( 'rtm_rad_inlwref' )
+!--           array of lw radiation falling to surface from reflections
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfinlwref_av(isurf) = surfinlwref_av(isurf) + &
+                                          surfinlw(isurf) - surfinlwdif(isurf)
+                  ENDIF
+              ENDDO
+
+          CASE ( 'rtm_rad_outsw' )
+!--           array of sw radiation emitted from surface after i-th reflection
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfoutsw_av(isurf) = surfoutsw_av(isurf) + surfoutsw(isurf)
+                  ENDIF
+              ENDDO
+
+          CASE ( 'rtm_rad_outlw' )
+!--           array of lw radiation emitted from surface after i-th reflection
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfoutlw_av(isurf) = surfoutlw_av(isurf) + surfoutlw(isurf)
+                  ENDIF
+              ENDDO
+
+          CASE ( 'rtm_rad_ressw' )
+!--           array of residua of sw radiation absorbed in surface after last reflection
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfins_av(isurf) = surfins_av(isurf) + surfins(isurf)
+                  ENDIF
+              ENDDO
+
+          CASE ( 'rtm_rad_reslw' )
+!--           array of residua of lw radiation absorbed in surface after last reflection
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfinl_av(isurf) = surfinl_av(isurf) + surfinl(isurf)
+                  ENDIF
+              ENDDO
+
+          CASE ( 'rtm_rad_pc_inlw' )
+              DO l = 1, npcbl
+                 pcbinlw_av(l) = pcbinlw_av(l) + pcbinlw(l)
+              ENDDO
+
+          CASE ( 'rtm_rad_pc_insw' )
+              DO l = 1, npcbl
+                 pcbinsw_av(l) = pcbinsw_av(l) + pcbinsw(l)
+              ENDDO
+
+          CASE ( 'rtm_rad_pc_inswdir' )
+              DO l = 1, npcbl
+                 pcbinswdir_av(l) = pcbinswdir_av(l) + pcbinswdir(l)
+              ENDDO
+
+          CASE ( 'rtm_rad_pc_inswdif' )
+              DO l = 1, npcbl
+                 pcbinswdif_av(l) = pcbinswdif_av(l) + pcbinswdif(l)
+              ENDDO
+
+          CASE ( 'rtm_rad_pc_inswref' )
+              DO l = 1, npcbl
+                 pcbinswref_av(l) = pcbinswref_av(l) + pcbinsw(l) - pcbinswdir(l) - pcbinswdif(l)
+              ENDDO
+
           CASE ( 'rad_mrt_sw' )
              IF ( ALLOCATED( mrtinsw_av ) )  THEN
@@ -8993,6 +9317,6 @@
     ELSEIF ( mode == 'average' )  THEN
 
-       SELECT CASE ( TRIM( variable ) )
-
+       SELECT CASE ( TRIM( var ) )
+!--       block of large scale (e.g. RRTMG) radiation output variables
           CASE ( 'rad_net*' )
              IF ( ALLOCATED( rad_net_av ) ) THEN
@@ -9141,8 +9465,125 @@
              ENDIF
 
-          CASE ( 'rad_mrt_sw' )
-             IF ( ALLOCATED( mrtinsw_av ) )  THEN
-                mrtinsw_av(:) = mrtinsw_av(:)  / REAL( average_count_3d, KIND=wp )
-             ENDIF
+!--       block of RTM output variables
+          CASE ( 'rtm_rad_net' )
+!--           array of complete radiation balance
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfradnet_av(isurf) = surfinsw_av(isurf) / REAL( average_count_3d, kind=wp )
+                  ENDIF
+              ENDDO
+
+          CASE ( 'rtm_rad_insw' )
+!--           array of sw radiation falling to surface after i-th reflection
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfinsw_av(isurf) = surfinsw_av(isurf) / REAL( average_count_3d, kind=wp )
+                  ENDIF
+              ENDDO
+
+          CASE ( 'rtm_rad_inlw' )
+!--           array of lw radiation falling to surface after i-th reflection
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfinlw_av(isurf) = surfinlw_av(isurf) / REAL( average_count_3d, kind=wp )
+                  ENDIF
+              ENDDO
+
+          CASE ( 'rtm_rad_inswdir' )
+!--           array of direct sw radiation falling to surface from sun
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfinswdir_av(isurf) = surfinswdir_av(isurf) / REAL( average_count_3d, kind=wp )
+                  ENDIF
+              ENDDO
+
+          CASE ( 'rtm_rad_inswdif' )
+!--           array of difusion sw radiation falling to surface from sky and borders of the domain
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfinswdif_av(isurf) = surfinswdif_av(isurf) / REAL( average_count_3d, kind=wp )
+                  ENDIF
+              ENDDO
+
+          CASE ( 'rtm_rad_inswref' )
+!--           array of sw radiation falling to surface from reflections
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfinswref_av(isurf) = surfinswref_av(isurf) / REAL( average_count_3d, kind=wp )
+                  ENDIF
+              ENDDO
+
+          CASE ( 'rtm_rad_inlwdif' )
+!--           array of sw radiation falling to surface after i-th reflection
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfinlwdif_av(isurf) = surfinlwdif_av(isurf) / REAL( average_count_3d, kind=wp )
+                  ENDIF
+              ENDDO
+
+          CASE ( 'rtm_rad_inlwref' )
+!--           array of lw radiation falling to surface from reflections
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfinlwref_av(isurf) = surfinlwref_av(isurf) / REAL( average_count_3d, kind=wp )
+                  ENDIF
+              ENDDO
+
+          CASE ( 'rtm_rad_outsw' )
+!--           array of sw radiation emitted from surface after i-th reflection
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfoutsw_av(isurf) = surfoutsw_av(isurf) / REAL( average_count_3d, kind=wp )
+                  ENDIF
+              ENDDO
+
+          CASE ( 'rtm_rad_outlw' )
+!--           array of lw radiation emitted from surface after i-th reflection
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfoutlw_av(isurf) = surfoutlw_av(isurf) / REAL( average_count_3d, kind=wp )
+                  ENDIF
+              ENDDO
+
+          CASE ( 'rtm_rad_ressw' )
+!--           array of residua of sw radiation absorbed in surface after last reflection
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfins_av(isurf) = surfins_av(isurf) / REAL( average_count_3d, kind=wp )
+                  ENDIF
+              ENDDO
+
+          CASE ( 'rtm_rad_reslw' )
+!--           array of residua of lw radiation absorbed in surface after last reflection
+              DO isurf = dirstart(ids), dirend(ids)
+                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+                      surfinl_av(isurf) = surfinl_av(isurf) / REAL( average_count_3d, kind=wp )
+                  ENDIF
+              ENDDO
+
+          CASE ( 'rtm_rad_pc_inlw' )
+              DO l = 1, npcbl
+                 pcbinlw_av(:) = pcbinlw_av(:) / REAL( average_count_3d, kind=wp )
+              ENDDO
+
+          CASE ( 'rtm_rad_pc_insw' )
+              DO l = 1, npcbl
+                 pcbinsw_av(:) = pcbinsw_av(:) / REAL( average_count_3d, kind=wp )
+              ENDDO
+
+          CASE ( 'rtm_rad_pc_inswdir' )
+              DO l = 1, npcbl
+                 pcbinswdir_av(:) = pcbinswdir_av(:) / REAL( average_count_3d, kind=wp )
+              ENDDO
+
+          CASE ( 'rtm_rad_pc_inswdif' )
+              DO l = 1, npcbl
+                 pcbinswdif_av(:) = pcbinswdif_av(:) / REAL( average_count_3d, kind=wp )
+              ENDDO
+
+          CASE ( 'rtm_rad_pc_inswref' )
+              DO l = 1, npcbl
+                 pcbinswref_av(:) = pcbinswref_av(:) / REAL( average_count_3d, kind=wp )
+              ENDDO
 
           CASE ( 'rad_mrt_lw' )
@@ -9170,9 +9611,9 @@
 !> It is called out from subroutine netcdf.
 !------------------------------------------------------------------------------!
-SUBROUTINE radiation_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
+SUBROUTINE radiation_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
     
     IMPLICIT NONE
 
-    CHARACTER (LEN=*), INTENT(IN)  ::  var         !< 
+    CHARACTER (LEN=*), INTENT(IN)  ::  variable    !<
     LOGICAL, INTENT(OUT)           ::  found       !< 
     CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !< 
@@ -9180,36 +9621,60 @@
     CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !< 
 
+    CHARACTER (len=varnamelength)  :: var
+
     found  = .TRUE.
 
 !
 !-- Check for the grid
-    SELECT CASE ( TRIM( var ) )
-
-       CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr',        &
-              'rad_lw_cs_hr_xy', 'rad_lw_hr_xy', 'rad_sw_cs_hr_xy',            &
-              'rad_sw_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_hr_xz',               &
-              'rad_sw_cs_hr_xz', 'rad_sw_hr_xz', 'rad_lw_cs_hr_yz',            &
-              'rad_lw_hr_yz', 'rad_sw_cs_hr_yz', 'rad_sw_hr_yz',               &
-              'rad_mrt', 'rad_mrt_sw', 'rad_mrt_lw' )
-          grid_x = 'x'
-          grid_y = 'y'
-          grid_z = 'zu'
-
-       CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out',            &
-              'rad_lw_in_xy', 'rad_lw_out_xy', 'rad_sw_in_xy','rad_sw_out_xy', &
-              'rad_lw_in_xz', 'rad_lw_out_xz', 'rad_sw_in_xz','rad_sw_out_xz', &
-              'rad_lw_in_yz', 'rad_lw_out_yz', 'rad_sw_in_yz','rad_sw_out_yz' )
-          grid_x = 'x'
-          grid_y = 'y'
-          grid_z = 'zw'
-
-
-       CASE DEFAULT
-          found  = .FALSE.
-          grid_x = 'none'
-          grid_y = 'none'
-          grid_z = 'none'
-
-        END SELECT
+    var = TRIM(variable)
+!-- RTM directional variables
+    IF ( var(1:12) == 'rtm_rad_net_'  .OR.  var(1:13) == 'rtm_rad_insw_'  .OR.          &
+         var(1:13) == 'rtm_rad_inlw_'  .OR.  var(1:16) == 'rtm_rad_inswdir_'  .OR.      &
+         var(1:16) == 'rtm_rad_inswdif_'  .OR.  var(1:16) == 'rtm_rad_inswref_'  .OR.   &
+         var(1:16) == 'rtm_rad_inlwdif_'  .OR.  var(1:16) == 'rtm_rad_inlwref_'  .OR.   &
+         var(1:14) == 'rtm_rad_outsw_'  .OR.  var(1:14) == 'rtm_rad_outlw_'  .OR.       &
+         var(1:14) == 'rtm_rad_ressw_'  .OR.  var(1:14) == 'rtm_rad_reslw_'  .OR.       &
+         var == 'rtm_rad_pc_inlw'  .OR.                                                 &
+         var == 'rtm_rad_pc_insw'  .OR.  var == 'rtm_rad_pc_inswdir'  .OR.              &
+         var == 'rtm_rad_pc_inswdif'  .OR.  var == 'rtm_rad_pc_inswref'  .OR.           &
+         var(1:7) == 'rtm_svf'  .OR.  var(1:7) == 'rtm_dif'  .OR.                       &
+         var(1:9) == 'rtm_skyvf' .OR. var(1:10) == 'rtm_skyvft'  .OR.                   &
+         var == 'rtm_mrt'  .OR.  var ==  'rtm_mrt_sw'  .OR.  var == 'rtm_mrt_lw' )  THEN
+
+         found = .TRUE.
+         grid_x = 'x'
+         grid_y = 'y'
+         grid_z = 'zu'
+    ELSE
+
+       SELECT CASE ( TRIM( var ) )
+
+          CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr',        &
+                 'rad_lw_cs_hr_xy', 'rad_lw_hr_xy', 'rad_sw_cs_hr_xy',            &
+                 'rad_sw_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_hr_xz',               &
+                 'rad_sw_cs_hr_xz', 'rad_sw_hr_xz', 'rad_lw_cs_hr_yz',            &
+                 'rad_lw_hr_yz', 'rad_sw_cs_hr_yz', 'rad_sw_hr_yz',               &
+                 'rad_mrt', 'rad_mrt_sw', 'rad_mrt_lw' )
+             grid_x = 'x'
+             grid_y = 'y'
+             grid_z = 'zu'
+
+          CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out',            &
+                 'rad_lw_in_xy', 'rad_lw_out_xy', 'rad_sw_in_xy','rad_sw_out_xy', &
+                 'rad_lw_in_xz', 'rad_lw_out_xz', 'rad_sw_in_xz','rad_sw_out_xz', &
+                 'rad_lw_in_yz', 'rad_lw_out_yz', 'rad_sw_in_yz','rad_sw_out_yz' )
+             grid_x = 'x'
+             grid_y = 'y'
+             grid_z = 'zw'
+
+
+          CASE DEFAULT
+             found  = .FALSE.
+             grid_x = 'none'
+             grid_y = 'none'
+             grid_z = 'none'
+
+           END SELECT
+       ENDIF
 
     END SUBROUTINE radiation_define_netcdf_grid
@@ -9649,9 +10114,48 @@
     REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
 
+    CHARACTER (len=varnamelength)                   :: var, surfid
+    INTEGER(iwp)                                    :: ids,idsint_u,idsint_l,isurf,isvf,isurfs,isurflt,ipcgb
+    INTEGER(iwp)                                    :: is, js, ks, istat
+
     found = .TRUE.
 
-
-    SELECT CASE ( TRIM( variable ) )
-
+    ids = -1
+    var = TRIM(variable)
+    DO i = 0, nd-1
+        k = len(TRIM(var))
+        j = len(TRIM(dirname(i)))
+        IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) )  THEN
+            ids = i
+            idsint_u = dirint_u(ids)
+            idsint_l = dirint_l(ids)
+            var = var(:k-j)
+            EXIT
+        ENDIF
+    ENDDO
+    IF ( ids == -1 )  THEN
+        var = TRIM(variable)
+    ENDIF
+
+    IF ( (var(1:8) == 'rtm_svf_'  .OR.  var(1:8) == 'rtm_dif_')  .AND.  len(TRIM(var)) >= 13 )  THEN
+!--     svf values to particular surface
+        surfid = var(9:)
+        i = index(surfid,'_')
+        j = index(surfid(i+1:),'_')
+        READ(surfid(1:i-1),*, iostat=istat ) is
+        IF ( istat == 0 )  THEN
+            READ(surfid(i+1:i+j-1),*, iostat=istat ) js
+        ENDIF
+        IF ( istat == 0 )  THEN
+            READ(surfid(i+j+1:),*, iostat=istat ) ks
+        ENDIF
+        IF ( istat == 0 )  THEN
+            var = var(1:7)
+        ENDIF
+    ENDIF
+
+    local_pf = fill_value
+
+    SELECT CASE ( TRIM( var ) )
+!--   block of large scale radiation model (e.g. RRTMG) output variables
       CASE ( 'rad_sw_in' )
          IF ( av == 0 )  THEN
@@ -9838,60 +10342,275 @@
          ENDIF
 
-      CASE ( 'rad_mrt_sw' )
+!--   block of RTM output variables
+!--   variables are intended mainly for debugging and detailed analyse purposes
+      CASE ( 'rtm_skyvf' )
+!--        sky view factor
+         DO isurf = dirstart(ids), dirend(ids)
+            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+               local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = skyvf(isurf)
+            ENDIF
+         ENDDO
+
+      CASE ( 'rtm_skyvft' )
+!--      sky view factor
+         DO isurf = dirstart(ids), dirend(ids)
+            IF ( surfl(id,isurf) == ids )  THEN
+               local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = skyvft(isurf)
+            ENDIF
+         ENDDO
+
+      CASE ( 'rtm_svf', 'rtm_dif' )
+!--      shape view factors or iradiance factors to selected surface
+         IF ( TRIM(var)=='rtm_svf' )  THEN
+             k = 1
+         ELSE
+             k = 2
+         ENDIF
+         DO isvf = 1, nsvfl
+            isurflt = svfsurf(1, isvf)
+            isurfs = svfsurf(2, isvf)
+
+            IF ( surf(ix,isurfs) == is  .AND.  surf(iy,isurfs) == js  .AND. surf(iz,isurfs) == ks  .AND. &
+                 (surf(id,isurfs) == idsint_u .OR. surfl(id,isurf) == idsint_l ) ) THEN
+!--            correct source surface
+               local_pf(surfl(ix,isurflt),surfl(iy,isurflt),surfl(iz,isurflt)) = svf(k,isvf)
+            ENDIF
+         ENDDO
+
+      CASE ( 'rtm_rad_net' )
+!--     array of complete radiation balance
+         DO isurf = dirstart(ids), dirend(ids)
+            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+               IF ( av == 0 )  THEN
+                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = &
+                         surfinsw(isurf) - surfoutsw(isurf) +  surfinlw(isurf) - surfoutlw(isurf)
+               ELSE
+                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfradnet_av(isurf)
+               ENDIF
+            ENDIF
+         ENDDO
+
+      CASE ( 'rtm_rad_insw' )
+!--      array of sw radiation falling to surface after i-th reflection
+         DO isurf = dirstart(ids), dirend(ids)
+            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+               IF ( av == 0 )  THEN
+                 local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinsw(isurf)
+               ELSE
+                 local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinsw_av(isurf)
+               ENDIF
+            ENDIF
+         ENDDO
+
+      CASE ( 'rtm_rad_inlw' )
+!--      array of lw radiation falling to surface after i-th reflection
+         DO isurf = dirstart(ids), dirend(ids)
+            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+               IF ( av == 0 )  THEN
+                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinlw(isurf)
+               ELSE
+                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinlw_av(isurf)
+               ENDIF
+             ENDIF
+         ENDDO
+
+      CASE ( 'rtm_rad_inswdir' )
+!--      array of direct sw radiation falling to surface from sun
+         DO isurf = dirstart(ids), dirend(ids)
+            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+               IF ( av == 0 )  THEN
+                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinswdir(isurf)
+               ELSE
+                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinswdir_av(isurf)
+               ENDIF
+            ENDIF
+         ENDDO
+
+      CASE ( 'rtm_rad_inswdif' )
+!--      array of difusion sw radiation falling to surface from sky and borders of the domain
+         DO isurf = dirstart(ids), dirend(ids)
+            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+               IF ( av == 0 )  THEN
+                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinswdif(isurf)
+               ELSE
+                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinswdif_av(isurf)
+               ENDIF
+            ENDIF
+         ENDDO
+
+      CASE ( 'rtm_rad_inswref' )
+!--      array of sw radiation falling to surface from reflections
+         DO isurf = dirstart(ids), dirend(ids)
+            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+               IF ( av == 0 )  THEN
+                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = &
+                    surfinsw(isurf) - surfinswdir(isurf) - surfinswdif(isurf)
+               ELSE
+                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinswref_av(isurf)
+               ENDIF
+            ENDIF
+         ENDDO
+
+      CASE ( 'rtm_rad_inlwdif' )
+!--      array of difusion lw radiation falling to surface from sky and borders of the domain
+         DO isurf = dirstart(ids), dirend(ids)
+            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+               IF ( av == 0 )  THEN
+                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinlwdif(isurf)
+               ELSE
+                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinlwdif_av(isurf)
+               ENDIF
+            ENDIF
+         ENDDO
+
+      CASE ( 'rtm_rad_inlwref' )
+!--      array of lw radiation falling to surface from reflections
+         DO isurf = dirstart(ids), dirend(ids)
+            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+               IF ( av == 0 )  THEN
+                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinlw(isurf) - surfinlwdif(isurf)
+               ELSE
+                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinlwref_av(isurf)
+               ENDIF
+            ENDIF
+         ENDDO
+
+      CASE ( 'rtm_rad_outsw' )
+!--      array of sw radiation emitted from surface after i-th reflection
+         DO isurf = dirstart(ids), dirend(ids)
+            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+               IF ( av == 0 )  THEN
+                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfoutsw(isurf)
+               ELSE
+                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfoutsw_av(isurf)
+               ENDIF
+            ENDIF
+         ENDDO
+
+      CASE ( 'rtm_rad_outlw' )
+!--      array of lw radiation emitted from surface after i-th reflection
+         DO isurf = dirstart(ids), dirend(ids)
+            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+               IF ( av == 0 )  THEN
+                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfoutlw(isurf)
+               ELSE
+                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfoutlw_av(isurf)
+               ENDIF
+            ENDIF
+         ENDDO
+
+      CASE ( 'rtm_rad_ressw' )
+!--      average of array of residua of sw radiation absorbed in surface after last reflection
+         DO isurf = dirstart(ids), dirend(ids)
+            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+               IF ( av == 0 )  THEN
+                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfins(isurf)
+               ELSE
+                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfins_av(isurf)
+               ENDIF
+            ENDIF
+         ENDDO
+
+      CASE ( 'rtm_rad_reslw' )
+!--      average of array of residua of lw radiation absorbed in surface after last reflection
+         DO isurf = dirstart(ids), dirend(ids)
+            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+               IF ( av == 0 )  THEN
+                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinl(isurf)
+               ELSE
+                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinl_av(isurf)
+               ENDIF
+            ENDIF
+         ENDDO
+
+      CASE ( 'rtm_rad_pc_inlw' )
+!--      array of lw radiation absorbed by plant canopy
+         DO ipcgb = 1, npcbl
+            IF ( av == 0 )  THEN
+               local_pf(pcbl(ix,ipcgb),pcbl(iy,ipcgb),pcbl(iz,ipcgb)) = pcbinlw(ipcgb)
+            ELSE
+               local_pf(pcbl(ix,ipcgb),pcbl(iy,ipcgb),pcbl(iz,ipcgb)) = pcbinlw_av(ipcgb)
+            ENDIF
+         ENDDO
+
+      CASE ( 'rtm_rad_pc_insw' )
+!--      array of sw radiation absorbed by plant canopy
+         DO ipcgb = 1, npcbl
+            IF ( av == 0 )  THEN
+              local_pf(pcbl(ix,ipcgb),pcbl(iy,ipcgb),pcbl(iz,ipcgb)) = pcbinsw(ipcgb)
+            ELSE
+              local_pf(pcbl(ix,ipcgb),pcbl(iy,ipcgb),pcbl(iz,ipcgb)) = pcbinsw_av(ipcgb)
+            ENDIF
+         ENDDO
+
+      CASE ( 'rtm_rad_pc_inswdir' )
+!--      array of direct sw radiation absorbed by plant canopy
+         DO ipcgb = 1, npcbl
+            IF ( av == 0 )  THEN
+               local_pf(pcbl(ix,ipcgb),pcbl(iy,ipcgb),pcbl(iz,ipcgb)) = pcbinswdir(ipcgb)
+            ELSE
+               local_pf(pcbl(ix,ipcgb),pcbl(iy,ipcgb),pcbl(iz,ipcgb)) = pcbinswdir_av(ipcgb)
+            ENDIF
+         ENDDO
+
+      CASE ( 'rtm_rad_pc_inswdif' )
+!--      array of diffuse sw radiation absorbed by plant canopy
+         DO ipcgb = 1, npcbl
+            IF ( av == 0 )  THEN
+               local_pf(pcbl(ix,ipcgb),pcbl(iy,ipcgb),pcbl(iz,ipcgb)) = pcbinswdif(ipcgb)
+            ELSE
+               local_pf(pcbl(ix,ipcgb),pcbl(iy,ipcgb),pcbl(iz,ipcgb)) = pcbinswdif_av(ipcgb)
+            ENDIF
+         ENDDO
+
+      CASE ( 'rtm_rad_pc_inswref' )
+!--      array of reflected sw radiation absorbed by plant canopy
+         DO ipcgb = 1, npcbl
+            IF ( av == 0 )  THEN
+               local_pf(pcbl(ix,ipcgb),pcbl(iy,ipcgb),pcbl(iz,ipcgb)) = &
+                                    pcbinsw(ipcgb) - pcbinswdir(ipcgb) - pcbinswdif(ipcgb)
+            ELSE
+               local_pf(pcbl(ix,ipcgb),pcbl(iy,ipcgb),pcbl(iz,ipcgb)) = pcbinswref_av(ipcgb)
+            ENDIF
+         ENDDO
+
+      CASE ( 'rtm_mrt_sw' )
          local_pf = REAL( fill_value, KIND = wp )
          IF ( av == 0 )  THEN
             DO  l = 1, nmrtbl
-               i = mrtbl(ix,l)
-               j = mrtbl(iy,l)
-               k = mrtbl(iz,l)
-               local_pf(i,j,k) = mrtinsw(l)
+               local_pf(mrtbl(ix,l),mrtbl(iy,l),mrtbl(iz,l)) = mrtinsw(l)
             ENDDO
          ELSE
             IF ( ALLOCATED( mrtinsw_av ) ) THEN
                DO  l = 1, nmrtbl
-                  i = mrtbl(ix,l)
-                  j = mrtbl(iy,l)
-                  k = mrtbl(iz,l)
-                  local_pf(i,j,k) = mrtinsw_av(l)
+                  local_pf(mrtbl(ix,l),mrtbl(iy,l),mrtbl(iz,l)) = mrtinsw_av(l)
                ENDDO
             ENDIF
          ENDIF
 
-      CASE ( 'rad_mrt_lw' )
+      CASE ( 'rtm_mrt_lw' )
          local_pf = REAL( fill_value, KIND = wp )
          IF ( av == 0 )  THEN
             DO  l = 1, nmrtbl
-               i = mrtbl(ix,l)
-               j = mrtbl(iy,l)
-               k = mrtbl(iz,l)
-               local_pf(i,j,k) = mrtinlw(l)
+               local_pf(mrtbl(ix,l),mrtbl(iy,l),mrtbl(iz,l)) = mrtinlw(l)
             ENDDO
          ELSE
             IF ( ALLOCATED( mrtinlw_av ) ) THEN
                DO  l = 1, nmrtbl
-                  i = mrtbl(ix,l)
-                  j = mrtbl(iy,l)
-                  k = mrtbl(iz,l)
-                  local_pf(i,j,k) = mrtinlw_av(l)
+                  local_pf(mrtbl(ix,l),mrtbl(iy,l),mrtbl(iz,l)) = mrtinlw_av(l)
                ENDDO
             ENDIF
          ENDIF
 
-      CASE ( 'rad_mrt' )
+      CASE ( 'rtm_mrt' )
          local_pf = REAL( fill_value, KIND = wp )
          IF ( av == 0 )  THEN
             DO  l = 1, nmrtbl
-               i = mrtbl(ix,l)
-               j = mrtbl(iy,l)
-               k = mrtbl(iz,l)
-               local_pf(i,j,k) = mrt(l)
+               local_pf(mrtbl(ix,l),mrtbl(iy,l),mrtbl(iz,l)) = mrt(l)
             ENDDO
          ELSE
             IF ( ALLOCATED( mrt_av ) ) THEN
                DO  l = 1, nmrtbl
-                  i = mrtbl(ix,l)
-                  j = mrtbl(iy,l)
-                  k = mrtbl(iz,l)
-                  local_pf(i,j,k) = mrt_av(l)
+                  local_pf(mrtbl(ix,l),mrtbl(iy,l),mrtbl(iz,l)) = mrt_av(l)
                ENDDO
             ENDIF
Index: palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3606)
+++ palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3607)
@@ -28,4 +28,7 @@
 ! -----------------
 ! $Id$
+! Output of radiation-related quantities migrated to radiation_model_mod.
+! 
+! 3597 2018-12-04 08:40:18Z maronga
 ! Fixed calculation method of near surface air potential temperature at 10 cm
 ! and moved to surface_layer_fluxes. Removed unnecessary _eb strings.
@@ -438,15 +441,9 @@
                radiation, rad_sw_in, rad_lw_in, rad_sw_out, rad_lw_out,        &
                sigma_sb, sun_direction, sun_dir_lat, sun_dir_lon,              &
-               force_radiation_call, surfinsw, surfinlw, surfinswdir,          &
-               surfinswdif, surfoutsw, surfoutlw, surfins,nsvfl, svf, svfsurf, &
-               surfinl, surfinlwdif, rad_sw_in_dir, rad_sw_in_diff,            &
-               rad_lw_in_diff, surfouts, surfoutl, surfoutsl, surfoutll, surf, &
-               surfl, nsurfl, pcbinsw, pcbinlw, pcbinswdir,                    &
-               pcbinswdif, iup_u, inorth_u, isouth_u, ieast_u, iwest_u, iup_l, &
-               inorth_l, isouth_l, ieast_l, iwest_l, id,                       &
+               force_radiation_call, iup_u, inorth_u, isouth_u, ieast_u,       &
+               iwest_u, iup_l, inorth_l, isouth_l, ieast_l, iwest_l, id,       &
                iz, iy, ix,  nsurf, idsvf, ndsvf,                               &
                idcsf, ndcsf, kdcsf, pct,                                       &
-               startland, endland, startwall, endwall, skyvf, skyvft, nzub,    &
-               nzut, npcbl, pcbl, unscheduled_radiation_calls
+               nzub, nzut, unscheduled_radiation_calls
 
     USE statistics,                                                            &
@@ -461,11 +458,11 @@
     IMPLICIT NONE
 
-    !
+!
 !-- USM model constants
 
     REAL(wp), PARAMETER ::                     &
               b_ch               = 6.04_wp,    & ! Clapp & Hornberger exponent
-              lambda_h_green_dry       = 0.19_wp,    & ! heat conductivity for dry soil   
-              lambda_h_green_sm        = 3.44_wp,    & ! heat conductivity of the soil matrix
+              lambda_h_green_dry = 0.19_wp,    & ! heat conductivity for dry soil   
+              lambda_h_green_sm  = 3.44_wp,    & ! heat conductivity of the soil matrix
               lambda_h_water     = 0.57_wp,    & ! heat conductivity of water
               psi_sat            = -0.388_wp,  & ! soil matrix potential at saturation
@@ -965,21 +962,4 @@
 !-- arrays for time averages
 !-- Attention: the variable rad_net_av is also used in the 3d field variable in radiation_model_mod.f90. It may be better to rename it 
-    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinsw_av      !< average of sw radiation falling to local surface including radiation from reflections
-    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlw_av      !< average of lw radiation falling to local surface including radiation from reflections
-    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdir_av   !< average of direct sw radiation falling to local surface
-    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdif_av   !< average of diffuse sw radiation from sky and model boundary falling to local surface
-    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwdif_av   !< average of diffuse lw radiation from sky and model boundary falling to local surface
-    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswref_av   !< average of sw radiation falling to surface from reflections
-    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwref_av   !< average of lw radiation falling to surface from reflections
-    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutsw_av     !< average of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
-    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutlw_av     !< average of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
-    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins_av       !< average of array of residua of sw radiation absorbed in surface after last reflection
-    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinl_av       !< average of array of residua of lw radiation absorbed in surface after last reflection
-    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinlw_av       !< Average of pcbinlw
-    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinsw_av       !< Average of pcbinsw
-    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdir_av    !< Average of pcbinswdir
-    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdif_av    !< Average of pcbinswdif
-    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswref_av    !< Average of pcbinswref
-    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfhf_av        !< average of total radiation flux incoming to minus outgoing from local surface  
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  wghf_eb_av       !< average of wghf_eb
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  wshf_eb_av       !< average of wshf_eb
@@ -1833,143 +1813,4 @@
            
            SELECT CASE ( TRIM( var ) )
-                
-                CASE ( 'usm_rad_net' )
-!--                 array of complete radiation balance
-                    IF ( l == -1 ) THEN
-                        IF ( .NOT.  ALLOCATED(surf_usm_h%rad_net_av) )  THEN
-                            ALLOCATE( surf_usm_h%rad_net_av(1:surf_usm_h%ns) )
-                            surf_usm_h%rad_net_av = 0.0_wp
-                        ENDIF
-                    ELSE
-                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%rad_net_av) )  THEN
-                            ALLOCATE( surf_usm_v(l)%rad_net_av(1:surf_usm_v(l)%ns) )
-                            surf_usm_v(l)%rad_net_av = 0.0_wp
-                        ENDIF
-                    ENDIF
-                    
-                CASE ( 'usm_rad_insw' )
-!--                 array of sw radiation falling to surface after i-th reflection
-                    IF ( .NOT.  ALLOCATED(surfinsw_av) )  THEN
-                        ALLOCATE( surfinsw_av(nsurfl) )
-                        surfinsw_av = 0.0_wp
-                    ENDIF
-
-                CASE ( 'usm_rad_inlw' )
-!--                 array of lw radiation falling to surface after i-th reflection
-                    IF ( .NOT.  ALLOCATED(surfinlw_av) )  THEN
-                        ALLOCATE( surfinlw_av(nsurfl) )
-                        surfinlw_av = 0.0_wp
-                    ENDIF
-
-                CASE ( 'usm_rad_inswdir' )
-!--                 array of direct sw radiation falling to surface from sun
-                    IF ( .NOT.  ALLOCATED(surfinswdir_av) )  THEN
-                        ALLOCATE( surfinswdir_av(nsurfl) )
-                        surfinswdir_av = 0.0_wp
-                    ENDIF
-
-                CASE ( 'usm_rad_inswdif' )
-!--                 array of difusion sw radiation falling to surface from sky and borders of the domain
-                    IF ( .NOT.  ALLOCATED(surfinswdif_av) )  THEN
-                        ALLOCATE( surfinswdif_av(nsurfl) )
-                        surfinswdif_av = 0.0_wp
-                    ENDIF
-
-                CASE ( 'usm_rad_inswref' )
-!--                 array of sw radiation falling to surface from reflections
-                    IF ( .NOT.  ALLOCATED(surfinswref_av) )  THEN
-                        ALLOCATE( surfinswref_av(nsurfl) )
-                        surfinswref_av = 0.0_wp
-                    ENDIF
-
-                CASE ( 'usm_rad_inlwdif' )
-!--                 array of sw radiation falling to surface after i-th reflection
-                   IF ( .NOT.  ALLOCATED(surfinlwdif_av) )  THEN
-                        ALLOCATE( surfinlwdif_av(nsurfl) )
-                        surfinlwdif_av = 0.0_wp
-                    ENDIF
-
-                CASE ( 'usm_rad_inlwref' )
-!--                 array of lw radiation falling to surface from reflections
-                    IF ( .NOT.  ALLOCATED(surfinlwref_av) )  THEN
-                        ALLOCATE( surfinlwref_av(nsurfl) )
-                        surfinlwref_av = 0.0_wp
-                    ENDIF
-
-                CASE ( 'usm_rad_outsw' )
-!--                 array of sw radiation emitted from surface after i-th reflection
-                    IF ( .NOT.  ALLOCATED(surfoutsw_av) )  THEN
-                        ALLOCATE( surfoutsw_av(nsurfl) )
-                        surfoutsw_av = 0.0_wp
-                    ENDIF
-
-                CASE ( 'usm_rad_outlw' )
-!--                 array of lw radiation emitted from surface after i-th reflection
-                    IF ( .NOT.  ALLOCATED(surfoutlw_av) )  THEN
-                        ALLOCATE( surfoutlw_av(nsurfl) )
-                        surfoutlw_av = 0.0_wp
-                    ENDIF
-                CASE ( 'usm_rad_ressw' )
-!--                 array of residua of sw radiation absorbed in surface after last reflection
-                    IF ( .NOT.  ALLOCATED(surfins_av) )  THEN
-                        ALLOCATE( surfins_av(nsurfl) )
-                        surfins_av = 0.0_wp
-                    ENDIF
-                                   
-                CASE ( 'usm_rad_reslw' )
-!--                 array of residua of lw radiation absorbed in surface after last reflection
-                    IF ( .NOT.  ALLOCATED(surfinl_av) )  THEN
-                        ALLOCATE( surfinl_av(nsurfl) )
-                        surfinl_av = 0.0_wp
-                    ENDIF
-                                    
-                CASE ( 'usm_rad_pc_inlw' )
-!--                 array of of lw radiation absorbed in plant canopy
-                    IF ( .NOT.  ALLOCATED(pcbinlw_av) )  THEN
-                        ALLOCATE( pcbinlw_av(1:npcbl) )
-                        pcbinlw_av = 0.0_wp
-                    ENDIF
-                                    
-                CASE ( 'usm_rad_pc_insw' )
-!--                 array of of sw radiation absorbed in plant canopy
-                    IF ( .NOT.  ALLOCATED(pcbinsw_av) )  THEN
-                        ALLOCATE( pcbinsw_av(1:npcbl) )
-                        pcbinsw_av = 0.0_wp
-                    ENDIF
-                                    
-                CASE ( 'usm_rad_pc_inswdir' )
-!--                 array of of direct sw radiation absorbed in plant canopy
-                    IF ( .NOT.  ALLOCATED(pcbinswdir_av) )  THEN
-                        ALLOCATE( pcbinswdir_av(1:npcbl) )
-                        pcbinswdir_av = 0.0_wp
-                    ENDIF
-                                    
-                CASE ( 'usm_rad_pc_inswdif' )
-!--                 array of of diffuse sw radiation absorbed in plant canopy
-                    IF ( .NOT.  ALLOCATED(pcbinswdif_av) )  THEN
-                        ALLOCATE( pcbinswdif_av(1:npcbl) )
-                        pcbinswdif_av = 0.0_wp
-                    ENDIF
-                                    
-                CASE ( 'usm_rad_pc_inswref' )
-!--                 array of of reflected sw radiation absorbed in plant canopy
-                    IF ( .NOT.  ALLOCATED(pcbinswref_av) )  THEN
-                        ALLOCATE( pcbinswref_av(1:npcbl) )
-                        pcbinswref_av = 0.0_wp
-                    ENDIF
-                                   
-                CASE ( 'usm_rad_hf' )
-!--                 array of heat flux from radiation for surfaces after i-th reflection
-                    IF ( l == -1 ) THEN
-                        IF ( .NOT.  ALLOCATED(surf_usm_h%surfhf_av) )  THEN
-                           ALLOCATE( surf_usm_h%surfhf_av(1:surf_usm_h%ns) )
-                           surf_usm_h%surfhf_av = 0.0_wp
-                        ENDIF
-                    ELSE
-                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%surfhf_av) )  THEN
-                           ALLOCATE( surf_usm_v(l)%surfhf_av(1:surf_usm_v(l)%ns) )
-                           surf_usm_v(l)%surfhf_av = 0.0_wp
-                       ENDIF
-                    ENDIF
 
                 CASE ( 'usm_wshf' )
@@ -2218,145 +2059,5 @@
            
            SELECT CASE ( TRIM( var ) )
-                
-                CASE ( 'usm_rad_net' )
-!--                 array of complete radiation balance
-                    IF ( l == -1 ) THEN
-                       DO  m = 1, surf_usm_h%ns
-                          surf_usm_h%rad_net_av(m) =                              &
-                                             surf_usm_h%rad_net_av(m) +           &
-                                             surf_usm_h%rad_net_l(m)
-                       ENDDO
-                    ELSE
-                       DO  m = 1, surf_usm_v(l)%ns
-                          surf_usm_v(l)%rad_net_av(m) =                        &
-                                          surf_usm_v(l)%rad_net_av(m) +        &
-                                          surf_usm_v(l)%rad_net_l(m)
-                       ENDDO
-                    ENDIF
-
-                CASE ( 'usm_rad_insw' )
-!--                 array of sw radiation falling to surface after i-th reflection
-                    DO l = 1, nsurfl
-                        IF ( surfl(id,l) == idsint )  THEN
-                            surfinsw_av(l) = surfinsw_av(l) + surfinsw(l)
-                        ENDIF
-                    ENDDO
-                             
-                CASE ( 'usm_rad_inlw' )
-!--                 array of lw radiation falling to surface after i-th reflection
-                    DO l = 1, nsurfl
-                        IF ( surfl(id,l) == idsint )  THEN
-                            surfinlw_av(l) = surfinlw_av(l) + surfinlw(l)
-                        ENDIF
-                    ENDDO
-                    
-                CASE ( 'usm_rad_inswdir' )
-!--                 array of direct sw radiation falling to surface from sun
-                    DO l = 1, nsurfl
-                        IF ( surfl(id,l) == idsint )  THEN
-                            surfinswdir_av(l) = surfinswdir_av(l) + surfinswdir(l)
-                        ENDIF
-                    ENDDO
-                    
-                CASE ( 'usm_rad_inswdif' )
-!--                 array of difusion sw radiation falling to surface from sky and borders of the domain
-                    DO l = 1, nsurfl
-                        IF ( surfl(id,l) == idsint )  THEN
-                            surfinswdif_av(l) = surfinswdif_av(l) + surfinswdif(l)
-                        ENDIF
-                    ENDDO
-                    
-                CASE ( 'usm_rad_inswref' )
-!--                 array of sw radiation falling to surface from reflections
-                    DO l = 1, nsurfl
-                        IF ( surfl(id,l) == idsint )  THEN
-                            surfinswref_av(l) = surfinswref_av(l) + surfinsw(l) - &
-                                                surfinswdir(l) - surfinswdif(l)
-                        ENDIF
-                    ENDDO
-
-                    
-                CASE ( 'usm_rad_inlwdif' )
-!--                 array of sw radiation falling to surface after i-th reflection
-                    DO l = 1, nsurfl
-                        IF ( surfl(id,l) == idsint )  THEN
-                            surfinlwdif_av(l) = surfinlwdif_av(l) + surfinlwdif(l)
-                        ENDIF
-                    ENDDO
-!                     
-                CASE ( 'usm_rad_inlwref' )
-!--                 array of lw radiation falling to surface from reflections
-                    DO l = 1, nsurfl
-                        IF ( surfl(id,l) == idsint )  THEN
-                            surfinlwref_av(l) = surfinlwref_av(l) + &
-                                                surfinlw(l) - surfinlwdif(l)
-                        ENDIF
-                    ENDDO
-                    
-                CASE ( 'usm_rad_outsw' )
-!--                 array of sw radiation emitted from surface after i-th reflection
-                    DO l = 1, nsurfl
-                        IF ( surfl(id,l) == idsint )  THEN
-                            surfoutsw_av(l) = surfoutsw_av(l) + surfoutsw(l)
-                        ENDIF
-                    ENDDO
-                    
-                CASE ( 'usm_rad_outlw' )
-!--                 array of lw radiation emitted from surface after i-th reflection
-                    DO l = 1, nsurfl
-                        IF ( surfl(id,l) == idsint )  THEN
-                            surfoutlw_av(l) = surfoutlw_av(l) + surfoutlw(l)
-                        ENDIF
-                    ENDDO
-                    
-                CASE ( 'usm_rad_ressw' )
-!--                 array of residua of sw radiation absorbed in surface after last reflection
-                    DO l = 1, nsurfl
-                        IF ( surfl(id,l) == idsint )  THEN
-                            surfins_av(l) = surfins_av(l) + surfins(l)
-                        ENDIF
-                    ENDDO
-                                    
-                CASE ( 'usm_rad_reslw' )
-!--                 array of residua of lw radiation absorbed in surface after last reflection
-                    DO l = 1, nsurfl
-                        IF ( surfl(id,l) == idsint )  THEN
-                            surfinl_av(l) = surfinl_av(l) + surfinl(l)
-                        ENDIF
-                    ENDDO
-                    
-                CASE ( 'usm_rad_pc_inlw' )
-                    pcbinlw_av(:) = pcbinlw_av(:) + pcbinlw(:)
-                    
-                CASE ( 'usm_rad_pc_insw' )
-                    pcbinsw_av(:) = pcbinsw_av(:) + pcbinsw(:)
-                    
-                CASE ( 'usm_rad_pc_inswdir' )
-                    pcbinswdir_av(:) = pcbinswdir_av(:) + pcbinswdir(:)
-                    
-                CASE ( 'usm_rad_pc_inswdif' )
-                    pcbinswdif_av(:) = pcbinswdif_av(:) + pcbinswdif(:)
-                    
-                CASE ( 'usm_rad_pc_inswref' )
-                    pcbinswref_av(:) = pcbinswref_av(:) + pcbinsw(:)     &
-                                                        - pcbinswdir(:)  &
-                                                        - pcbinswdif(:)
-                    
-                CASE ( 'usm_rad_hf' )
-!--                 array of heat flux from radiation for surfaces after i-th reflection
-                    IF ( l == -1 ) THEN
-                       DO  m = 1, surf_usm_h%ns
-                          surf_usm_h%surfhf_av(m) =                               &
-                                             surf_usm_h%surfhf_av(m) +            &
-                                             surf_usm_h%surfhf(m)
-                       ENDDO
-                    ELSE
-                       DO  m = 1, surf_usm_v(l)%ns
-                          surf_usm_v(l)%surfhf_av(m) =                         &
-                                          surf_usm_v(l)%surfhf_av(m) +         &
-                                          surf_usm_v(l)%surfhf(m)
-                       ENDDO
-                    ENDIF
-                    
+
                 CASE ( 'usm_wshf' )
 !--                 array of sensible heat flux from surfaces (land, roof, wall)
@@ -2639,140 +2340,5 @@
            
            SELECT CASE ( TRIM( var ) )
-                
-                CASE ( 'usm_rad_net' )
-!--                 array of complete radiation balance
-                    IF ( l == -1 ) THEN
-                       DO  m = 1, surf_usm_h%ns
-                          surf_usm_h%rad_net_av(m) =                              &
-                                             surf_usm_h%rad_net_av(m) /           &
-                                             REAL( average_count_3d, kind=wp )
-                       ENDDO
-                    ELSE
-                       DO  m = 1, surf_usm_v(l)%ns
-                          surf_usm_v(l)%rad_net_av(m) =                        &
-                                          surf_usm_v(l)%rad_net_av(m) /        &
-                                          REAL( average_count_3d, kind=wp )
-                       ENDDO
-                    ENDIF
-                    
-                CASE ( 'usm_rad_insw' )
-!--                 array of sw radiation falling to surface after i-th reflection
-                    DO l = 1, nsurfl
-                        IF ( surfl(id,l) == idsint )  THEN
-                            surfinsw_av(l) = surfinsw_av(l) / REAL( average_count_3d, kind=wp )
-                        ENDIF
-                    ENDDO
-                             
-                CASE ( 'usm_rad_inlw' )
-!--                 array of lw radiation falling to surface after i-th reflection
-                    DO l = 1, nsurfl
-                        IF ( surfl(id,l) == idsint )  THEN
-                            surfinlw_av(l) = surfinlw_av(l) / REAL( average_count_3d, kind=wp )
-                        ENDIF
-                    ENDDO
-                    
-                CASE ( 'usm_rad_inswdir' )
-!--                 array of direct sw radiation falling to surface from sun
-                    DO l = 1, nsurfl
-                        IF ( surfl(id,l) == idsint )  THEN
-                            surfinswdir_av(l) = surfinswdir_av(l) / REAL( average_count_3d, kind=wp )
-                        ENDIF
-                    ENDDO
-                    
-                CASE ( 'usm_rad_inswdif' )
-!--                 array of difusion sw radiation falling to surface from sky and borders of the domain
-                    DO l = 1, nsurfl
-                        IF ( surfl(id,l) == idsint )  THEN
-                            surfinswdif_av(l) = surfinswdif_av(l) / REAL( average_count_3d, kind=wp )
-                        ENDIF
-                    ENDDO
-                    
-                CASE ( 'usm_rad_inswref' )
-!--                 array of sw radiation falling to surface from reflections
-                    DO l = 1, nsurfl
-                        IF ( surfl(id,l) == idsint )  THEN
-                            surfinswref_av(l) = surfinswref_av(l) / REAL( average_count_3d, kind=wp )
-                        ENDIF
-                    ENDDO
-                    
-                CASE ( 'usm_rad_inlwdif' )
-!--                 array of sw radiation falling to surface after i-th reflection
-                    DO l = 1, nsurfl
-                        IF ( surfl(id,l) == idsint )  THEN
-                            surfinlwdif_av(l) = surfinlwdif_av(l) / REAL( average_count_3d, kind=wp )
-                        ENDIF
-                    ENDDO
-                    
-                CASE ( 'usm_rad_inlwref' )
-!--                 array of lw radiation falling to surface from reflections
-                    DO l = 1, nsurfl
-                        IF ( surfl(id,l) == idsint )  THEN
-                            surfinlwref_av(l) = surfinlwref_av(l) / REAL( average_count_3d, kind=wp )
-                        ENDIF
-                    ENDDO
-                    
-                CASE ( 'usm_rad_outsw' )
-!--                 array of sw radiation emitted from surface after i-th reflection
-                    DO l = 1, nsurfl
-                        IF ( surfl(id,l) == idsint )  THEN
-                            surfoutsw_av(l) = surfoutsw_av(l) / REAL( average_count_3d, kind=wp )
-                        ENDIF
-                    ENDDO
-                    
-                CASE ( 'usm_rad_outlw' )
-!--                 array of lw radiation emitted from surface after i-th reflection
-                    DO l = 1, nsurfl
-                        IF ( surfl(id,l) == idsint )  THEN
-                            surfoutlw_av(l) = surfoutlw_av(l) / REAL( average_count_3d, kind=wp )
-                        ENDIF
-                    ENDDO
-                    
-                CASE ( 'usm_rad_ressw' )
-!--                 array of residua of sw radiation absorbed in surface after last reflection
-                    DO l = 1, nsurfl
-                        IF ( surfl(id,l) == idsint )  THEN
-                            surfins_av(l) = surfins_av(l) / REAL( average_count_3d, kind=wp )
-                        ENDIF
-                    ENDDO
-                                    
-                CASE ( 'usm_rad_reslw' )
-!--                 array of residua of lw radiation absorbed in surface after last reflection
-                    DO l = 1, nsurfl
-                        IF ( surfl(id,l) == idsint )  THEN
-                            surfinl_av(l) = surfinl_av(l) / REAL( average_count_3d, kind=wp )
-                        ENDIF
-                    ENDDO
-                    
-                CASE ( 'usm_rad_pc_inlw' )
-                    pcbinlw_av(:) = pcbinlw_av(:) / REAL( average_count_3d, kind=wp )
-                    
-                CASE ( 'usm_rad_pc_insw' )
-                    pcbinsw_av(:) = pcbinsw_av(:) / REAL( average_count_3d, kind=wp )
-                    
-                CASE ( 'usm_rad_pc_inswdir' )
-                    pcbinswdir_av(:) = pcbinswdir_av(:) / REAL( average_count_3d, kind=wp )
-                    
-                CASE ( 'usm_rad_pc_inswdif' )
-                    pcbinswdif_av(:) = pcbinswdif_av(:) / REAL( average_count_3d, kind=wp )
-                    
-                CASE ( 'usm_rad_pc_inswref' )
-                    pcbinswref_av(:) = pcbinswref_av(:) / REAL( average_count_3d, kind=wp )
-                    
-                CASE ( 'usm_rad_hf' )
-!--                 array of heat flux from radiation for surfaces after i-th reflection
-                    IF ( l == -1 ) THEN
-                       DO  m = 1, surf_usm_h%ns
-                          surf_usm_h%surfhf_av(m) =                               &
-                                             surf_usm_h%surfhf_av(m) /            &
-                                             REAL( average_count_3d, kind=wp )
-                       ENDDO
-                    ELSE
-                       DO  m = 1, surf_usm_v(l)%ns
-                          surf_usm_v(l)%surfhf_av(m) =                         &
-                                          surf_usm_v(l)%surfhf_av(m) /         &
-                                          REAL( average_count_3d, kind=wp )
-                       ENDDO
-                    ENDIF
-                    
+
                 CASE ( 'usm_wshf' )
 !--                 array of sensible heat flux from surfaces (land, roof, wall)
@@ -3115,20 +2681,7 @@
         CHARACTER(LEN=2)                              :: ls
         CHARACTER(LEN=varnamelength)                  :: var          !< TRIM(variable)
-        INTEGER(iwp), PARAMETER                       :: nl1 = 31     !< number of directional usm variables
+        INTEGER(iwp), PARAMETER                       :: nl1 = 16     !< 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_wshf                      ', &
                     'usm_wghf                      ', &
                     'usm_wghf_window               ', &
@@ -3145,20 +2698,8 @@
                     'usm_t_surf_green              ', &
                     'usm_t_green                   ', &
-                    'usm_theta_10cm                ', &
-                    'usm_skyvf                     ', &
-                    'usm_skyvft                    '/)
-
-        INTEGER(iwp), PARAMETER                       :: nl2 = 7      !< number of other variables
-        CHARACTER(LEN=varnamelength), DIMENSION(nl2)  :: varlist2 = & !< list of other usm variables
-                  (/'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                       :: nl3 = 3      !< number of directional layer usm variables
-        CHARACTER(LEN=varnamelength), DIMENSION(nl3)  :: varlist3 = & !< list of directional layer usm variables
+                    'usm_theta_10cm                '/)
+
+        INTEGER(iwp), PARAMETER                       :: nl2 = 3      !< number of directional layer usm variables
+        CHARACTER(LEN=varnamelength), DIMENSION(nl2)  :: varlist2 = & !< list of directional layer usm variables
                   (/'usm_t_wall                    ', &
                     'usm_t_window                  ', &
@@ -3188,9 +2729,9 @@
         IF ( lfound ) GOTO 10
 !       directional layer variables
-        DO i = 1, nl3
+        DO i = 1, nl2
            DO j = 1, nd
               DO l = nzb_wall, nzt_wall
                  WRITE(ls,'(A1,I1)') '_',l
-                 IF ( TRIM(var) == TRIM(varlist3(i))//TRIM(ls)//TRIM(dirname(j)) ) THEN
+                 IF ( TRIM(var) == TRIM(varlist2(i))//TRIM(ls)//TRIM(dirname(j)) ) THEN
                     lfound = .TRUE.
                     EXIT
@@ -3200,12 +2741,4 @@
            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'
@@ -3214,12 +2747,5 @@
 10      CONTINUE
 
-        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.    &
-             var(1:16) == 'usm_rad_inswdif_'  .OR.  var(1:16) == 'usm_rad_inswref_'  .OR. &
-             var(1:16) == 'usm_rad_inlwdif_'  .OR.  var(1:16) == 'usm_rad_inlwref_'  .OR. &
-             var(1:14) == 'usm_rad_outsw_'  .OR.  var(1:14) == 'usm_rad_outlw_'  .OR.     &
-             var(1:14) == 'usm_rad_ressw_'  .OR.  var(1:14) == 'usm_rad_reslw_'  .OR.     &
-             var(1:11) == 'usm_rad_hf_'  .OR.                                             &
-             var(1:9)  == 'usm_wshf_'  .OR.  var(1:9) == 'usm_wghf_' .OR.                 &
+        IF ( var(1:9)  == 'usm_wshf_'  .OR.  var(1:9) == 'usm_wghf_' .OR.                 &
              var(1:16) == 'usm_wghf_window_' .OR. var(1:15) == 'usm_wghf_green_' .OR.     &
              var(1:10) == 'usm_iwghf_' .OR. var(1:17) == 'usm_iwghf_window_'    .OR.      &
@@ -3234,12 +2760,6 @@
                   var(1:14) == 'usm_theta_10cm' )  THEN
             unit = 'K'
-        ELSE IF ( var == 'usm_rad_pc_inlw'  .OR.  var == 'usm_rad_pc_insw'  .OR.          &
-                  var == 'usm_rad_pc_inswdir'  .OR.  var == 'usm_rad_pc_inswdif'  .OR.    &
-                  var == 'usm_rad_pc_inswref' )  THEN
-            unit = 'W'
-        ELSE IF ( var(1:9) == 'usm_surfz'  .OR.  var(1:7) == 'usm_svf'  .OR.              & 
-                  var(1:7) == 'usm_dif'  .OR.  var(1:11) == 'usm_surfcat'  .OR.           &
-                  var(1:11) == 'usm_surfalb'  .OR.  var(1:12) == 'usm_surfemis'  .OR.     &
-                  var(1:9) == 'usm_skyvf' .OR. var(1:9) == 'usm_skyvft' )  THEN
+        ELSE IF ( var(1:9) == 'usm_surfz'  .OR.  var(1:11) == 'usm_surfcat'  .OR.         &
+                  var(1:11) == 'usm_surfalb'  .OR.  var(1:12) == 'usm_surfemis'  )  THEN
             unit = '1'
         ELSE
@@ -3336,11 +2856,6 @@
         INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER             :: diridx =  (/       -1,        1,        0,        3,        2 /)
                                                                      !< index for surf_*_v: 0:3 = (North, South, East, West)
-        INTEGER(iwp), DIMENSION(0:nd-1)                        :: dirstart
-        INTEGER(iwp), DIMENSION(0:nd-1)                        :: dirend
         INTEGER(iwp)                                           :: ids,idsint,idsidx,isurf,isvf,isurfs,isurflt,ipcgb
         INTEGER(iwp)                                           :: is,js,ks,i,j,k,iwl,istat, l, m
-
-        dirstart = (/ startland, startwall, startwall, startwall, startwall /)
-        dirend = (/ endland, endwall, endwall, endwall, endwall /)
 
         found = .TRUE.
@@ -3388,20 +2903,4 @@
             READ(var(9:9), '(I1)', iostat=istat ) iwl
             IF ( istat == 0  .AND.  iwl >= nzb_wall  .AND.  iwl <= nzt_wall )  THEN
-                var = var(1:7)
-            ENDIF
-        ENDIF
-        IF ( (var(1:8) == 'usm_svf_'  .OR.  var(1:8) == 'usm_dif_')  .AND.  len(TRIM(var)) >= 13 )  THEN
-!--         svf values to particular surface
-            surfid = var(9:)
-            i = index(surfid,'_')
-            j = index(surfid(i+1:),'_')
-            READ(surfid(1:i-1),*, iostat=istat ) is
-            IF ( istat == 0 )  THEN
-                READ(surfid(i+1:i+j-1),*, iostat=istat ) js
-            ENDIF
-            IF ( istat == 0 )  THEN
-                READ(surfid(i+j+1:),*, iostat=istat ) ks
-            ENDIF
-            IF ( istat == 0 )  THEN
                 var = var(1:7)
             ENDIF
@@ -3525,301 +3024,4 @@
               ENDIF
 
-          CASE ( 'usm_skyvf' )
-!--           sky view factor
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == idsint )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = skyvf(isurf)
-                 ENDIF
-              ENDDO
-              
-          CASE ( 'usm_skyvft' )
-!--           sky view factor
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == ids )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = skyvft(isurf)
-                 ENDIF
-              ENDDO
-
-!
-!-- Not adjusted so far              
-          CASE ( 'usm_svf', 'usm_dif' )
-!--           shape view factors or iradiance factors to selected surface
-              IF ( TRIM(var)=='usm_svf' )  THEN
-                  k = 1
-              ELSE
-                  k = 2
-              ENDIF
-              DO isvf = 1, nsvfl
-                  isurflt = svfsurf(1, isvf)
-                  isurfs = svfsurf(2, isvf)
-                             
-                  IF ( surf(ix,isurfs) == is  .AND.  surf(iy,isurfs) == js  .AND.       &
-                       surf(iz,isurfs) == ks  .AND.  surf(id,isurfs) == idsint )  THEN
-  !--                 correct source surface
-                      temp_pf(surfl(iz,isurflt),surfl(iy,isurflt),surfl(ix,isurflt)) = svf(k,isvf)
-                  ENDIF
-              ENDDO
-
-          CASE ( 'usm_rad_net' )
-!--           array of complete radiation balance
-              IF ( av == 0 )  THEN
-                 IF ( idsint == iup_u )  THEN
-                    DO  m = 1, surf_usm_h%ns
-                       i = surf_usm_h%i(m)
-                       j = surf_usm_h%j(m)
-                       k = surf_usm_h%k(m)
-                       temp_pf(k,j,i) = surf_usm_h%rad_net_l(m)
-                    ENDDO
-                 ELSE
-                    l = idsidx
-                    DO  m = 1, surf_usm_v(l)%ns
-                       i = surf_usm_v(l)%i(m)
-                       j = surf_usm_v(l)%j(m)
-                       k = surf_usm_v(l)%k(m)
-                       temp_pf(k,j,i) = surf_usm_v(l)%rad_net_l(m)
-                    ENDDO
-                 ENDIF
-              ELSE
-                 IF ( idsint == iup_u )  THEN
-                    DO  m = 1, surf_usm_h%ns
-                       i = surf_usm_h%i(m)
-                       j = surf_usm_h%j(m)
-                       k = surf_usm_h%k(m)
-                       temp_pf(k,j,i) = surf_usm_h%rad_net_av(m)
-                    ENDDO
-                 ELSE
-                    l = idsidx
-                    DO  m = 1, surf_usm_v(l)%ns
-                       i = surf_usm_v(l)%i(m)
-                       j = surf_usm_v(l)%j(m)
-                       k = surf_usm_v(l)%k(m)
-                       temp_pf(k,j,i) = surf_usm_v(l)%rad_net_av(m)
-                    ENDDO
-                 ENDIF
-              ENDIF
-
-          CASE ( 'usm_rad_insw' )
-!--           array of sw radiation falling to surface after i-th reflection
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == idsint )  THEN
-                   IF ( av == 0 )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinsw(isurf)
-                   ELSE
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinsw_av(isurf)
-                   ENDIF
-                 ENDIF
-              ENDDO
-
-          CASE ( 'usm_rad_inlw' )
-!--           array of lw radiation falling to surface after i-th reflection
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == idsint )  THEN
-                   IF ( av == 0 )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw(isurf)
-                   ELSE
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw_av(isurf)
-                   ENDIF
-                 ENDIF
-              ENDDO
-
-          CASE ( 'usm_rad_inswdir' )
-!--           array of direct sw radiation falling to surface from sun
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == idsint )  THEN
-                   IF ( av == 0 )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdir(isurf)
-                   ELSE
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdir_av(isurf)
-                   ENDIF
-                 ENDIF
-              ENDDO
-
-          CASE ( 'usm_rad_inswdif' )
-!--           array of difusion sw radiation falling to surface from sky and borders of the domain
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == idsint )  THEN
-                   IF ( av == 0 )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdif(isurf)
-                   ELSE
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdif_av(isurf)
-                   ENDIF
-                 ENDIF
-              ENDDO
-
-          CASE ( 'usm_rad_inswref' )
-!--           array of sw radiation falling to surface from reflections
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == idsint )  THEN
-                   IF ( av == 0 )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = &
-                       surfinsw(isurf) - surfinswdir(isurf) - surfinswdif(isurf)
-                   ELSE
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswref_av(isurf)
-                   ENDIF
-                 ENDIF
-              ENDDO
-
-          CASE ( 'usm_rad_inlwdif' )
-!--           array of difusion lw radiation falling to surface from sky and borders of the domain
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == idsint )  THEN
-                   IF ( av == 0 )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwdif(isurf)
-                   ELSE
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwdif_av(isurf)
-                   ENDIF
-                 ENDIF
-              ENDDO
-
-          CASE ( 'usm_rad_inlwref' )
-!--           array of lw radiation falling to surface from reflections
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == idsint )  THEN
-                   IF ( av == 0 )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw(isurf) - surfinlwdif(isurf)
-                   ELSE
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwref_av(isurf)
-                   ENDIF
-                 ENDIF
-              ENDDO
-
-          CASE ( 'usm_rad_outsw' )
-!--           array of sw radiation emitted from surface after i-th reflection
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == idsint )  THEN
-                   IF ( av == 0 )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutsw(isurf)
-                   ELSE
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutsw_av(isurf)
-                   ENDIF
-                 ENDIF
-              ENDDO
-
-          CASE ( 'usm_rad_outlw' )
-!--           array of lw radiation emitted from surface after i-th reflection
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == idsint )  THEN
-                   IF ( av == 0 )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutlw(isurf)
-                   ELSE
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutlw_av(isurf)
-                   ENDIF
-                 ENDIF
-              ENDDO
-
-          CASE ( 'usm_rad_ressw' )
-!--           average of array of residua of sw radiation absorbed in surface after last reflection
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == idsint )  THEN
-                   IF ( av == 0 )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfins(isurf)
-                   ELSE
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfins_av(isurf)
-                   ENDIF
-                 ENDIF
-              ENDDO
-
-          CASE ( 'usm_rad_reslw' )
-!--           average of array of residua of lw radiation absorbed in surface after last reflection
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == idsint )  THEN
-                   IF ( av == 0 )  THEN
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinl(isurf)
-                   ELSE
-                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinl_av(isurf)
-                   ENDIF
-                 ENDIF
-              ENDDO
-
-          CASE ( 'usm_rad_pc_inlw' )
-!--           array of lw radiation absorbed by plant canopy
-              DO ipcgb = 1, npcbl
-                  IF ( av == 0 )  THEN
-                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinlw(ipcgb)
-                  ELSE
-                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinlw_av(ipcgb)
-                  ENDIF
-              ENDDO
-
-          CASE ( 'usm_rad_pc_insw' )
-!--           array of sw radiation absorbed by plant canopy
-              DO ipcgb = 1, npcbl
-                  IF ( av == 0 )  THEN
-                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinsw(ipcgb)
-                  ELSE
-                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinsw_av(ipcgb)
-                  ENDIF
-              ENDDO
-
-          CASE ( 'usm_rad_pc_inswdir' )
-!--           array of direct sw radiation absorbed by plant canopy
-              DO ipcgb = 1, npcbl
-                  IF ( av == 0 )  THEN
-                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswdir(ipcgb)
-                  ELSE
-                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswdir_av(ipcgb)
-                  ENDIF
-              ENDDO
-
-          CASE ( 'usm_rad_pc_inswdif' )
-!--           array of diffuse sw radiation absorbed by plant canopy
-              DO ipcgb = 1, npcbl
-                  IF ( av == 0 )  THEN
-                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswdif(ipcgb)
-                  ELSE
-                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswdif_av(ipcgb)
-                  ENDIF
-              ENDDO
-
-          CASE ( 'usm_rad_pc_inswref' )
-!--           array of reflected sw radiation absorbed by plant canopy
-              DO ipcgb = 1, npcbl
-                  IF ( av == 0 )  THEN
-                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinsw(ipcgb)      &
-                                                                              - pcbinswdir(ipcgb) &
-                                                                              - pcbinswdif(ipcgb)
-                  ELSE
-                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswref_av(ipcgb)
-                  ENDIF
-              ENDDO
- 
-          CASE ( 'usm_rad_hf' )
-!--           array of heat flux from radiation for surfaces after all reflections
-              IF ( av == 0 )  THEN
-                 IF ( idsint == iup_u )  THEN
-                    DO  m = 1, surf_usm_h%ns
-                       i = surf_usm_h%i(m)
-                       j = surf_usm_h%j(m)
-                       k = surf_usm_h%k(m)
-                       temp_pf(k,j,i) = surf_usm_h%surfhf(m)
-                    ENDDO
-                 ELSE
-                    l = idsidx
-                    DO  m = 1, surf_usm_v(l)%ns
-                       i = surf_usm_v(l)%i(m)
-                       j = surf_usm_v(l)%j(m)
-                       k = surf_usm_v(l)%k(m)
-                       temp_pf(k,j,i) = surf_usm_v(l)%surfhf(m)
-                    ENDDO
-                 ENDIF
-              ELSE
-                 IF ( idsint == iup_u )  THEN
-                    DO  m = 1, surf_usm_h%ns
-                       i = surf_usm_h%i(m)
-                       j = surf_usm_h%j(m)
-                       k = surf_usm_h%k(m)
-                       temp_pf(k,j,i) = surf_usm_h%surfhf_av(m)
-                    ENDDO
-                 ELSE
-                    l = idsidx
-                    DO  m = 1, surf_usm_v(l)%ns
-                       i = surf_usm_v(l)%i(m)
-                       j = surf_usm_v(l)%j(m)
-                       k = surf_usm_v(l)%k(m)
-                       temp_pf(k,j,i) = surf_usm_v(l)%surfhf_av(m)
-                    ENDDO
-                 ENDIF
-              ENDIF
- 
           CASE ( 'usm_wshf' )
 !--           array of sensible heat flux from surfaces
@@ -3975,7 +3177,4 @@
               ENDIF
 
-
-
-
           CASE ( 'usm_wghf' )
 !--           array of heat flux from ground (land, wall, roof)
@@ -4527,14 +3726,5 @@
 
         var = TRIM(variable)
-        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.      &
-             var(1:16) == 'usm_rad_inswdif_'  .OR.  var(1:16) == 'usm_rad_inswref_'  .OR.   &
-             var(1:16) == 'usm_rad_inlwdif_'  .OR.  var(1:16) == 'usm_rad_inlwref_'  .OR.   &
-             var(1:14) == 'usm_rad_outsw_'  .OR.  var(1:14) == 'usm_rad_outlw_'  .OR.       &
-             var(1:14) == 'usm_rad_ressw_'  .OR.  var(1:14) == 'usm_rad_reslw_'  .OR.       &
-             var(1:11) == 'usm_rad_hf_'  .OR.  var == 'usm_rad_pc_inlw'  .OR.               &
-             var == 'usm_rad_pc_insw'  .OR.  var == 'usm_rad_pc_inswdir'  .OR.              &
-             var == 'usm_rad_pc_inswdif'  .OR.  var == 'usm_rad_pc_inswref'  .OR.           &
-             var(1:9) == 'usm_wshf_'  .OR.  var(1:9) == 'usm_wghf_'  .OR.                   &
+        IF ( var(1:9) == 'usm_wshf_'  .OR.  var(1:9) == 'usm_wghf_'  .OR.                   &
              var(1:16) == 'usm_wghf_window_'  .OR. var(1:15) == 'usm_wghf_green_' .OR.      &
              var(1:10) == 'usm_iwghf_'  .OR. var(1:17) == 'usm_iwghf_window_' .OR.          &
@@ -4545,9 +3735,7 @@
              var(1:16) == 'usm_t_surf_green'  .OR. var(1:11) == 'usm_t_green' .OR.          &
              var(1:15) == 'usm_theta_10cm' .OR.                                             &
-             var(1:9) == 'usm_surfz'  .OR.  var(1:7) == 'usm_svf'  .OR.                     & 
-             var(1:7) == 'usm_dif'  .OR.  var(1:11) == 'usm_surfcat'  .OR.                  &
+             var(1:9) == 'usm_surfz'  .OR.  var(1:11) == 'usm_surfcat'  .OR.                  &
              var(1:11) == 'usm_surfalb'  .OR.  var(1:12) == 'usm_surfemis'  .OR.            &
-             var(1:16) == 'usm_surfwintrans'  .OR. var(1:7) == 'usm_swc' .OR.               &
-             var(1:9) == 'usm_skyvf' .OR. var(1:9) == 'usm_skyvft' ) THEN
+             var(1:16) == 'usm_surfwintrans'  .OR. var(1:7) == 'usm_swc' ) THEN
 
             found = .TRUE.
@@ -10176,5 +9364,5 @@
           CALL wrd_write_string( 't_green_v(' // dum // ')' )
           WRITE ( 14 )  t_green_v(l)%t
-!        
+       
        ENDDO
 
