Index: palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4668)
+++ palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4671)
@@ -28,4 +28,11 @@
 ! -----------------
 ! $Id$
+! Radiative transfer model RTM version 4.1
+! - Implementation of downward facing USM and LSM surfaces
+! - Removal of deprecated CSV inputs
+! - Bugfixes
+! Authors: J. Resler, P. Krc  (Institute of Computer Science, Prague)
+!
+! 4668 2020-09-09 13:00:16Z pavelkrc
 ! Improve debug messages during timestepping
 !
@@ -758,7 +765,10 @@
 
 !-- 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 = (/ iup, isouth, inorth, iwest, ieast /)
+    INTEGER(iwp), PARAMETER                        ::  nd = 6                                  !< number of directions
+    CHARACTER(LEN=6), DIMENSION(0:nd-1), PARAMETER ::  dirname = (/ '_up   ', '_down ', '_south', '_north', '_west ', '_east ' /)
+    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     ::  dirint = (/ iup, idown, isouth, inorth, iwest, ieast /) !< direction integers
+    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     ::  diridx =  (/  0,     1,      1,      0,     3,     2 /) !< mapping to surf_h and surf_v
+    INTEGER(iwp), DIMENSION(0:1), PARAMETER        ::  dirint_h = (/ iup, idown /)                             !< mapping to surf_h directions
+    INTEGER(iwp), DIMENSION(0:3), PARAMETER        ::  dirint_v = (/ inorth, isouth, ieast, iwest /)           !< mapping to surf_v directions
 
 !-- indices and sizes of urban and land surface models
@@ -798,5 +808,5 @@
     REAL(wp), PARAMETER                            ::  ext_coef = 0.6_wp                  !< extinction coefficient (a.k.a. alpha)
     INTEGER(iwp), PARAMETER                        ::  rad_version_len = 10               !< length of identification string of rad version
-    CHARACTER(rad_version_len), PARAMETER          ::  rad_version = 'RAD v. 4.0'         !< identification of version of binary svf and restart files
+    CHARACTER(rad_version_len), PARAMETER          ::  rad_version = 'RAD v. 4.1'         !< identification of version of binary svf and restart files
     INTEGER(iwp)                                   ::  raytrace_discrete_elevs = 40       !< number of discretization steps for elevation (nadir to zenith)
     INTEGER(iwp)                                   ::  raytrace_discrete_azims = 80       !< number of discretization steps for azimuth (out of 360 degrees)
@@ -1107,5 +1117,6 @@
            cos_zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon,   &
            idir, jdir, kdir, id, iz, iy, ix,                                   &
-           iup, inorth, isouth, ieast, iwest,                                  &
+           iup, idown, inorth, isouth, ieast, iwest,                           &
+           nd, dirname, diridx, dirint,                                        &
            nsurf_type, nz_urban_b, nz_urban_t, nz_urban, pch, nsurf,           &
            idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct,                             &
@@ -1183,4 +1194,5 @@
        INTEGER(iwp)                 ::  ilast_word
        INTEGER(iwp)                 ::  ilen
+       INTEGER(iwp)                 ::  id
 
        var = TRIM(variable)
@@ -1188,19 +1200,14 @@
 !--    Identify directional variables
        ilast_word = SCAN(var, '_', back=.TRUE.)
+       directional = .FALSE.
        IF ( ilast_word > 0 )  THEN
-          SELECT CASE ( var(ilast_word:) )
-             CASE ( '_roof', '_south', '_north', '_west', '_east' )
+          DO id = 0, nd-1
+             IF ( TRIM(var(ilast_word:)) == TRIM(dirname(id)) ) THEN
                 directional = .TRUE.
-                write(9, *) 'vardir', var
-                flush(9)
                 var = var(1:ilast_word-1)
-             CASE DEFAULT
-                directional = .FALSE.
-                write(9, *) 'varnd', var
-                flush(9)
-          END SELECT
-       ELSE
-          directional = .FALSE.
-       END IF
+                EXIT
+             ENDIF
+          ENDDO
+       ENDIF
 
        IF ( directional )  THEN
@@ -1213,5 +1220,5 @@
              unit = '1'
           ELSE
-             SELECT CASE ( var )
+             SELECT CASE ( TRIM(var) )
                 CASE ( 'rtm_rad_net', 'rtm_rad_insw', 'rtm_rad_inlw',             &
                        'rtm_rad_inswdir', 'rtm_rad_inswdif', 'rtm_rad_inswref',   &
@@ -1689,14 +1696,16 @@
 !
 !--    Allocate array for storing the surface net radiation
-       IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_net )  .AND.                      &
-                  surf_lsm_h%ns > 0  )   THEN
-          ALLOCATE( surf_lsm_h%rad_net(1:surf_lsm_h%ns) )
-          surf_lsm_h%rad_net = 0.0_wp
-       ENDIF
-       IF ( .NOT. ALLOCATED ( surf_usm_h%rad_net )  .AND.                      &
-                  surf_usm_h%ns > 0  )  THEN
-          ALLOCATE( surf_usm_h%rad_net(1:surf_usm_h%ns) )
-          surf_usm_h%rad_net = 0.0_wp
-       ENDIF
+       DO  l = 0, 1
+          IF ( .NOT. ALLOCATED ( surf_lsm_h(l)%rad_net )  .AND.                      &
+                     surf_lsm_h(l)%ns > 0  )   THEN
+             ALLOCATE( surf_lsm_h(l)%rad_net(1:surf_lsm_h(l)%ns) )
+             surf_lsm_h(l)%rad_net = 0.0_wp
+          ENDIF
+          IF ( .NOT. ALLOCATED ( surf_usm_h(l)%rad_net )  .AND.                      &
+                     surf_usm_h(l)%ns > 0  )  THEN
+             ALLOCATE( surf_usm_h(l)%rad_net(1:surf_usm_h(l)%ns) )
+             surf_usm_h(l)%rad_net = 0.0_wp
+          ENDIF
+       ENDDO
        DO  l = 0, 3
           IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_net )  .AND.                &
@@ -1715,14 +1724,16 @@
 !
 !--    Allocate array for storing the surface longwave (out) radiation change
-       IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_lw_out_change_0 )  .AND.          &
-                  surf_lsm_h%ns > 0  )   THEN
-          ALLOCATE( surf_lsm_h%rad_lw_out_change_0(1:surf_lsm_h%ns) )
-          surf_lsm_h%rad_lw_out_change_0 = 0.0_wp
-       ENDIF
-       IF ( .NOT. ALLOCATED ( surf_usm_h%rad_lw_out_change_0 )  .AND.          &
-                  surf_usm_h%ns > 0  )  THEN
-          ALLOCATE( surf_usm_h%rad_lw_out_change_0(1:surf_usm_h%ns) )
-          surf_usm_h%rad_lw_out_change_0 = 0.0_wp
-       ENDIF
+       DO  l = 0, 1
+          IF ( .NOT. ALLOCATED ( surf_lsm_h(l)%rad_lw_out_change_0 )  .AND.          &
+                     surf_lsm_h(l)%ns > 0  )   THEN
+             ALLOCATE( surf_lsm_h(l)%rad_lw_out_change_0(1:surf_lsm_h(l)%ns) )
+             surf_lsm_h(l)%rad_lw_out_change_0 = 0.0_wp
+          ENDIF
+          IF ( .NOT. ALLOCATED ( surf_usm_h(l)%rad_lw_out_change_0 )  .AND.          &
+                     surf_usm_h(l)%ns > 0  )  THEN
+             ALLOCATE( surf_usm_h(l)%rad_lw_out_change_0(1:surf_usm_h(l)%ns) )
+             surf_usm_h(l)%rad_lw_out_change_0 = 0.0_wp
+          ENDIF
+       ENDDO
        DO  l = 0, 3
           IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_lw_out_change_0 )  .AND.    &
@@ -1740,54 +1751,56 @@
 !
 !--    Allocate surface arrays for incoming/outgoing short/longwave radiation
-       IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_sw_in )  .AND.                    &
-                  surf_lsm_h%ns > 0  )   THEN
-          ALLOCATE( surf_lsm_h%rad_sw_in(1:surf_lsm_h%ns)  )
-          ALLOCATE( surf_lsm_h%rad_sw_out(1:surf_lsm_h%ns) )
-          ALLOCATE( surf_lsm_h%rad_sw_dir(1:surf_lsm_h%ns) )
-          ALLOCATE( surf_lsm_h%rad_sw_dif(1:surf_lsm_h%ns) )
-          ALLOCATE( surf_lsm_h%rad_sw_ref(1:surf_lsm_h%ns) )
-          ALLOCATE( surf_lsm_h%rad_sw_res(1:surf_lsm_h%ns) )
-          ALLOCATE( surf_lsm_h%rad_lw_in(1:surf_lsm_h%ns)  )
-          ALLOCATE( surf_lsm_h%rad_lw_out(1:surf_lsm_h%ns) )
-          ALLOCATE( surf_lsm_h%rad_lw_dif(1:surf_lsm_h%ns) )
-          ALLOCATE( surf_lsm_h%rad_lw_ref(1:surf_lsm_h%ns) )
-          ALLOCATE( surf_lsm_h%rad_lw_res(1:surf_lsm_h%ns) )
-          surf_lsm_h%rad_sw_in  = 0.0_wp
-          surf_lsm_h%rad_sw_out = 0.0_wp
-          surf_lsm_h%rad_sw_dir = 0.0_wp
-          surf_lsm_h%rad_sw_dif = 0.0_wp
-          surf_lsm_h%rad_sw_ref = 0.0_wp
-          surf_lsm_h%rad_sw_res = 0.0_wp
-          surf_lsm_h%rad_lw_in  = 0.0_wp
-          surf_lsm_h%rad_lw_out = 0.0_wp
-          surf_lsm_h%rad_lw_dif = 0.0_wp
-          surf_lsm_h%rad_lw_ref = 0.0_wp
-          surf_lsm_h%rad_lw_res = 0.0_wp
-       ENDIF
-       IF ( .NOT. ALLOCATED ( surf_usm_h%rad_sw_in )  .AND.                    &
-                  surf_usm_h%ns > 0  )  THEN
-          ALLOCATE( surf_usm_h%rad_sw_in(1:surf_usm_h%ns)  )
-          ALLOCATE( surf_usm_h%rad_sw_out(1:surf_usm_h%ns) )
-          ALLOCATE( surf_usm_h%rad_sw_dir(1:surf_usm_h%ns) )
-          ALLOCATE( surf_usm_h%rad_sw_dif(1:surf_usm_h%ns) )
-          ALLOCATE( surf_usm_h%rad_sw_ref(1:surf_usm_h%ns) )
-          ALLOCATE( surf_usm_h%rad_sw_res(1:surf_usm_h%ns) )
-          ALLOCATE( surf_usm_h%rad_lw_in(1:surf_usm_h%ns)  )
-          ALLOCATE( surf_usm_h%rad_lw_out(1:surf_usm_h%ns) )
-          ALLOCATE( surf_usm_h%rad_lw_dif(1:surf_usm_h%ns) )
-          ALLOCATE( surf_usm_h%rad_lw_ref(1:surf_usm_h%ns) )
-          ALLOCATE( surf_usm_h%rad_lw_res(1:surf_usm_h%ns) )
-          surf_usm_h%rad_sw_in  = 0.0_wp
-          surf_usm_h%rad_sw_out = 0.0_wp
-          surf_usm_h%rad_sw_dir = 0.0_wp
-          surf_usm_h%rad_sw_dif = 0.0_wp
-          surf_usm_h%rad_sw_ref = 0.0_wp
-          surf_usm_h%rad_sw_res = 0.0_wp
-          surf_usm_h%rad_lw_in  = 0.0_wp
-          surf_usm_h%rad_lw_out = 0.0_wp
-          surf_usm_h%rad_lw_dif = 0.0_wp
-          surf_usm_h%rad_lw_ref = 0.0_wp
-          surf_usm_h%rad_lw_res = 0.0_wp
-       ENDIF
+       DO  l = 0, 1
+          IF ( .NOT. ALLOCATED ( surf_lsm_h(l)%rad_sw_in )  .AND.                    &
+                     surf_lsm_h(l)%ns > 0  )   THEN
+             ALLOCATE( surf_lsm_h(l)%rad_sw_in(1:surf_lsm_h(l)%ns)  )
+             ALLOCATE( surf_lsm_h(l)%rad_sw_out(1:surf_lsm_h(l)%ns) )
+             ALLOCATE( surf_lsm_h(l)%rad_sw_dir(1:surf_lsm_h(l)%ns) )
+             ALLOCATE( surf_lsm_h(l)%rad_sw_dif(1:surf_lsm_h(l)%ns) )
+             ALLOCATE( surf_lsm_h(l)%rad_sw_ref(1:surf_lsm_h(l)%ns) )
+             ALLOCATE( surf_lsm_h(l)%rad_sw_res(1:surf_lsm_h(l)%ns) )
+             ALLOCATE( surf_lsm_h(l)%rad_lw_in(1:surf_lsm_h(l)%ns)  )
+             ALLOCATE( surf_lsm_h(l)%rad_lw_out(1:surf_lsm_h(l)%ns) )
+             ALLOCATE( surf_lsm_h(l)%rad_lw_dif(1:surf_lsm_h(l)%ns) )
+             ALLOCATE( surf_lsm_h(l)%rad_lw_ref(1:surf_lsm_h(l)%ns) )
+             ALLOCATE( surf_lsm_h(l)%rad_lw_res(1:surf_lsm_h(l)%ns) )
+             surf_lsm_h(l)%rad_sw_in  = 0.0_wp
+             surf_lsm_h(l)%rad_sw_out = 0.0_wp
+             surf_lsm_h(l)%rad_sw_dir = 0.0_wp
+             surf_lsm_h(l)%rad_sw_dif = 0.0_wp
+             surf_lsm_h(l)%rad_sw_ref = 0.0_wp
+             surf_lsm_h(l)%rad_sw_res = 0.0_wp
+             surf_lsm_h(l)%rad_lw_in  = 0.0_wp
+             surf_lsm_h(l)%rad_lw_out = 0.0_wp
+             surf_lsm_h(l)%rad_lw_dif = 0.0_wp
+             surf_lsm_h(l)%rad_lw_ref = 0.0_wp
+             surf_lsm_h(l)%rad_lw_res = 0.0_wp
+          ENDIF
+          IF ( .NOT. ALLOCATED ( surf_usm_h(l)%rad_sw_in )  .AND.                    &
+                     surf_usm_h(l)%ns > 0  )  THEN
+             ALLOCATE( surf_usm_h(l)%rad_sw_in(1:surf_usm_h(l)%ns)  )
+             ALLOCATE( surf_usm_h(l)%rad_sw_out(1:surf_usm_h(l)%ns) )
+             ALLOCATE( surf_usm_h(l)%rad_sw_dir(1:surf_usm_h(l)%ns) )
+             ALLOCATE( surf_usm_h(l)%rad_sw_dif(1:surf_usm_h(l)%ns) )
+             ALLOCATE( surf_usm_h(l)%rad_sw_ref(1:surf_usm_h(l)%ns) )
+             ALLOCATE( surf_usm_h(l)%rad_sw_res(1:surf_usm_h(l)%ns) )
+             ALLOCATE( surf_usm_h(l)%rad_lw_in(1:surf_usm_h(l)%ns)  )
+             ALLOCATE( surf_usm_h(l)%rad_lw_out(1:surf_usm_h(l)%ns) )
+             ALLOCATE( surf_usm_h(l)%rad_lw_dif(1:surf_usm_h(l)%ns) )
+             ALLOCATE( surf_usm_h(l)%rad_lw_ref(1:surf_usm_h(l)%ns) )
+             ALLOCATE( surf_usm_h(l)%rad_lw_res(1:surf_usm_h(l)%ns) )
+             surf_usm_h(l)%rad_sw_in  = 0.0_wp
+             surf_usm_h(l)%rad_sw_out = 0.0_wp
+             surf_usm_h(l)%rad_sw_dir = 0.0_wp
+             surf_usm_h(l)%rad_sw_dif = 0.0_wp
+             surf_usm_h(l)%rad_sw_ref = 0.0_wp
+             surf_usm_h(l)%rad_sw_res = 0.0_wp
+             surf_usm_h(l)%rad_lw_in  = 0.0_wp
+             surf_usm_h(l)%rad_lw_out = 0.0_wp
+             surf_usm_h(l)%rad_lw_dif = 0.0_wp
+             surf_usm_h(l)%rad_lw_ref = 0.0_wp
+             surf_usm_h(l)%rad_lw_res = 0.0_wp
+          ENDIF
+       ENDDO
        DO  l = 0, 3
           IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_sw_in )  .AND.              &
@@ -1848,8 +1861,10 @@
 !--    Fix net radiation in case of radiation_scheme = 'constant'
        IF ( radiation_scheme == 'constant' )  THEN
-          IF ( ALLOCATED( surf_lsm_h%rad_net ) )                               &
-             surf_lsm_h%rad_net    = net_radiation
-          IF ( ALLOCATED( surf_usm_h%rad_net ) )                               &
-             surf_usm_h%rad_net    = net_radiation
+          DO  l = 0, 1
+             IF ( ALLOCATED( surf_lsm_h(l)%rad_net ) )                               &
+                surf_lsm_h(l)%rad_net    = net_radiation
+             IF ( ALLOCATED( surf_usm_h(l)%rad_net ) )                               &
+                surf_usm_h(l)%rad_net    = net_radiation
+          ENDDO
 !
 !--       Todo: weight with inclination angle
@@ -1912,13 +1927,14 @@
 !--       Allocate arrays for broadband albedo, and level 1 initialization
 !--       via namelist paramter, unless not already allocated.
-          IF ( .NOT. ALLOCATED(surf_lsm_h%albedo) )  THEN
-             ALLOCATE( surf_lsm_h%albedo(1:surf_lsm_h%ns,0:2)     )
-             surf_lsm_h%albedo    = albedo
-          ENDIF
-          IF ( .NOT. ALLOCATED(surf_usm_h%albedo) )  THEN
-             ALLOCATE( surf_usm_h%albedo(1:surf_usm_h%ns,0:2)     )
-             surf_usm_h%albedo    = albedo
-          ENDIF
-
+          DO  l = 0, 1
+             IF ( .NOT. ALLOCATED(surf_lsm_h(l)%albedo) )  THEN
+                ALLOCATE( surf_lsm_h(l)%albedo(1:surf_lsm_h(l)%ns,0:2)     )
+                surf_lsm_h(l)%albedo    = albedo
+             ENDIF
+             IF ( .NOT. ALLOCATED(surf_usm_h(l)%albedo) )  THEN
+                ALLOCATE( surf_usm_h(l)%albedo(1:surf_usm_h(l)%ns,0:2)     )
+                surf_usm_h(l)%albedo    = albedo
+             ENDIF
+          ENDDO
           DO  l = 0, 3
              IF ( .NOT. ALLOCATED( surf_lsm_v(l)%albedo ) )  THEN
@@ -1936,27 +1952,28 @@
 !--       input data is read from ASCII file, albedo_type will be zero, so that
 !--       albedo won't be overwritten.
-          DO  m = 1, surf_lsm_h%ns
-             IF ( surf_lsm_h%albedo_type(m,ind_veg_wall) /= 0 )                &
-                surf_lsm_h%albedo(m,ind_veg_wall) =                            &
-                           albedo_pars(0,surf_lsm_h%albedo_type(m,ind_veg_wall))
-             IF ( surf_lsm_h%albedo_type(m,ind_pav_green) /= 0 )               &
-                surf_lsm_h%albedo(m,ind_pav_green) =                           &
-                           albedo_pars(0,surf_lsm_h%albedo_type(m,ind_pav_green))
-             IF ( surf_lsm_h%albedo_type(m,ind_wat_win) /= 0 )                 &
-                surf_lsm_h%albedo(m,ind_wat_win) =                             &
-                           albedo_pars(0,surf_lsm_h%albedo_type(m,ind_wat_win))
+          DO  l = 0, 1
+             DO  m = 1, surf_lsm_h(l)%ns
+                IF ( surf_lsm_h(l)%albedo_type(m,ind_veg_wall) /= 0 )                &
+                   surf_lsm_h(l)%albedo(m,ind_veg_wall) =                            &
+                              albedo_pars(0,surf_lsm_h(l)%albedo_type(m,ind_veg_wall))
+                IF ( surf_lsm_h(l)%albedo_type(m,ind_pav_green) /= 0 )               &
+                   surf_lsm_h(l)%albedo(m,ind_pav_green) =                           &
+                              albedo_pars(0,surf_lsm_h(l)%albedo_type(m,ind_pav_green))
+                IF ( surf_lsm_h(l)%albedo_type(m,ind_wat_win) /= 0 )                 &
+                   surf_lsm_h(l)%albedo(m,ind_wat_win) =                             &
+                              albedo_pars(0,surf_lsm_h(l)%albedo_type(m,ind_wat_win))
+             ENDDO
+             DO  m = 1, surf_usm_h(l)%ns
+                IF ( surf_usm_h(l)%albedo_type(m,ind_veg_wall) /= 0 )                &
+                   surf_usm_h(l)%albedo(m,ind_veg_wall) =                            &
+                              albedo_pars(0,surf_usm_h(l)%albedo_type(m,ind_veg_wall))
+                IF ( surf_usm_h(l)%albedo_type(m,ind_pav_green) /= 0 )               &
+                   surf_usm_h(l)%albedo(m,ind_pav_green) =                           &
+                              albedo_pars(0,surf_usm_h(l)%albedo_type(m,ind_pav_green))
+                IF ( surf_usm_h(l)%albedo_type(m,ind_wat_win) /= 0 )                 &
+                   surf_usm_h(l)%albedo(m,ind_wat_win) =                             &
+                              albedo_pars(0,surf_usm_h(l)%albedo_type(m,ind_wat_win))
+             ENDDO
           ENDDO
-          DO  m = 1, surf_usm_h%ns
-             IF ( surf_usm_h%albedo_type(m,ind_veg_wall) /= 0 )                &
-                surf_usm_h%albedo(m,ind_veg_wall) =                            &
-                           albedo_pars(0,surf_usm_h%albedo_type(m,ind_veg_wall))
-             IF ( surf_usm_h%albedo_type(m,ind_pav_green) /= 0 )               &
-                surf_usm_h%albedo(m,ind_pav_green) =                           &
-                           albedo_pars(0,surf_usm_h%albedo_type(m,ind_pav_green))
-             IF ( surf_usm_h%albedo_type(m,ind_wat_win) /= 0 )                 &
-                surf_usm_h%albedo(m,ind_wat_win) =                             &
-                           albedo_pars(0,surf_usm_h%albedo_type(m,ind_wat_win))
-          ENDDO
-
           DO  l = 0, 3
              DO  m = 1, surf_lsm_v(l)%ns
@@ -1991,21 +2008,23 @@
 !
 !--          Horizontal surfaces
-             DO  m = 1, surf_lsm_h%ns
-                i = surf_lsm_h%i(m)
-                j = surf_lsm_h%j(m)
-                IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
-                   surf_lsm_h%albedo(m,ind_veg_wall)  = albedo_pars_f%pars_xy(0,j,i)
-                   surf_lsm_h%albedo(m,ind_pav_green) = albedo_pars_f%pars_xy(0,j,i)
-                   surf_lsm_h%albedo(m,ind_wat_win)   = albedo_pars_f%pars_xy(0,j,i)
-                ENDIF
-             ENDDO
-             DO  m = 1, surf_usm_h%ns
-                i = surf_usm_h%i(m)
-                j = surf_usm_h%j(m)
-                IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
-                   surf_usm_h%albedo(m,ind_veg_wall)  = albedo_pars_f%pars_xy(0,j,i)
-                   surf_usm_h%albedo(m,ind_pav_green) = albedo_pars_f%pars_xy(0,j,i)
-                   surf_usm_h%albedo(m,ind_wat_win)   = albedo_pars_f%pars_xy(0,j,i)
-                ENDIF
+             DO  l = 0, 1
+                DO  m = 1, surf_lsm_h(l)%ns
+                   i = surf_lsm_h(l)%i(m)
+                   j = surf_lsm_h(l)%j(m)
+                   IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
+                      surf_lsm_h(l)%albedo(m,ind_veg_wall)  = albedo_pars_f%pars_xy(0,j,i)
+                      surf_lsm_h(l)%albedo(m,ind_pav_green) = albedo_pars_f%pars_xy(0,j,i)
+                      surf_lsm_h(l)%albedo(m,ind_wat_win)   = albedo_pars_f%pars_xy(0,j,i)
+                   ENDIF
+                ENDDO
+                DO  m = 1, surf_usm_h(l)%ns
+                   i = surf_usm_h(l)%i(m)
+                   j = surf_usm_h(l)%j(m)
+                   IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
+                      surf_usm_h(l)%albedo(m,ind_veg_wall)  = albedo_pars_f%pars_xy(0,j,i)
+                      surf_usm_h(l)%albedo(m,ind_pav_green) = albedo_pars_f%pars_xy(0,j,i)
+                      surf_usm_h(l)%albedo(m,ind_wat_win)   = albedo_pars_f%pars_xy(0,j,i)
+                   ENDIF
+                ENDDO
              ENDDO
 !
@@ -2044,38 +2063,40 @@
 !--       to zero in order to take effect.
           IF ( building_surface_pars_f%from_file )  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)
-!
-!--             Iterate over surfaces in column, check height and orientation
-                DO  is = building_surface_pars_f%index_ji(1,j,i), &
-                         building_surface_pars_f%index_ji(2,j,i)
-                   IF ( building_surface_pars_f%coords(4,is) == -surf_usm_h%koff .AND. &
-                        building_surface_pars_f%coords(1,is) == k )  THEN
-
-                      IF ( building_surface_pars_f%pars(ind_s_alb_b_wall,is) /=      &
-                           building_surface_pars_f%fill )  THEN
-                         surf_usm_h%albedo(m,ind_veg_wall) =                         &
-                                  building_surface_pars_f%pars(ind_s_alb_b_wall,is)
-                         surf_usm_h%albedo_type(m,ind_veg_wall) = 0
+             DO l = 0, 1
+                DO  m = 1, surf_usm_h(l)%ns
+                   i = surf_usm_h(l)%i(m)
+                   j = surf_usm_h(l)%j(m)
+                   k = surf_usm_h(l)%k(m)
+!
+!--                Iterate over surfaces in column, check height and orientation
+                   DO  is = building_surface_pars_f%index_ji(1,j,i), &
+                            building_surface_pars_f%index_ji(2,j,i)
+                      IF ( building_surface_pars_f%coords(4,is) == -surf_usm_h(l)%koff .AND. &
+                           building_surface_pars_f%coords(1,is) == k )  THEN
+
+                         IF ( building_surface_pars_f%pars(ind_s_alb_b_wall,is) /=           &
+                              building_surface_pars_f%fill )  THEN
+                            surf_usm_h(l)%albedo(m,ind_veg_wall) =                           &
+                                     building_surface_pars_f%pars(ind_s_alb_b_wall,is)
+                            surf_usm_h(l)%albedo_type(m,ind_veg_wall) = 0
+                         ENDIF
+
+                         IF ( building_surface_pars_f%pars(ind_s_alb_b_win,is) /=            &
+                              building_surface_pars_f%fill )  THEN
+                            surf_usm_h(l)%albedo(m,ind_wat_win) =                            &
+                                     building_surface_pars_f%pars(ind_s_alb_b_win,is)
+                            surf_usm_h(l)%albedo_type(m,ind_wat_win) = 0
+                         ENDIF
+
+                         IF ( building_surface_pars_f%pars(ind_s_alb_b_green,is) /=          &
+                              building_surface_pars_f%fill )  THEN
+                            surf_usm_h(l)%albedo(m,ind_pav_green) =                          &
+                                     building_surface_pars_f%pars(ind_s_alb_b_green,is)
+                            surf_usm_h(l)%albedo_type(m,ind_pav_green) = 0
+                         ENDIF
+
+                         EXIT ! surface was found and processed
                       ENDIF
-
-                      IF ( building_surface_pars_f%pars(ind_s_alb_b_win,is) /=       &
-                           building_surface_pars_f%fill )  THEN
-                         surf_usm_h%albedo(m,ind_wat_win) =                          &
-                                  building_surface_pars_f%pars(ind_s_alb_b_win,is)
-                         surf_usm_h%albedo_type(m,ind_wat_win) = 0
-                      ENDIF
-
-                      IF ( building_surface_pars_f%pars(ind_s_alb_b_green,is) /=     &
-                           building_surface_pars_f%fill )  THEN
-                         surf_usm_h%albedo(m,ind_pav_green) =                        &
-                                  building_surface_pars_f%pars(ind_s_alb_b_green,is)
-                         surf_usm_h%albedo_type(m,ind_pav_green) = 0
-                      ENDIF
-
-                      EXIT ! surface was found and processed
-                   ENDIF
+                   ENDDO
                 ENDDO
              ENDDO
@@ -2129,30 +2150,31 @@
 !--       for wall/green/window (USM) or vegetation/pavement/water surfaces
 !--       (LSM).
-          ALLOCATE ( surf_lsm_h%aldif(1:surf_lsm_h%ns,0:2)       )
-          ALLOCATE ( surf_lsm_h%aldir(1:surf_lsm_h%ns,0:2)       )
-          ALLOCATE ( surf_lsm_h%asdif(1:surf_lsm_h%ns,0:2)       )
-          ALLOCATE ( surf_lsm_h%asdir(1:surf_lsm_h%ns,0:2)       )
-          ALLOCATE ( surf_lsm_h%rrtm_aldif(1:surf_lsm_h%ns,0:2)  )
-          ALLOCATE ( surf_lsm_h%rrtm_aldir(1:surf_lsm_h%ns,0:2)  )
-          ALLOCATE ( surf_lsm_h%rrtm_asdif(1:surf_lsm_h%ns,0:2)  )
-          ALLOCATE ( surf_lsm_h%rrtm_asdir(1:surf_lsm_h%ns,0:2)  )
-
-          ALLOCATE ( surf_usm_h%aldif(1:surf_usm_h%ns,0:2)       )
-          ALLOCATE ( surf_usm_h%aldir(1:surf_usm_h%ns,0:2)       )
-          ALLOCATE ( surf_usm_h%asdif(1:surf_usm_h%ns,0:2)       )
-          ALLOCATE ( surf_usm_h%asdir(1:surf_usm_h%ns,0:2)       )
-          ALLOCATE ( surf_usm_h%rrtm_aldif(1:surf_usm_h%ns,0:2)  )
-          ALLOCATE ( surf_usm_h%rrtm_aldir(1:surf_usm_h%ns,0:2)  )
-          ALLOCATE ( surf_usm_h%rrtm_asdif(1:surf_usm_h%ns,0:2)  )
-          ALLOCATE ( surf_usm_h%rrtm_asdir(1:surf_usm_h%ns,0:2)  )
-
-!
-!--       Allocate broadband albedo (temporary for the current radiation
-!--       implementations)
-          IF ( .NOT. ALLOCATED(surf_lsm_h%albedo) )                            &
-             ALLOCATE( surf_lsm_h%albedo(1:surf_lsm_h%ns,0:2)     )
-          IF ( .NOT. ALLOCATED(surf_usm_h%albedo) )                            &
-             ALLOCATE( surf_usm_h%albedo(1:surf_usm_h%ns,0:2)     )
-
+          DO  l = 0, 1
+             ALLOCATE ( surf_lsm_h(l)%aldif(1:surf_lsm_h(l)%ns,0:2)       )
+             ALLOCATE ( surf_lsm_h(l)%aldir(1:surf_lsm_h(l)%ns,0:2)       )
+             ALLOCATE ( surf_lsm_h(l)%asdif(1:surf_lsm_h(l)%ns,0:2)       )
+             ALLOCATE ( surf_lsm_h(l)%asdir(1:surf_lsm_h(l)%ns,0:2)       )
+             ALLOCATE ( surf_lsm_h(l)%rrtm_aldif(1:surf_lsm_h(l)%ns,0:2)  )
+             ALLOCATE ( surf_lsm_h(l)%rrtm_aldir(1:surf_lsm_h(l)%ns,0:2)  )
+             ALLOCATE ( surf_lsm_h(l)%rrtm_asdif(1:surf_lsm_h(l)%ns,0:2)  )
+             ALLOCATE ( surf_lsm_h(l)%rrtm_asdir(1:surf_lsm_h(l)%ns,0:2)  )
+
+             ALLOCATE ( surf_usm_h(l)%aldif(1:surf_usm_h(l)%ns,0:2)       )
+             ALLOCATE ( surf_usm_h(l)%aldir(1:surf_usm_h(l)%ns,0:2)       )
+             ALLOCATE ( surf_usm_h(l)%asdif(1:surf_usm_h(l)%ns,0:2)       )
+             ALLOCATE ( surf_usm_h(l)%asdir(1:surf_usm_h(l)%ns,0:2)       )
+             ALLOCATE ( surf_usm_h(l)%rrtm_aldif(1:surf_usm_h(l)%ns,0:2)  )
+             ALLOCATE ( surf_usm_h(l)%rrtm_aldir(1:surf_usm_h(l)%ns,0:2)  )
+             ALLOCATE ( surf_usm_h(l)%rrtm_asdif(1:surf_usm_h(l)%ns,0:2)  )
+             ALLOCATE ( surf_usm_h(l)%rrtm_asdir(1:surf_usm_h(l)%ns,0:2)  )
+
+!
+!--          Allocate broadband albedo (temporary for the current radiation
+!--          implementations)
+             IF ( .NOT. ALLOCATED(surf_lsm_h(l)%albedo) )                            &
+                ALLOCATE( surf_lsm_h(l)%albedo(1:surf_lsm_h(l)%ns,0:2)     )
+             IF ( .NOT. ALLOCATED(surf_usm_h(l)%albedo) )                            &
+                ALLOCATE( surf_usm_h(l)%albedo(1:surf_usm_h(l)%ns,0:2)     )
+          ENDDO
 !
 !--       Allocate albedos for short/longwave radiation, vertical surfaces
@@ -2191,26 +2213,27 @@
 !--       paramters. Please note, this case all surface tiles are initialized
 !--       the same.
-          IF ( surf_lsm_h%ns > 0 )  THEN
-             surf_lsm_h%aldif  = albedo_lw_dif
-             surf_lsm_h%aldir  = albedo_lw_dir
-             surf_lsm_h%asdif  = albedo_sw_dif
-             surf_lsm_h%asdir  = albedo_sw_dir
-             surf_lsm_h%albedo = albedo_sw_dif
-          ENDIF
-          IF ( surf_usm_h%ns > 0 )  THEN
-             IF ( surf_usm_h%albedo_from_ascii )  THEN
-                surf_usm_h%aldif  = surf_usm_h%albedo
-                surf_usm_h%aldir  = surf_usm_h%albedo
-                surf_usm_h%asdif  = surf_usm_h%albedo
-                surf_usm_h%asdir  = surf_usm_h%albedo
-             ELSE
-                surf_usm_h%aldif  = albedo_lw_dif
-                surf_usm_h%aldir  = albedo_lw_dir
-                surf_usm_h%asdif  = albedo_sw_dif
-                surf_usm_h%asdir  = albedo_sw_dir
-                surf_usm_h%albedo = albedo_sw_dif
+          DO  l = 0, 1
+             IF ( surf_lsm_h(l)%ns > 0 )  THEN
+                surf_lsm_h(l)%aldif  = albedo_lw_dif
+                surf_lsm_h(l)%aldir  = albedo_lw_dir
+                surf_lsm_h(l)%asdif  = albedo_sw_dif
+                surf_lsm_h(l)%asdir  = albedo_sw_dir
+                surf_lsm_h(l)%albedo = albedo_sw_dif
              ENDIF
-          ENDIF
-
+             IF ( surf_usm_h(l)%ns > 0 )  THEN
+                IF ( surf_usm_h(l)%albedo_from_ascii )  THEN
+                   surf_usm_h(l)%aldif  = surf_usm_h(l)%albedo
+                   surf_usm_h(l)%aldir  = surf_usm_h(l)%albedo
+                   surf_usm_h(l)%asdif  = surf_usm_h(l)%albedo
+                   surf_usm_h(l)%asdir  = surf_usm_h(l)%albedo
+                ELSE
+                   surf_usm_h(l)%aldif  = albedo_lw_dif
+                   surf_usm_h(l)%aldir  = albedo_lw_dir
+                   surf_usm_h(l)%asdif  = albedo_sw_dif
+                   surf_usm_h(l)%asdir  = albedo_sw_dir
+                   surf_usm_h(l)%albedo = albedo_sw_dif
+                ENDIF
+             ENDIF
+          ENDDO
           DO  l = 0, 3
 
@@ -2243,50 +2266,50 @@
 !--       is applied so that the resulting albedo is calculated via the weighted
 !--       average of respective surface fractions.
-          DO  m = 1, surf_lsm_h%ns
-!
-!--          Spectral albedos for vegetation/pavement/water surfaces
-             DO  ind_type = 0, 2
-                IF ( surf_lsm_h%albedo_type(m,ind_type) /= 0 )  THEN
-                   surf_lsm_h%aldif(m,ind_type) =                              &
-                               albedo_pars(1,surf_lsm_h%albedo_type(m,ind_type))
-                   surf_lsm_h%asdif(m,ind_type) =                              &
-                               albedo_pars(2,surf_lsm_h%albedo_type(m,ind_type))
-                   surf_lsm_h%aldir(m,ind_type) =                              &
-                               albedo_pars(1,surf_lsm_h%albedo_type(m,ind_type))
-                   surf_lsm_h%asdir(m,ind_type) =                              &
-                               albedo_pars(2,surf_lsm_h%albedo_type(m,ind_type))
-                   surf_lsm_h%albedo(m,ind_type) =                             &
-                               albedo_pars(0,surf_lsm_h%albedo_type(m,ind_type))
-                ENDIF
-             ENDDO
-
-          ENDDO
-!
-!--       For urban surface only if albedo has not been already initialized
-!--       in the urban-surface model via the ASCII file.
-          IF ( .NOT. surf_usm_h%albedo_from_ascii )  THEN
-             DO  m = 1, surf_usm_h%ns
-!
-!--             Spectral albedos for wall/green/window surfaces
+          DO l = 0, 1
+             DO  m = 1, surf_lsm_h(l)%ns
+!
+!--             Spectral albedos for vegetation/pavement/water surfaces
                 DO  ind_type = 0, 2
-                   IF ( surf_usm_h%albedo_type(m,ind_type) /= 0 )  THEN
-                      surf_usm_h%aldif(m,ind_type) =                           &
-                               albedo_pars(1,surf_usm_h%albedo_type(m,ind_type))
-                      surf_usm_h%asdif(m,ind_type) =                           &
-                               albedo_pars(2,surf_usm_h%albedo_type(m,ind_type))
-                      surf_usm_h%aldir(m,ind_type) =                           &
-                               albedo_pars(1,surf_usm_h%albedo_type(m,ind_type))
-                      surf_usm_h%asdir(m,ind_type) =                           &
-                               albedo_pars(2,surf_usm_h%albedo_type(m,ind_type))
-                      surf_usm_h%albedo(m,ind_type) =                          &
-                               albedo_pars(0,surf_usm_h%albedo_type(m,ind_type))
+                   IF ( surf_lsm_h(l)%albedo_type(m,ind_type) /= 0 )  THEN
+                      surf_lsm_h(l)%aldif(m,ind_type) =                              &
+                                  albedo_pars(1,surf_lsm_h(l)%albedo_type(m,ind_type))
+                      surf_lsm_h(l)%asdif(m,ind_type) =                              &
+                                  albedo_pars(2,surf_lsm_h(l)%albedo_type(m,ind_type))
+                      surf_lsm_h(l)%aldir(m,ind_type) =                              &
+                                  albedo_pars(1,surf_lsm_h(l)%albedo_type(m,ind_type))
+                      surf_lsm_h(l)%asdir(m,ind_type) =                              &
+                                  albedo_pars(2,surf_lsm_h(l)%albedo_type(m,ind_type))
+                      surf_lsm_h(l)%albedo(m,ind_type) =                             &
+                                  albedo_pars(0,surf_lsm_h(l)%albedo_type(m,ind_type))
                    ENDIF
                 ENDDO
 
              ENDDO
-          ENDIF
+!
+!--          For urban surface only if albedo has not been already initialized
+!--          in the urban-surface model via the ASCII file.
+             IF ( .NOT. surf_usm_h(l)%albedo_from_ascii )  THEN
+                DO  m = 1, surf_usm_h(l)%ns
+!
+!--                Spectral albedos for wall/green/window surfaces
+                   DO  ind_type = 0, 2
+                      IF ( surf_usm_h(l)%albedo_type(m,ind_type) /= 0 )  THEN
+                         surf_usm_h(l)%aldif(m,ind_type) =                           &
+                                  albedo_pars(1,surf_usm_h(l)%albedo_type(m,ind_type))
+                         surf_usm_h(l)%asdif(m,ind_type) =                           &
+                                  albedo_pars(2,surf_usm_h(l)%albedo_type(m,ind_type))
+                         surf_usm_h(l)%aldir(m,ind_type) =                           &
+                                  albedo_pars(1,surf_usm_h(l)%albedo_type(m,ind_type))
+                         surf_usm_h(l)%asdir(m,ind_type) =                           &
+                                  albedo_pars(2,surf_usm_h(l)%albedo_type(m,ind_type))
+                         surf_usm_h(l)%albedo(m,ind_type) =                          &
+                                  albedo_pars(0,surf_usm_h(l)%albedo_type(m,ind_type))
+                      ENDIF
+                   ENDDO
+                ENDDO
+             ENDIF
+          ENDDO
 
           DO l = 0, 3
-
              DO  m = 1, surf_lsm_v(l)%ns
 !
@@ -2338,86 +2361,88 @@
 !
 !--          Horizontal
-             DO  m = 1, surf_lsm_h%ns
-                i = surf_lsm_h%i(m)
-                j = surf_lsm_h%j(m)
-!
-!--             Spectral albedos for vegetation/pavement/water surfaces
-                DO  ind_type = 0, 2
-                   IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )   &
-                      surf_lsm_h%albedo(m,ind_type) =                          &
-                                             albedo_pars_f%pars_xy(0,j,i)
-                   IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )   &
-                      surf_lsm_h%aldir(m,ind_type) =                           &
-                                             albedo_pars_f%pars_xy(1,j,i)
-                   IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )   &
-                      surf_lsm_h%aldif(m,ind_type) =                           &
-                                             albedo_pars_f%pars_xy(1,j,i)
-                   IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill )   &
-                      surf_lsm_h%asdir(m,ind_type) =                           &
-                                             albedo_pars_f%pars_xy(2,j,i)
-                   IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill )   &
-                      surf_lsm_h%asdif(m,ind_type) =                           &
-                                             albedo_pars_f%pars_xy(2,j,i)
+             DO l = 0, 1
+                DO  m = 1, surf_lsm_h(l)%ns
+                   i = surf_lsm_h(l)%i(m)
+                   j = surf_lsm_h(l)%j(m)
+   !
+   !--             Spectral albedos for vegetation/pavement/water surfaces
+                   DO  ind_type = 0, 2
+                      IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )   &
+                         surf_lsm_h(l)%albedo(m,ind_type) =                          &
+                                                albedo_pars_f%pars_xy(0,j,i)
+                      IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )   &
+                         surf_lsm_h(l)%aldir(m,ind_type) =                           &
+                                                albedo_pars_f%pars_xy(1,j,i)
+                      IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )   &
+                         surf_lsm_h(l)%aldif(m,ind_type) =                           &
+                                                albedo_pars_f%pars_xy(1,j,i)
+                      IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill )   &
+                         surf_lsm_h(l)%asdir(m,ind_type) =                           &
+                                                albedo_pars_f%pars_xy(2,j,i)
+                      IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill )   &
+                         surf_lsm_h(l)%asdif(m,ind_type) =                           &
+                                                albedo_pars_f%pars_xy(2,j,i)
+                   ENDDO
                 ENDDO
+!
+!--             For urban surface only if albedo has not been already initialized
+!--             in the urban-surface model via the ASCII file.
+                IF ( .NOT. surf_usm_h(l)%albedo_from_ascii )  THEN
+                   DO  m = 1, surf_usm_h(l)%ns
+                      i = surf_usm_h(l)%i(m)
+                      j = surf_usm_h(l)%j(m)
+!
+!--                   Broadband albedos for wall/green/window surfaces
+                      DO  ind_type = 0, 2
+                         IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )&
+                            surf_usm_h(l)%albedo(m,ind_type) =                       &
+                                                albedo_pars_f%pars_xy(0,j,i)
+                      ENDDO
+!
+!--                   Spectral albedos especially for building wall surfaces
+                      IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )  THEN
+                         surf_usm_h(l)%aldir(m,ind_veg_wall) =                       &
+                                                   albedo_pars_f%pars_xy(1,j,i)
+                         surf_usm_h(l)%aldif(m,ind_veg_wall) =                       &
+                                                   albedo_pars_f%pars_xy(1,j,i)
+                      ENDIF
+                      IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill )  THEN
+                         surf_usm_h(l)%asdir(m,ind_veg_wall) =                       &
+                                                   albedo_pars_f%pars_xy(2,j,i)
+                         surf_usm_h(l)%asdif(m,ind_veg_wall) =                       &
+                                                   albedo_pars_f%pars_xy(2,j,i)
+                      ENDIF
+!
+!--                   Spectral albedos especially for building green surfaces
+                      IF ( albedo_pars_f%pars_xy(3,j,i) /= albedo_pars_f%fill )  THEN
+                         surf_usm_h(l)%aldir(m,ind_pav_green) =                      &
+                                                   albedo_pars_f%pars_xy(3,j,i)
+                         surf_usm_h(l)%aldif(m,ind_pav_green) =                      &
+                                                   albedo_pars_f%pars_xy(3,j,i)
+                      ENDIF
+                      IF ( albedo_pars_f%pars_xy(4,j,i) /= albedo_pars_f%fill )  THEN
+                         surf_usm_h(l)%asdir(m,ind_pav_green) =                      &
+                                                   albedo_pars_f%pars_xy(4,j,i)
+                         surf_usm_h(l)%asdif(m,ind_pav_green) =                      &
+                                                   albedo_pars_f%pars_xy(4,j,i)
+                      ENDIF
+!
+!--                   Spectral albedos especially for building window surfaces
+                      IF ( albedo_pars_f%pars_xy(5,j,i) /= albedo_pars_f%fill )  THEN
+                         surf_usm_h(l)%aldir(m,ind_wat_win) =                        &
+                                                   albedo_pars_f%pars_xy(5,j,i)
+                         surf_usm_h(l)%aldif(m,ind_wat_win) =                        &
+                                                   albedo_pars_f%pars_xy(5,j,i)
+                      ENDIF
+                      IF ( albedo_pars_f%pars_xy(6,j,i) /= albedo_pars_f%fill )  THEN
+                         surf_usm_h(l)%asdir(m,ind_wat_win) =                        &
+                                                   albedo_pars_f%pars_xy(6,j,i)
+                         surf_usm_h(l)%asdif(m,ind_wat_win) =                        &
+                                                   albedo_pars_f%pars_xy(6,j,i)
+                      ENDIF
+
+                   ENDDO
+                ENDIF
              ENDDO
-!
-!--          For urban surface only if albedo has not been already initialized
-!--          in the urban-surface model via the ASCII file.
-             IF ( .NOT. surf_usm_h%albedo_from_ascii )  THEN
-                DO  m = 1, surf_usm_h%ns
-                   i = surf_usm_h%i(m)
-                   j = surf_usm_h%j(m)
-!
-!--                Broadband albedos for wall/green/window surfaces
-                   DO  ind_type = 0, 2
-                      IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )&
-                         surf_usm_h%albedo(m,ind_type) =                       &
-                                             albedo_pars_f%pars_xy(0,j,i)
-                   ENDDO
-!
-!--                Spectral albedos especially for building wall surfaces
-                   IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )  THEN
-                      surf_usm_h%aldir(m,ind_veg_wall) =                       &
-                                                albedo_pars_f%pars_xy(1,j,i)
-                      surf_usm_h%aldif(m,ind_veg_wall) =                       &
-                                                albedo_pars_f%pars_xy(1,j,i)
-                   ENDIF
-                   IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill )  THEN
-                      surf_usm_h%asdir(m,ind_veg_wall) =                       &
-                                                albedo_pars_f%pars_xy(2,j,i)
-                      surf_usm_h%asdif(m,ind_veg_wall) =                       &
-                                                albedo_pars_f%pars_xy(2,j,i)
-                   ENDIF
-!
-!--                Spectral albedos especially for building green surfaces
-                   IF ( albedo_pars_f%pars_xy(3,j,i) /= albedo_pars_f%fill )  THEN
-                      surf_usm_h%aldir(m,ind_pav_green) =                      &
-                                                albedo_pars_f%pars_xy(3,j,i)
-                      surf_usm_h%aldif(m,ind_pav_green) =                      &
-                                                albedo_pars_f%pars_xy(3,j,i)
-                   ENDIF
-                   IF ( albedo_pars_f%pars_xy(4,j,i) /= albedo_pars_f%fill )  THEN
-                      surf_usm_h%asdir(m,ind_pav_green) =                      &
-                                                albedo_pars_f%pars_xy(4,j,i)
-                      surf_usm_h%asdif(m,ind_pav_green) =                      &
-                                                albedo_pars_f%pars_xy(4,j,i)
-                   ENDIF
-!
-!--                Spectral albedos especially for building window surfaces
-                   IF ( albedo_pars_f%pars_xy(5,j,i) /= albedo_pars_f%fill )  THEN
-                      surf_usm_h%aldir(m,ind_wat_win) =                        &
-                                                albedo_pars_f%pars_xy(5,j,i)
-                      surf_usm_h%aldif(m,ind_wat_win) =                        &
-                                                albedo_pars_f%pars_xy(5,j,i)
-                   ENDIF
-                   IF ( albedo_pars_f%pars_xy(6,j,i) /= albedo_pars_f%fill )  THEN
-                      surf_usm_h%asdir(m,ind_wat_win) =                        &
-                                                albedo_pars_f%pars_xy(6,j,i)
-                      surf_usm_h%asdif(m,ind_wat_win) =                        &
-                                                albedo_pars_f%pars_xy(6,j,i)
-                   ENDIF
-
-                ENDDO
-             ENDIF
 !
 !--          Vertical
@@ -2530,95 +2555,96 @@
 !--       to zero in order to take effect.
           IF ( building_surface_pars_f%from_file )  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)
-!
-!--             Iterate over surfaces in column, check height and orientation
-                DO  is = building_surface_pars_f%index_ji(1,j,i), &
-                         building_surface_pars_f%index_ji(2,j,i)
-                   IF ( building_surface_pars_f%coords(4,is) == -surf_usm_h%koff .AND. &
-                        building_surface_pars_f%coords(1,is) == k )  THEN
-
-                      IF ( building_surface_pars_f%pars(ind_s_alb_b_wall,is) /=      &
-                           building_surface_pars_f%fill )  THEN
-                         surf_usm_h%albedo(m,ind_veg_wall) =                         &
-                                  building_surface_pars_f%pars(ind_s_alb_b_wall,is)
-                         surf_usm_h%albedo_type(m,ind_veg_wall) = 0
+             DO l = 0, 1
+                DO  m = 1, surf_usm_h(l)%ns
+                   i = surf_usm_h(l)%i(m)
+                   j = surf_usm_h(l)%j(m)
+                   k = surf_usm_h(l)%k(m)
+   !
+   !--             Iterate over surfaces in column, check height and orientation
+                   DO  is = building_surface_pars_f%index_ji(1,j,i), &
+                            building_surface_pars_f%index_ji(2,j,i)
+                      IF ( building_surface_pars_f%coords(4,is) == -surf_usm_h(l)%koff .AND. &
+                           building_surface_pars_f%coords(1,is) == k )  THEN
+
+                         IF ( building_surface_pars_f%pars(ind_s_alb_b_wall,is) /=      &
+                              building_surface_pars_f%fill )  THEN
+                            surf_usm_h(l)%albedo(m,ind_veg_wall) =                         &
+                                     building_surface_pars_f%pars(ind_s_alb_b_wall,is)
+                            surf_usm_h(l)%albedo_type(m,ind_veg_wall) = 0
+                         ENDIF
+
+                         IF ( building_surface_pars_f%pars(ind_s_alb_l_wall,is) /=      &
+                              building_surface_pars_f%fill )  THEN
+                            surf_usm_h(l)%aldir(m,ind_veg_wall) =                          &
+                                     building_surface_pars_f%pars(ind_s_alb_l_wall,is)
+                            surf_usm_h(l)%aldif(m,ind_veg_wall) =                          &
+                                     building_surface_pars_f%pars(ind_s_alb_l_wall,is)
+                            surf_usm_h(l)%albedo_type(m,ind_veg_wall) = 0
+                         ENDIF
+
+                         IF ( building_surface_pars_f%pars(ind_s_alb_s_wall,is) /=      &
+                              building_surface_pars_f%fill )  THEN
+                            surf_usm_h(l)%asdir(m,ind_veg_wall) =                          &
+                                     building_surface_pars_f%pars(ind_s_alb_s_wall,is)
+                            surf_usm_h(l)%asdif(m,ind_veg_wall) =                          &
+                                     building_surface_pars_f%pars(ind_s_alb_s_wall,is)
+                            surf_usm_h(l)%albedo_type(m,ind_veg_wall) = 0
+                         ENDIF
+
+                         IF ( building_surface_pars_f%pars(ind_s_alb_b_win,is) /=       &
+                              building_surface_pars_f%fill )  THEN
+                            surf_usm_h(l)%albedo(m,ind_wat_win) =                          &
+                                     building_surface_pars_f%pars(ind_s_alb_b_win,is)
+                            surf_usm_h(l)%albedo_type(m,ind_wat_win) = 0
+                         ENDIF
+
+                         IF ( building_surface_pars_f%pars(ind_s_alb_l_win,is) /=       &
+                              building_surface_pars_f%fill )  THEN
+                            surf_usm_h(l)%aldir(m,ind_wat_win) =                           &
+                                     building_surface_pars_f%pars(ind_s_alb_l_win,is)
+                            surf_usm_h(l)%aldif(m,ind_wat_win) =                           &
+                                     building_surface_pars_f%pars(ind_s_alb_l_win,is)
+                            surf_usm_h(l)%albedo_type(m,ind_wat_win) = 0
+                         ENDIF
+
+                         IF ( building_surface_pars_f%pars(ind_s_alb_s_win,is) /=       &
+                              building_surface_pars_f%fill )  THEN
+                            surf_usm_h(l)%asdir(m,ind_wat_win) =                           &
+                                     building_surface_pars_f%pars(ind_s_alb_s_win,is)
+                            surf_usm_h(l)%asdif(m,ind_wat_win) =                           &
+                                     building_surface_pars_f%pars(ind_s_alb_s_win,is)
+                            surf_usm_h(l)%albedo_type(m,ind_wat_win) = 0
+                         ENDIF
+
+                         IF ( building_surface_pars_f%pars(ind_s_alb_b_green,is) /=     &
+                              building_surface_pars_f%fill )  THEN
+                            surf_usm_h(l)%albedo(m,ind_pav_green) =                        &
+                                     building_surface_pars_f%pars(ind_s_alb_b_green,is)
+                            surf_usm_h(l)%albedo_type(m,ind_pav_green) = 0
+                         ENDIF
+
+                         IF ( building_surface_pars_f%pars(ind_s_alb_l_green,is) /=     &
+                              building_surface_pars_f%fill )  THEN
+                            surf_usm_h(l)%aldir(m,ind_pav_green) =                         &
+                                     building_surface_pars_f%pars(ind_s_alb_l_green,is)
+                            surf_usm_h(l)%aldif(m,ind_pav_green) =                         &
+                                     building_surface_pars_f%pars(ind_s_alb_l_green,is)
+                            surf_usm_h(l)%albedo_type(m,ind_pav_green) = 0
+                         ENDIF
+
+                         IF ( building_surface_pars_f%pars(ind_s_alb_s_green,is) /=     &
+                              building_surface_pars_f%fill )  THEN
+                            surf_usm_h(l)%asdir(m,ind_pav_green) =                         &
+                                     building_surface_pars_f%pars(ind_s_alb_s_green,is)
+                            surf_usm_h(l)%asdif(m,ind_pav_green) =                         &
+                                     building_surface_pars_f%pars(ind_s_alb_s_green,is)
+                            surf_usm_h(l)%albedo_type(m,ind_pav_green) = 0
+                         ENDIF
+
+                         EXIT ! surface was found and processed
                       ENDIF
-
-                      IF ( building_surface_pars_f%pars(ind_s_alb_l_wall,is) /=      &
-                           building_surface_pars_f%fill )  THEN
-                         surf_usm_h%aldir(m,ind_veg_wall) =                          &
-                                  building_surface_pars_f%pars(ind_s_alb_l_wall,is)
-                         surf_usm_h%aldif(m,ind_veg_wall) =                          &
-                                  building_surface_pars_f%pars(ind_s_alb_l_wall,is)
-                         surf_usm_h%albedo_type(m,ind_veg_wall) = 0
-                      ENDIF
-
-                      IF ( building_surface_pars_f%pars(ind_s_alb_s_wall,is) /=      &
-                           building_surface_pars_f%fill )  THEN
-                         surf_usm_h%asdir(m,ind_veg_wall) =                          &
-                                  building_surface_pars_f%pars(ind_s_alb_s_wall,is)
-                         surf_usm_h%asdif(m,ind_veg_wall) =                          &
-                                  building_surface_pars_f%pars(ind_s_alb_s_wall,is)
-                         surf_usm_h%albedo_type(m,ind_veg_wall) = 0
-                      ENDIF
-
-                      IF ( building_surface_pars_f%pars(ind_s_alb_b_win,is) /=       &
-                           building_surface_pars_f%fill )  THEN
-                         surf_usm_h%albedo(m,ind_wat_win) =                          &
-                                  building_surface_pars_f%pars(ind_s_alb_b_win,is)
-                         surf_usm_h%albedo_type(m,ind_wat_win) = 0
-                      ENDIF
-
-                      IF ( building_surface_pars_f%pars(ind_s_alb_l_win,is) /=       &
-                           building_surface_pars_f%fill )  THEN
-                         surf_usm_h%aldir(m,ind_wat_win) =                           &
-                                  building_surface_pars_f%pars(ind_s_alb_l_win,is)
-                         surf_usm_h%aldif(m,ind_wat_win) =                           &
-                                  building_surface_pars_f%pars(ind_s_alb_l_win,is)
-                         surf_usm_h%albedo_type(m,ind_wat_win) = 0
-                      ENDIF
-
-                      IF ( building_surface_pars_f%pars(ind_s_alb_s_win,is) /=       &
-                           building_surface_pars_f%fill )  THEN
-                         surf_usm_h%asdir(m,ind_wat_win) =                           &
-                                  building_surface_pars_f%pars(ind_s_alb_s_win,is)
-                         surf_usm_h%asdif(m,ind_wat_win) =                           &
-                                  building_surface_pars_f%pars(ind_s_alb_s_win,is)
-                         surf_usm_h%albedo_type(m,ind_wat_win) = 0
-                      ENDIF
-
-                      IF ( building_surface_pars_f%pars(ind_s_alb_b_green,is) /=     &
-                           building_surface_pars_f%fill )  THEN
-                         surf_usm_h%albedo(m,ind_pav_green) =                        &
-                                  building_surface_pars_f%pars(ind_s_alb_b_green,is)
-                         surf_usm_h%albedo_type(m,ind_pav_green) = 0
-                      ENDIF
-
-                      IF ( building_surface_pars_f%pars(ind_s_alb_l_green,is) /=     &
-                           building_surface_pars_f%fill )  THEN
-                         surf_usm_h%aldir(m,ind_pav_green) =                         &
-                                  building_surface_pars_f%pars(ind_s_alb_l_green,is)
-                         surf_usm_h%aldif(m,ind_pav_green) =                         &
-                                  building_surface_pars_f%pars(ind_s_alb_l_green,is)
-                         surf_usm_h%albedo_type(m,ind_pav_green) = 0
-                      ENDIF
-
-                      IF ( building_surface_pars_f%pars(ind_s_alb_s_green,is) /=     &
-                           building_surface_pars_f%fill )  THEN
-                         surf_usm_h%asdir(m,ind_pav_green) =                         &
-                                  building_surface_pars_f%pars(ind_s_alb_s_green,is)
-                         surf_usm_h%asdif(m,ind_pav_green) =                         &
-                                  building_surface_pars_f%pars(ind_s_alb_s_green,is)
-                         surf_usm_h%albedo_type(m,ind_pav_green) = 0
-                      ENDIF
-
-                      EXIT ! surface was found and processed
-                   ENDIF
+                   ENDDO
                 ENDDO
              ENDDO
-
              DO  l = 0, 3
                 DO  m = 1, surf_usm_v(l)%ns
@@ -2729,6 +2755,8 @@
 !
 !--          Horizontally aligned natural and urban surfaces
-             CALL calc_albedo( surf_lsm_h )
-             CALL calc_albedo( surf_usm_h )
+             DO l = 0, 1
+                CALL calc_albedo( surf_lsm_h(l) )
+                CALL calc_albedo( surf_usm_h(l) )
+             ENDDO
 !
 !--          Vertically aligned natural and urban surfaces
@@ -2742,16 +2770,18 @@
 !--          Initialize sun-inclination independent spectral albedos
 !--          Horizontal surfaces
-             IF ( surf_lsm_h%ns > 0 )  THEN
-                surf_lsm_h%rrtm_aldir = surf_lsm_h%aldir
-                surf_lsm_h%rrtm_asdir = surf_lsm_h%asdir
-                surf_lsm_h%rrtm_aldif = surf_lsm_h%aldif
-                surf_lsm_h%rrtm_asdif = surf_lsm_h%asdif
-             ENDIF
-             IF ( surf_usm_h%ns > 0 )  THEN
-                surf_usm_h%rrtm_aldir = surf_usm_h%aldir
-                surf_usm_h%rrtm_asdir = surf_usm_h%asdir
-                surf_usm_h%rrtm_aldif = surf_usm_h%aldif
-                surf_usm_h%rrtm_asdif = surf_usm_h%asdif
-             ENDIF
+             DO  l = 0, 1
+                IF ( surf_lsm_h(l)%ns > 0 )  THEN
+                   surf_lsm_h(l)%rrtm_aldir = surf_lsm_h(l)%aldir
+                   surf_lsm_h(l)%rrtm_asdir = surf_lsm_h(l)%asdir
+                   surf_lsm_h(l)%rrtm_aldif = surf_lsm_h(l)%aldif
+                   surf_lsm_h(l)%rrtm_asdif = surf_lsm_h(l)%asdif
+                ENDIF
+                IF ( surf_usm_h(l)%ns > 0 )  THEN
+                   surf_usm_h(l)%rrtm_aldir = surf_usm_h(l)%aldir
+                   surf_usm_h(l)%rrtm_asdir = surf_usm_h(l)%asdir
+                   surf_usm_h(l)%rrtm_aldif = surf_usm_h(l)%aldif
+                   surf_usm_h(l)%rrtm_asdif = surf_usm_h(l)%asdif
+                ENDIF
+             ENDDO
 !
 !--          Vertical surfaces
@@ -3287,8 +3317,10 @@
 !--    First, horizontal surfaces
        horizontal = .TRUE.
-       surf => surf_lsm_h
-       CALL radiation_external_surf
-       surf => surf_usm_h
-       CALL radiation_external_surf
+       DO  l = 0, 1
+          surf => surf_lsm_h(l)
+          CALL radiation_external_surf
+          surf => surf_usm_h(l)
+          CALL radiation_external_surf
+       ENDDO
        horizontal = .FALSE.
 !
@@ -3567,8 +3599,10 @@
 !--    First, horizontal surfaces
        horizontal = .TRUE.
-       surf => surf_lsm_h
-       CALL radiation_clearsky_surf
-       surf => surf_usm_h
-       CALL radiation_clearsky_surf
+       DO  l = 0, 1
+          surf => surf_lsm_h(l)
+          CALL radiation_clearsky_surf
+          surf => surf_usm_h(l)
+          CALL radiation_clearsky_surf
+       ENDDO
        horizontal = .FALSE.
 !
@@ -3752,8 +3786,10 @@
 !--    First, horizontal surfaces
        horizontal = .TRUE.
-       surf => surf_lsm_h
-       CALL radiation_constant_surf
-       surf => surf_usm_h
-       CALL radiation_constant_surf
+       DO  l = 0, 1
+          surf => surf_lsm_h(l)
+          CALL radiation_constant_surf
+          surf => surf_usm_h(l)
+          CALL radiation_constant_surf
+       ENDDO
        horizontal = .FALSE.
 !
@@ -4125,6 +4161,8 @@
 !
 !--       Horizontally aligned default, natural and urban surfaces
-          CALL calc_albedo( surf_lsm_h    )
-          CALL calc_albedo( surf_usm_h    )
+          DO  l = 0, 1
+             CALL calc_albedo( surf_lsm_h(l) )
+             CALL calc_albedo( surf_usm_h(l) )
+          ENDDO
 !
 !--       Vertically aligned default, natural and urban surfaces
@@ -4472,9 +4510,9 @@
 !
 !--             Prepare profiles of temperature and H2O volume mixing ratio
-                DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
-                   rrtm_tlev(0,nzb+1) = surf_lsm_h%pt_surface(m) * exner(nzb)
+                DO  m = surf_lsm_h(0)%start_index(j,i), surf_lsm_h(0)%end_index(j,i)
+                   rrtm_tlev(0,nzb+1) = surf_lsm_h(0)%pt_surface(m) * exner(nzb)
                 ENDDO
-                DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
-                   rrtm_tlev(0,nzb+1) = surf_usm_h%pt_surface(m) * exner(nzb)
+                DO  m = surf_usm_h(0)%start_index(j,i), surf_usm_h(0)%end_index(j,i)
+                   rrtm_tlev(0,nzb+1) = surf_usm_h(0)%pt_surface(m) * exner(nzb)
                 ENDDO
 
@@ -4608,6 +4646,6 @@
 !--             surface element on RRTMG-shaped array.
 !--             Please note, as RRTMG is a single column model, surface attributes
-!--             are only obtained from horizontally aligned surfaces (for
-!--             simplicity). Taking surface attributes from horizontal and
+!--             are only obtained from upward facing horizontally aligned surfaces
+!--             (for simplicity). Taking surface attributes from horizontal and
 !--             vertical walls would lead to multiple solutions.
 !--             Moreover, for natural- and urban-type surfaces, several surface
@@ -4615,21 +4653,21 @@
 !--             To obtain bulk parameters, apply a weighted average for these
 !--             surfaces.
-                DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
-                   rrtm_emis = surf_lsm_h%frac(m,ind_veg_wall)  *              &
-                               surf_lsm_h%emissivity(m,ind_veg_wall)  +        &
-                               surf_lsm_h%frac(m,ind_pav_green) *              &
-                               surf_lsm_h%emissivity(m,ind_pav_green) +        &
-                               surf_lsm_h%frac(m,ind_wat_win)   *              &
-                               surf_lsm_h%emissivity(m,ind_wat_win)
-                   rrtm_tsfc = surf_lsm_h%pt_surface(m) * exner(nzb)
+                DO  m = surf_lsm_h(0)%start_index(j,i), surf_lsm_h(0)%end_index(j,i)
+                   rrtm_emis = surf_lsm_h(0)%frac(m,ind_veg_wall)  *              &
+                               surf_lsm_h(0)%emissivity(m,ind_veg_wall)  +        &
+                               surf_lsm_h(0)%frac(m,ind_pav_green) *              &
+                               surf_lsm_h(0)%emissivity(m,ind_pav_green) +        &
+                               surf_lsm_h(0)%frac(m,ind_wat_win)   *              &
+                               surf_lsm_h(0)%emissivity(m,ind_wat_win)
+                   rrtm_tsfc = surf_lsm_h(0)%pt_surface(m) * exner(nzb)
                 ENDDO
-                DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
-                   rrtm_emis = surf_usm_h%frac(m,ind_veg_wall)  *              &
-                               surf_usm_h%emissivity(m,ind_veg_wall)  +        &
-                               surf_usm_h%frac(m,ind_pav_green) *              &
-                               surf_usm_h%emissivity(m,ind_pav_green) +        &
-                               surf_usm_h%frac(m,ind_wat_win)   *              &
-                               surf_usm_h%emissivity(m,ind_wat_win)
-                   rrtm_tsfc = surf_usm_h%pt_surface(m) * exner(nzb)
+                DO  m = surf_usm_h(0)%start_index(j,i), surf_usm_h(0)%end_index(j,i)
+                   rrtm_emis = surf_usm_h(0)%frac(m,ind_veg_wall)  *              &
+                               surf_usm_h(0)%emissivity(m,ind_veg_wall)  +        &
+                               surf_usm_h(0)%frac(m,ind_pav_green) *              &
+                               surf_usm_h(0)%emissivity(m,ind_pav_green) +        &
+                               surf_usm_h(0)%frac(m,ind_wat_win)   *              &
+                               surf_usm_h(0)%emissivity(m,ind_wat_win)
+                   rrtm_tsfc = surf_usm_h(0)%pt_surface(m) * exner(nzb)
                 ENDDO
 !
@@ -4713,15 +4751,17 @@
 !--                onto respective surface elements
 !--                Horizontal surfaces
-                   DO  m = surf_lsm_h%start_index(j,i),                        &
-                           surf_lsm_h%end_index(j,i)
-                      surf_lsm_h%rad_lw_in(m)           = rrtm_lwdflx(0,k_topo)
-                      surf_lsm_h%rad_lw_out(m)          = rrtm_lwuflx(0,k_topo)
-                      surf_lsm_h%rad_lw_out_change_0(m) = rrtm_lwuflx_dt(0,k_topo)
-                   ENDDO
-                   DO  m = surf_usm_h%start_index(j,i),                        &
-                           surf_usm_h%end_index(j,i)
-                      surf_usm_h%rad_lw_in(m)           = rrtm_lwdflx(0,k_topo)
-                      surf_usm_h%rad_lw_out(m)          = rrtm_lwuflx(0,k_topo)
-                      surf_usm_h%rad_lw_out_change_0(m) = rrtm_lwuflx_dt(0,k_topo)
+                   DO l = 0, 1
+                      DO  m = surf_lsm_h(l)%start_index(j,i),                        &
+                              surf_lsm_h(l)%end_index(j,i)
+                         surf_lsm_h(l)%rad_lw_in(m)           = rrtm_lwdflx(0,k_topo)
+                         surf_lsm_h(l)%rad_lw_out(m)          = rrtm_lwuflx(0,k_topo)
+                         surf_lsm_h(l)%rad_lw_out_change_0(m) = rrtm_lwuflx_dt(0,k_topo)
+                      ENDDO
+                      DO  m = surf_usm_h(l)%start_index(j,i),                        &
+                              surf_usm_h(l)%end_index(j,i)
+                         surf_usm_h(l)%rad_lw_in(m)           = rrtm_lwdflx(0,k_topo)
+                         surf_usm_h(l)%rad_lw_out(m)          = rrtm_lwuflx(0,k_topo)
+                         surf_usm_h(l)%rad_lw_out_change_0(m) = rrtm_lwuflx_dt(0,k_topo)
+                      ENDDO
                    ENDDO
 !
@@ -4751,29 +4791,28 @@
 !--                Get albedo for direct/diffusive long/shortwave radiation at
 !--                current (y,x)-location from surface variables.
-!--                Only obtain it from horizontal surfaces, as RRTMG is a single
-!--                column model
-!--                (Please note, only one loop will entered, controlled by
-!--                start-end index.)
-                   DO  m = surf_lsm_h%start_index(j,i),                        &
-                           surf_lsm_h%end_index(j,i)
-                      rrtm_asdir(1)  = SUM( surf_lsm_h%frac(m,:) *             &
-                                            surf_lsm_h%rrtm_asdir(m,:) )
-                      rrtm_asdif(1)  = SUM( surf_lsm_h%frac(m,:) *             &
-                                            surf_lsm_h%rrtm_asdif(m,:) )
-                      rrtm_aldir(1)  = SUM( surf_lsm_h%frac(m,:) *             &
-                                            surf_lsm_h%rrtm_aldir(m,:) )
-                      rrtm_aldif(1)  = SUM( surf_lsm_h%frac(m,:) *             &
-                                            surf_lsm_h%rrtm_aldif(m,:) )
+!--                Only obtain it from upward facing horizontal surfaces,
+!--                as RRTMG is a single column model. (Please note, only
+!--                one loop will entered, controlled by start-end index.)
+                   DO  m = surf_lsm_h(0)%start_index(j,i),                        &
+                           surf_lsm_h(0)%end_index(j,i)
+                      rrtm_asdir(1)  = SUM( surf_lsm_h(0)%frac(m,:) *             &
+                                            surf_lsm_h(0)%rrtm_asdir(m,:) )
+                      rrtm_asdif(1)  = SUM( surf_lsm_h(0)%frac(m,:) *             &
+                                            surf_lsm_h(0)%rrtm_asdif(m,:) )
+                      rrtm_aldir(1)  = SUM( surf_lsm_h(0)%frac(m,:) *             &
+                                            surf_lsm_h(0)%rrtm_aldir(m,:) )
+                      rrtm_aldif(1)  = SUM( surf_lsm_h(0)%frac(m,:) *             &
+                                            surf_lsm_h(0)%rrtm_aldif(m,:) )
                    ENDDO
-                   DO  m = surf_usm_h%start_index(j,i),                        &
-                           surf_usm_h%end_index(j,i)
-                      rrtm_asdir(1)  = SUM( surf_usm_h%frac(m,:) *             &
-                                            surf_usm_h%rrtm_asdir(m,:) )
-                      rrtm_asdif(1)  = SUM( surf_usm_h%frac(m,:) *             &
-                                            surf_usm_h%rrtm_asdif(m,:) )
-                      rrtm_aldir(1)  = SUM( surf_usm_h%frac(m,:) *             &
-                                            surf_usm_h%rrtm_aldir(m,:) )
-                      rrtm_aldif(1)  = SUM( surf_usm_h%frac(m,:) *             &
-                                            surf_usm_h%rrtm_aldif(m,:) )
+                   DO  m = surf_usm_h(0)%start_index(j,i),                        &
+                           surf_usm_h(0)%end_index(j,i)
+                      rrtm_asdir(1)  = SUM( surf_usm_h(0)%frac(m,:) *             &
+                                            surf_usm_h(0)%rrtm_asdir(m,:) )
+                      rrtm_asdif(1)  = SUM( surf_usm_h(0)%frac(m,:) *             &
+                                            surf_usm_h(0)%rrtm_asdif(m,:) )
+                      rrtm_aldir(1)  = SUM( surf_usm_h(0)%frac(m,:) *             &
+                                            surf_usm_h(0)%rrtm_aldir(m,:) )
+                      rrtm_aldif(1)  = SUM( surf_usm_h(0)%frac(m,:) *             &
+                                            surf_usm_h(0)%rrtm_aldif(m,:) )
                    ENDDO
 !
@@ -4876,13 +4915,15 @@
 !--                Save surface radiative fluxes onto respective surface elements
 !--                Horizontal surfaces
-                   DO  m = surf_lsm_h%start_index(j,i),                        &
-                           surf_lsm_h%end_index(j,i)
-                      surf_lsm_h%rad_sw_in(m)     = rrtm_swdflx(0,k_topo)
-                      surf_lsm_h%rad_sw_out(m)    = rrtm_swuflx(0,k_topo)
-                   ENDDO
-                   DO  m = surf_usm_h%start_index(j,i),                        &
-                           surf_usm_h%end_index(j,i)
-                      surf_usm_h%rad_sw_in(m)     = rrtm_swdflx(0,k_topo)
-                      surf_usm_h%rad_sw_out(m)    = rrtm_swuflx(0,k_topo)
+                   DO l = 0, 1
+                      DO  m = surf_lsm_h(l)%start_index(j,i),                        &
+                              surf_lsm_h(l)%end_index(j,i)
+                         surf_lsm_h(l)%rad_sw_in(m)     = rrtm_swdflx(0,k_topo)
+                         surf_lsm_h(l)%rad_sw_out(m)    = rrtm_swuflx(0,k_topo)
+                      ENDDO
+                      DO  m = surf_usm_h(l)%start_index(j,i),                        &
+                              surf_usm_h(l)%end_index(j,i)
+                         surf_usm_h(l)%rad_sw_in(m)     = rrtm_swdflx(0,k_topo)
+                         surf_usm_h(l)%rad_sw_out(m)    = rrtm_swuflx(0,k_topo)
+                      ENDDO
                    ENDDO
 !
@@ -4908,17 +4949,19 @@
                    rad_sw_in  = 0.0_wp
                    rad_sw_out = 0.0_wp
-!--             !!!!!!!! ATTENSION !!!!!!!!!!!!!!!
+!--             !!!!!!!! ATTENTION !!!!!!!!!!!!!!!
 !--             Surface radiative fluxes should be also set to zero here
 !--                Save surface radiative fluxes onto respective surface elements
 !--                Horizontal surfaces
-                   DO  m = surf_lsm_h%start_index(j,i),                        &
-                           surf_lsm_h%end_index(j,i)
-                      surf_lsm_h%rad_sw_in(m)     = 0.0_wp
-                      surf_lsm_h%rad_sw_out(m)    = 0.0_wp
-                   ENDDO
-                   DO  m = surf_usm_h%start_index(j,i),                        &
-                           surf_usm_h%end_index(j,i)
-                      surf_usm_h%rad_sw_in(m)     = 0.0_wp
-                      surf_usm_h%rad_sw_out(m)    = 0.0_wp
+                   DO l = 0, 1
+                      DO  m = surf_lsm_h(l)%start_index(j,i),                  &
+                              surf_lsm_h(l)%end_index(j,i)
+                         surf_lsm_h(l)%rad_sw_in(m)     = 0.0_wp
+                         surf_lsm_h(l)%rad_sw_out(m)    = 0.0_wp
+                      ENDDO
+                      DO  m = surf_usm_h(l)%start_index(j,i),                  &
+                              surf_usm_h(l)%end_index(j,i)
+                         surf_usm_h(l)%rad_sw_in(m)     = 0.0_wp
+                         surf_usm_h(l)%rad_sw_out(m)    = 0.0_wp
+                      ENDDO
                    ENDDO
 !
@@ -4949,15 +4992,17 @@
        IF (  .NOT.  radiation_interactions  ) THEN
 !--       First, for horizontal surfaces
-          DO  m = 1, surf_lsm_h%ns
-             surf_lsm_h%rad_net(m) = surf_lsm_h%rad_sw_in(m)                   &
-                                   - surf_lsm_h%rad_sw_out(m)                  &
-                                   + surf_lsm_h%rad_lw_in(m)                   &
-                                   - surf_lsm_h%rad_lw_out(m)
-          ENDDO
-          DO  m = 1, surf_usm_h%ns
-             surf_usm_h%rad_net(m) = surf_usm_h%rad_sw_in(m)                   &
-                                   - surf_usm_h%rad_sw_out(m)                  &
-                                   + surf_usm_h%rad_lw_in(m)                   &
-                                   - surf_usm_h%rad_lw_out(m)
+          DO l = 0, 1
+             DO  m = 1, surf_lsm_h(l)%ns
+                surf_lsm_h(l)%rad_net(m) = surf_lsm_h(l)%rad_sw_in(m)                &
+                                      - surf_lsm_h(l)%rad_sw_out(m)                  &
+                                      + surf_lsm_h(l)%rad_lw_in(m)                   &
+                                      - surf_lsm_h(l)%rad_lw_out(m)
+             ENDDO
+             DO  m = 1, surf_usm_h(l)%ns
+                surf_usm_h(l)%rad_net(m) = surf_usm_h(l)%rad_sw_in(m)                &
+                                      - surf_usm_h(l)%rad_sw_out(m)                  &
+                                      + surf_usm_h(l)%rad_lw_in(m)                   &
+                                      - surf_usm_h(l)%rad_lw_out(m)
+             ENDDO
           ENDDO
 !
@@ -5986,40 +6031,30 @@
 !
 !--        Horizontal walls
-!--        urban
-           DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
-              surfoutll(mm) = SUM ( surf_usm_h%frac(m,:) *                  &
-                                    surf_usm_h%emissivity(m,:) )            &
-                                  * sigma_sb                                &
-                                  * surf_usm_h%pt_surface(m)**4
-              albedo_surf(mm) = SUM ( surf_usm_h%frac(m,:) *                &
-                                      surf_usm_h%albedo(m,:) )
-              emiss_surf(mm)  = SUM ( surf_usm_h%frac(m,:) *                &
-                                      surf_usm_h%emissivity(m,:) )
-              mm = mm + 1
-           ENDDO
-!
-!--        land
-           DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
-              surfoutll(mm) = SUM ( surf_lsm_h%frac(m,:) *                  &
-                                    surf_lsm_h%emissivity(m,:) )            &
-                                  * sigma_sb                                &
-                                  * surf_lsm_h%pt_surface(m)**4
-              albedo_surf(mm) = SUM ( surf_lsm_h%frac(m,:) *                &
-                                      surf_lsm_h%albedo(m,:) )
-              emiss_surf(mm)  = SUM ( surf_lsm_h%frac(m,:) *                &
-                                      surf_lsm_h%emissivity(m,:) )
-              mm = mm + 1
-           ENDDO
-!
-!--        Downward facing default surfaces
-!--        TODO: Downward-facing surfaces are only available as default surfaces which lack
-!--        temperature and heat balance, therefore they have to be simulated as neutral to
-!--        radiation, that means having sw/lw reflectivity=1, so they neither absorb nor emit
-!--        radiation and they have zero net heat flux.
-           DO  m = surf_def_h(1)%start_index(j,i), surf_def_h(1)%end_index(j,i)
-              surfoutll(mm) = 0.0_wp
-              albedo_surf(mm) = 1.0_wp
-              emiss_surf(mm)  = 0.0_wp
-              mm = mm + 1
+           DO l = 0, 1
+!--           urban
+              DO  m = surf_usm_h(l)%start_index(j,i), surf_usm_h(l)%end_index(j,i)
+                 surfoutll(mm) = SUM ( surf_usm_h(l)%frac(m,:) *                  &
+                                       surf_usm_h(l)%emissivity(m,:) )            &
+                                     * sigma_sb                                   &
+                                     * surf_usm_h(l)%pt_surface(m)**4
+                 albedo_surf(mm) = SUM ( surf_usm_h(l)%frac(m,:) *                &
+                                         surf_usm_h(l)%albedo(m,:) )
+                 emiss_surf(mm)  = SUM ( surf_usm_h(l)%frac(m,:) *                &
+                                         surf_usm_h(l)%emissivity(m,:) )
+                 mm = mm + 1
+              ENDDO
+!
+!--           land
+              DO  m = surf_lsm_h(l)%start_index(j,i), surf_lsm_h(l)%end_index(j,i)
+                 surfoutll(mm) = SUM ( surf_lsm_h(l)%frac(m,:) *                  &
+                                       surf_lsm_h(l)%emissivity(m,:) )            &
+                                     * sigma_sb                                   &
+                                     * surf_lsm_h(l)%pt_surface(m)**4
+                 albedo_surf(mm) = SUM ( surf_lsm_h(l)%frac(m,:) *                &
+                                         surf_lsm_h(l)%albedo(m,:) )
+                 emiss_surf(mm)  = SUM ( surf_lsm_h(l)%frac(m,:) *                &
+                                         surf_lsm_h(l)%emissivity(m,:) )
+                 mm = mm + 1
+              ENDDO
            ENDDO
 !
@@ -6428,50 +6463,43 @@
 !
 !--        Horizontal walls
-!--        urban
-           DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
-              surf_usm_h%rad_sw_in(m)  = surfinsw(mm)
-              surf_usm_h%rad_sw_out(m) = surfoutsw(mm)
-              surf_usm_h%rad_sw_dir(m) = surfinswdir(mm)
-              surf_usm_h%rad_sw_dif(m) = surfinswdif(mm)
-              surf_usm_h%rad_sw_ref(m) = surfinsw(mm) - surfinswdir(mm) -        &
-                                         surfinswdif(mm)
-              surf_usm_h%rad_sw_res(m) = surfins(mm)
-              surf_usm_h%rad_lw_in(m)  = surfinlw(mm)
-              surf_usm_h%rad_lw_out(m) = surfoutlw(mm)
-              surf_usm_h%rad_net(m)    = surfinsw(mm) - surfoutsw(mm) +          &
-                                         surfinlw(mm) - surfoutlw(mm)
-              surf_usm_h%rad_net_l(m)  = surf_usm_h%rad_net(m)
-              surf_usm_h%rad_lw_dif(m) = surfinlwdif(mm)
-              surf_usm_h%rad_lw_ref(m) = surfinlw(mm) - surfinlwdif(mm)
-              surf_usm_h%rad_lw_res(m) = surfinl(mm)
-              mm = mm + 1
-           ENDDO
-!
-!--        land
-           DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
-              surf_lsm_h%rad_sw_in(m)  = surfinsw(mm)
-              surf_lsm_h%rad_sw_out(m) = surfoutsw(mm)
-              surf_lsm_h%rad_sw_dir(m) = surfinswdir(mm)
-              surf_lsm_h%rad_sw_dif(m) = surfinswdif(mm)
-              surf_lsm_h%rad_sw_ref(m) = surfinsw(mm) - surfinswdir(mm) -        &
+           DO l = 0, 1
+!--           urban
+              DO  m = surf_usm_h(l)%start_index(j,i), surf_usm_h(l)%end_index(j,i)
+                 surf_usm_h(l)%rad_sw_in(m)  = surfinsw(mm)
+                 surf_usm_h(l)%rad_sw_out(m) = surfoutsw(mm)
+                 surf_usm_h(l)%rad_sw_dir(m) = surfinswdir(mm)
+                 surf_usm_h(l)%rad_sw_dif(m) = surfinswdif(mm)
+                 surf_usm_h(l)%rad_sw_ref(m) = surfinsw(mm) - surfinswdir(mm) -        &
                                             surfinswdif(mm)
-              surf_lsm_h%rad_sw_res(m) = surfins(mm)
-              surf_lsm_h%rad_lw_in(m)  = surfinlw(mm)
-              surf_lsm_h%rad_lw_out(m) = surfoutlw(mm)
-              surf_lsm_h%rad_net(m)    = surfinsw(mm) - surfoutsw(mm) +          &
-                                         surfinlw(mm) - surfoutlw(mm)
-              surf_lsm_h%rad_lw_dif(m) = surfinlwdif(mm)
-              surf_lsm_h%rad_lw_ref(m) = surfinlw(mm) - surfinlwdif(mm)
-              surf_lsm_h%rad_lw_res(m) = surfinl(mm)
-              mm = mm + 1
-           ENDDO
-!
-!--        Downward facing default surfaces
-!--        TODO: Downward-facing surfaces are only available as default surfaces which lack
-!--        temperature and heat balance, therefore they have to be simulated as neutral to
-!--        radiation, that means having sw/lw reflectivity=1, so they neither absorb nor emit
-!--        radiation and they have zero net heat flux.
-           DO  m = surf_def_h(1)%start_index(j,i), surf_def_h(1)%end_index(j,i)
-              mm = mm + 1
+                 surf_usm_h(l)%rad_sw_res(m) = surfins(mm)
+                 surf_usm_h(l)%rad_lw_in(m)  = surfinlw(mm)
+                 surf_usm_h(l)%rad_lw_out(m) = surfoutlw(mm)
+                 surf_usm_h(l)%rad_net(m)    = surfinsw(mm) - surfoutsw(mm) +          &
+                                            surfinlw(mm) - surfoutlw(mm)
+                 surf_usm_h(l)%rad_net_l(m)  = surf_usm_h(l)%rad_net(m)
+                 surf_usm_h(l)%rad_lw_dif(m) = surfinlwdif(mm)
+                 surf_usm_h(l)%rad_lw_ref(m) = surfinlw(mm) - surfinlwdif(mm)
+                 surf_usm_h(l)%rad_lw_res(m) = surfinl(mm)
+                 mm = mm + 1
+              ENDDO
+!
+!--           land
+              DO  m = surf_lsm_h(l)%start_index(j,i), surf_lsm_h(l)%end_index(j,i)
+                 surf_lsm_h(l)%rad_sw_in(m)  = surfinsw(mm)
+                 surf_lsm_h(l)%rad_sw_out(m) = surfoutsw(mm)
+                 surf_lsm_h(l)%rad_sw_dir(m) = surfinswdir(mm)
+                 surf_lsm_h(l)%rad_sw_dif(m) = surfinswdif(mm)
+                 surf_lsm_h(l)%rad_sw_ref(m) = surfinsw(mm) - surfinswdir(mm) -        &
+                                               surfinswdif(mm)
+                 surf_lsm_h(l)%rad_sw_res(m) = surfins(mm)
+                 surf_lsm_h(l)%rad_lw_in(m)  = surfinlw(mm)
+                 surf_lsm_h(l)%rad_lw_out(m) = surfoutlw(mm)
+                 surf_lsm_h(l)%rad_net(m)    = surfinsw(mm) - surfoutsw(mm) +          &
+                                            surfinlw(mm) - surfoutlw(mm)
+                 surf_lsm_h(l)%rad_lw_dif(m) = surfinlwdif(mm)
+                 surf_lsm_h(l)%rad_lw_ref(m) = surfinlw(mm) - surfinlwdif(mm)
+                 surf_lsm_h(l)%rad_lw_res(m) = surfinl(mm)
+                 mm = mm + 1
+              ENDDO
            ENDDO
 !
@@ -6538,17 +6566,18 @@
      ENDDO
 
-     DO  m = 1, surf_usm_h%ns
-        surf_usm_h%surfhf(m) = surf_usm_h%rad_sw_in(m)  +                   &
-                               surf_usm_h%rad_lw_in(m)  -                   &
-                               surf_usm_h%rad_sw_out(m) -                   &
-                               surf_usm_h%rad_lw_out(m)
+     DO l = 0, 1
+        DO  m = 1, surf_usm_h(l)%ns
+           surf_usm_h(l)%surfhf(m) = surf_usm_h(l)%rad_sw_in(m)  +          &
+                                  surf_usm_h(l)%rad_lw_in(m)  -             &
+                                  surf_usm_h(l)%rad_sw_out(m) -             &
+                                  surf_usm_h(l)%rad_lw_out(m)
+        ENDDO
+        DO  m = 1, surf_lsm_h(l)%ns
+           surf_lsm_h(l)%surfhf(m) = surf_lsm_h(l)%rad_sw_in(m)  +          &
+                                  surf_lsm_h(l)%rad_lw_in(m)  -             &
+                                  surf_lsm_h(l)%rad_sw_out(m) -             &
+                                  surf_lsm_h(l)%rad_lw_out(m)
+        ENDDO
      ENDDO
-     DO  m = 1, surf_lsm_h%ns
-        surf_lsm_h%surfhf(m) = surf_lsm_h%rad_sw_in(m)  +                   &
-                               surf_lsm_h%rad_lw_in(m)  -                   &
-                               surf_lsm_h%rad_sw_out(m) -                   &
-                               surf_lsm_h%rad_lw_out(m)
-     ENDDO
-
      DO  l = 0, 3
 !--     urban
@@ -7020,10 +7049,12 @@
 !--    Number of horizontal surfaces including land- and roof surfaces in both USM and LSM. Note that
 !--    All horizontal surface elements are already counted in surface_mod.
-       nsurfl    = surf_usm_h%ns + surf_lsm_h%ns + surf_def_h(1)%ns
+       DO l = 0, 1
+          nsurfl = nsurfl + surf_usm_h(l)%ns + surf_lsm_h(l)%ns
+       ENDDO
 !
 !--    Number of vertical surfaces in both USM and LSM. Note that all vertical surface elements are
 !--    already counted in surface_mod.
-       DO  i = 0,3
-          nsurfl = nsurfl + surf_usm_v(i)%ns + surf_lsm_v(i)%ns
+       DO  l = 0,3
+          nsurfl = nsurfl + surf_usm_v(l)%ns + surf_lsm_v(l)%ns
        ENDDO
 
@@ -7144,82 +7175,47 @@
 !
 !--           Horizontal surfaces
-              DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
-                 k = surf_usm_h%k(m)
-                 isurf = isurf + 1
-                 surfl(:,isurf) = (/iup,k,j,i/)
+              DO l = 0, 1
+                 DO  m = surf_usm_h(l)%start_index(j,i), surf_usm_h(l)%end_index(j,i)
+                    k = surf_usm_h(l)%k(m)
+                    isurf = isurf + 1
+                    surfl(:,isurf) = (/dirint_h(l),k,j,i/)
+                 ENDDO
+                 DO  m = surf_lsm_h(l)%start_index(j,i), surf_lsm_h(l)%end_index(j,i)
+                    k = surf_lsm_h(l)%k(m)
+                    isurf = isurf + 1
+                    surfl(:,isurf) = (/dirint_h(l),k,j,i/)
+                 ENDDO
               ENDDO
-
-              DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
-                 k = surf_lsm_h%k(m)
-                 isurf = isurf + 1
-                 surfl(:,isurf) = (/iup,k,j,i/)
+!
+!--           Vertical surfaces
+              DO l = 0, 3
+                 DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
+                    k = surf_usm_v(l)%k(m)
+                    isurf = isurf + 1
+                    surfl(:,isurf) = (/dirint_v(l),k,j,i/)
+                 ENDDO
+                 DO  m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i)
+                    k = surf_lsm_v(l)%k(m)
+                    isurf = isurf + 1
+                    surfl(:,isurf) = (/dirint_v(l),k,j,i/)
+                 ENDDO
               ENDDO
 
-              DO  m = surf_def_h(1)%start_index(j,i), surf_def_h(1)%end_index(j,i)
-                 k = surf_def_h(1)%k(m)
-                 isurf = isurf + 1
-                 surfl(:,isurf) = (/idown,k,j,i/)
-              ENDDO
-!
-!--           Vertical surfaces
-              l = 0
-              DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
-                 k = surf_usm_v(l)%k(m)
-                 isurf = isurf + 1
-                 surfl(:,isurf) = (/inorth,k,j,i/)
-              ENDDO
-              DO  m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i)
-                 k = surf_lsm_v(l)%k(m)
-                 isurf = isurf + 1
-                 surfl(:,isurf) = (/inorth,k,j,i/)
-              ENDDO
-
-              l = 1
-              DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
-                 k = surf_usm_v(l)%k(m)
-                 isurf = isurf + 1
-                 surfl(:,isurf) = (/isouth,k,j,i/)
-              ENDDO
-              DO  m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i)
-                 k = surf_lsm_v(l)%k(m)
-                 isurf = isurf + 1
-                 surfl(:,isurf) = (/isouth,k,j,i/)
-              ENDDO
-
-              l = 2
-              DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
-                 k = surf_usm_v(l)%k(m)
-                 isurf = isurf + 1
-                 surfl(:,isurf) = (/ieast,k,j,i/)
-              ENDDO
-              DO  m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i)
-                 k = surf_lsm_v(l)%k(m)
-                 isurf = isurf + 1
-                 surfl(:,isurf) = (/ieast,k,j,i/)
-              ENDDO
-
-              l = 3
-              DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
-                 k = surf_usm_v(l)%k(m)
-                 isurf = isurf + 1
-                 surfl(:,isurf) = (/iwest,k,j,i/)
-              ENDDO
-              DO  m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i)
-                 k = surf_lsm_v(l)%k(m)
-                 isurf = isurf + 1
-                 surfl(:,isurf) = (/iwest,k,j,i/)
-              ENDDO
            ENDDO
        ENDDO
 !
-!--    Add local MRT boxes for specified number of levels
+!--    Add local MRT boxes for the specified number of levels
+!--    !!!! NEEDS TO RETHINK AGAIN - with full 3D structure, only the one
+!--    !!!! of the upward faced horizontal surfaces should be taken (the lowest one = ground?)
+!--    !!!! mrt_nlevels number of air grid boxes might not be available
+!--    !!!! in case of overhanging structures!
        nmrtbl = 0
        IF ( mrt_nlevels > 0 )  THEN
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
+                DO  m = surf_usm_h(0)%start_index(j,i), surf_usm_h(0)%end_index(j,i)
 !
 !--                Skip roof if requested
-                   IF ( mrt_skip_roof  .AND.  surf_usm_h%isroof_surf(m) )  CYCLE
+                   IF ( mrt_skip_roof  .AND.  surf_usm_h(0)%isroof_surf(m) )  CYCLE
 !
 !--                Cycle over specified no of levels
@@ -7228,5 +7224,5 @@
 !
 !--             Dtto for LSM
-                DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
+                DO  m = surf_lsm_h(0)%start_index(j,i), surf_lsm_h(0)%end_index(j,i)
                    nmrtbl = nmrtbl + mrt_nlevels
                 ENDDO
@@ -7240,11 +7236,11 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
+                DO  m = surf_usm_h(0)%start_index(j,i), surf_usm_h(0)%end_index(j,i)
 !
 !--                Skip roof if requested
-                   IF ( mrt_skip_roof  .AND.  surf_usm_h%isroof_surf(m) )  CYCLE
+                   IF ( mrt_skip_roof  .AND.  surf_usm_h(0)%isroof_surf(m) )  CYCLE
 !
 !--                Cycle over specified no of levels
-                   l = surf_usm_h%k(m)
+                   l = surf_usm_h(0)%k(m)
                    DO  k = l, l + mrt_nlevels - 1
                       imrt = imrt + 1
@@ -7254,6 +7250,6 @@
 !
 !--             Dtto for LSM
-                DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
-                   l = surf_lsm_h%k(m)
+                DO  m = surf_lsm_h(0)%start_index(j,i), surf_lsm_h(0)%end_index(j,i)
+                   l = surf_lsm_h(0)%k(m)
                    DO  k = l, l + mrt_nlevels - 1
                       imrt = imrt + 1
@@ -10298,21 +10294,22 @@
        SELECT CASE ( TRIM( var ) )
 !--       block of large scale (e.g. RRTMG) radiation output variables
+!--       only upward faced horizontal surfaces are considered here
           CASE ( 'rad_net*' )
              IF ( ALLOCATED( rad_net_av ) ) THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
-                      match_lsm = surf_lsm_h%start_index(j,i) <=               &
-                                  surf_lsm_h%end_index(j,i)
-                      match_usm = surf_usm_h%start_index(j,i) <=               &
-                                  surf_usm_h%end_index(j,i)
+                      match_lsm = surf_lsm_h(0)%start_index(j,i) <=               &
+                                  surf_lsm_h(0)%end_index(j,i)
+                      match_usm = surf_usm_h(0)%start_index(j,i) <=               &
+                                  surf_usm_h(0)%end_index(j,i)
 
                      IF ( match_lsm  .AND.  .NOT. match_usm )  THEN
-                         m = surf_lsm_h%end_index(j,i)
+                         m = surf_lsm_h(0)%end_index(j,i)
                          rad_net_av(j,i) = rad_net_av(j,i) +                   &
-                                         surf_lsm_h%rad_net(m)
+                                         surf_lsm_h(0)%rad_net(m)
                       ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h%end_index(j,i)
+                         m = surf_usm_h(0)%end_index(j,i)
                          rad_net_av(j,i) = rad_net_av(j,i) +                   &
-                                         surf_usm_h%rad_net(m)
+                                         surf_usm_h(0)%rad_net(m)
                       ENDIF
                    ENDDO
@@ -10324,17 +10321,17 @@
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
-                      match_lsm = surf_lsm_h%start_index(j,i) <=               &
-                                  surf_lsm_h%end_index(j,i)
-                      match_usm = surf_usm_h%start_index(j,i) <=               &
-                                  surf_usm_h%end_index(j,i)
+                      match_lsm = surf_lsm_h(0)%start_index(j,i) <=               &
+                                  surf_lsm_h(0)%end_index(j,i)
+                      match_usm = surf_usm_h(0)%start_index(j,i) <=               &
+                                  surf_usm_h(0)%end_index(j,i)
 
                       IF ( match_lsm  .AND.  .NOT. match_usm )  THEN
-                         m = surf_lsm_h%end_index(j,i)
+                         m = surf_lsm_h(0)%end_index(j,i)
                          rad_lw_in_xy_av(j,i) = rad_lw_in_xy_av(j,i) +         &
-                                         surf_lsm_h%rad_lw_in(m)
+                                         surf_lsm_h(0)%rad_lw_in(m)
                       ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h%end_index(j,i)
+                         m = surf_usm_h(0)%end_index(j,i)
                          rad_lw_in_xy_av(j,i) = rad_lw_in_xy_av(j,i) +         &
-                                         surf_usm_h%rad_lw_in(m)
+                                         surf_usm_h(0)%rad_lw_in(m)
                       ENDIF
                    ENDDO
@@ -10346,17 +10343,17 @@
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
-                      match_lsm = surf_lsm_h%start_index(j,i) <=               &
-                                  surf_lsm_h%end_index(j,i)
-                      match_usm = surf_usm_h%start_index(j,i) <=               &
-                                  surf_usm_h%end_index(j,i)
+                      match_lsm = surf_lsm_h(0)%start_index(j,i) <=               &
+                                  surf_lsm_h(0)%end_index(j,i)
+                      match_usm = surf_usm_h(0)%start_index(j,i) <=               &
+                                  surf_usm_h(0)%end_index(j,i)
 
                       IF ( match_lsm  .AND.  .NOT. match_usm )  THEN
-                         m = surf_lsm_h%end_index(j,i)
+                         m = surf_lsm_h(0)%end_index(j,i)
                          rad_lw_out_xy_av(j,i) = rad_lw_out_xy_av(j,i) +       &
-                                                 surf_lsm_h%rad_lw_out(m)
+                                                 surf_lsm_h(0)%rad_lw_out(m)
                       ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h%end_index(j,i)
+                         m = surf_usm_h(0)%end_index(j,i)
                          rad_lw_out_xy_av(j,i) = rad_lw_out_xy_av(j,i) +       &
-                                                 surf_usm_h%rad_lw_out(m)
+                                                 surf_usm_h(0)%rad_lw_out(m)
                       ENDIF
                    ENDDO
@@ -10368,17 +10365,17 @@
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
-                      match_lsm = surf_lsm_h%start_index(j,i) <=               &
-                                  surf_lsm_h%end_index(j,i)
-                      match_usm = surf_usm_h%start_index(j,i) <=               &
-                                  surf_usm_h%end_index(j,i)
+                      match_lsm = surf_lsm_h(0)%start_index(j,i) <=               &
+                                  surf_lsm_h(0)%end_index(j,i)
+                      match_usm = surf_usm_h(0)%start_index(j,i) <=               &
+                                  surf_usm_h(0)%end_index(j,i)
 
                       IF ( match_lsm  .AND.  .NOT. match_usm )  THEN
-                         m = surf_lsm_h%end_index(j,i)
+                         m = surf_lsm_h(0)%end_index(j,i)
                          rad_sw_in_xy_av(j,i) = rad_sw_in_xy_av(j,i) +         &
-                                                surf_lsm_h%rad_sw_in(m)
+                                                surf_lsm_h(0)%rad_sw_in(m)
                       ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h%end_index(j,i)
+                         m = surf_usm_h(0)%end_index(j,i)
                          rad_sw_in_xy_av(j,i) = rad_sw_in_xy_av(j,i) +         &
-                                                surf_usm_h%rad_sw_in(m)
+                                                surf_usm_h(0)%rad_sw_in(m)
                       ENDIF
                    ENDDO
@@ -10390,17 +10387,17 @@
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
-                      match_lsm = surf_lsm_h%start_index(j,i) <=               &
-                                  surf_lsm_h%end_index(j,i)
-                      match_usm = surf_usm_h%start_index(j,i) <=               &
-                                  surf_usm_h%end_index(j,i)
+                      match_lsm = surf_lsm_h(0)%start_index(j,i) <=               &
+                                  surf_lsm_h(0)%end_index(j,i)
+                      match_usm = surf_usm_h(0)%start_index(j,i) <=               &
+                                  surf_usm_h(0)%end_index(j,i)
 
                       IF ( match_lsm  .AND.  .NOT. match_usm )  THEN
-                         m = surf_lsm_h%end_index(j,i)
+                         m = surf_lsm_h(0)%end_index(j,i)
                          rad_sw_out_xy_av(j,i) = rad_sw_out_xy_av(j,i) +       &
-                                                 surf_lsm_h%rad_sw_out(m)
+                                                 surf_lsm_h(0)%rad_sw_out(m)
                       ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h%end_index(j,i)
+                         m = surf_usm_h(0)%end_index(j,i)
                          rad_sw_out_xy_av(j,i) = rad_sw_out_xy_av(j,i) +       &
-                                                 surf_usm_h%rad_sw_out(m)
+                                                 surf_usm_h(0)%rad_sw_out(m)
                       ENDIF
                    ENDDO
@@ -11072,14 +11069,15 @@
 !
 !--                Obtain rad_net from its respective surface type
+!--                Only upward faced horizontal outputs are considered here
 !--                Natural-type surfaces
-                   DO  m = surf_lsm_h%start_index(j,i),                        &
-                           surf_lsm_h%end_index(j,i)
-                      local_pf(i,j,nzb+1) = surf_lsm_h%rad_net(m)
+                   DO  m = surf_lsm_h(0)%start_index(j,i),                        &
+                           surf_lsm_h(0)%end_index(j,i)
+                      local_pf(i,j,nzb+1) = surf_lsm_h(0)%rad_net(m)
                    ENDDO
 !
 !--                Urban-type surfaces
-                   DO  m = surf_usm_h%start_index(j,i),                        &
-                           surf_usm_h%end_index(j,i)
-                      local_pf(i,j,nzb+1) = surf_usm_h%rad_net(m)
+                   DO  m = surf_usm_h(0)%start_index(j,i),                        &
+                           surf_usm_h(0)%end_index(j,i)
+                      local_pf(i,j,nzb+1) = surf_usm_h(0)%rad_net(m)
                    ENDDO
                 ENDDO
@@ -11106,13 +11104,13 @@
 !--                Obtain rad_net from its respective surface type
 !--                Natural-type surfaces
-                   DO  m = surf_lsm_h%start_index(j,i),                        &
-                           surf_lsm_h%end_index(j,i)
-                      local_pf(i,j,nzb+1) = surf_lsm_h%rad_lw_in(m)
+                   DO  m = surf_lsm_h(0)%start_index(j,i),                        &
+                           surf_lsm_h(0)%end_index(j,i)
+                      local_pf(i,j,nzb+1) = surf_lsm_h(0)%rad_lw_in(m)
                    ENDDO
 !
 !--                Urban-type surfaces
-                   DO  m = surf_usm_h%start_index(j,i),                        &
-                           surf_usm_h%end_index(j,i)
-                      local_pf(i,j,nzb+1) = surf_usm_h%rad_lw_in(m)
+                   DO  m = surf_usm_h(0)%start_index(j,i),                        &
+                           surf_usm_h(0)%end_index(j,i)
+                      local_pf(i,j,nzb+1) = surf_usm_h(0)%rad_lw_in(m)
                    ENDDO
                 ENDDO
@@ -11139,13 +11137,13 @@
 !--                Obtain rad_net from its respective surface type
 !--                Natural-type surfaces
-                   DO  m = surf_lsm_h%start_index(j,i),                        &
-                           surf_lsm_h%end_index(j,i)
-                      local_pf(i,j,nzb+1) = surf_lsm_h%rad_lw_out(m)
+                   DO  m = surf_lsm_h(0)%start_index(j,i),                        &
+                           surf_lsm_h(0)%end_index(j,i)
+                      local_pf(i,j,nzb+1) = surf_lsm_h(0)%rad_lw_out(m)
                    ENDDO
 !
 !--                Urban-type surfaces
-                   DO  m = surf_usm_h%start_index(j,i),                        &
-                           surf_usm_h%end_index(j,i)
-                      local_pf(i,j,nzb+1) = surf_usm_h%rad_lw_out(m)
+                   DO  m = surf_usm_h(0)%start_index(j,i),                        &
+                           surf_usm_h(0)%end_index(j,i)
+                      local_pf(i,j,nzb+1) = surf_usm_h(0)%rad_lw_out(m)
                    ENDDO
                 ENDDO
@@ -11172,13 +11170,13 @@
 !--                Obtain rad_net from its respective surface type
 !--                Natural-type surfaces
-                   DO  m = surf_lsm_h%start_index(j,i),                        &
-                           surf_lsm_h%end_index(j,i)
-                      local_pf(i,j,nzb+1) = surf_lsm_h%rad_sw_in(m)
+                   DO  m = surf_lsm_h(0)%start_index(j,i),                        &
+                           surf_lsm_h(0)%end_index(j,i)
+                      local_pf(i,j,nzb+1) = surf_lsm_h(0)%rad_sw_in(m)
                    ENDDO
 !
 !--                Urban-type surfaces
-                   DO  m = surf_usm_h%start_index(j,i),                        &
-                           surf_usm_h%end_index(j,i)
-                      local_pf(i,j,nzb+1) = surf_usm_h%rad_sw_in(m)
+                   DO  m = surf_usm_h(0)%start_index(j,i),                        &
+                           surf_usm_h(0)%end_index(j,i)
+                      local_pf(i,j,nzb+1) = surf_usm_h(0)%rad_sw_in(m)
                    ENDDO
                 ENDDO
@@ -11205,13 +11203,13 @@
 !--                Obtain rad_net from its respective surface type
 !--                Natural-type surfaces
-                   DO  m = surf_lsm_h%start_index(j,i),                        &
-                           surf_lsm_h%end_index(j,i)
-                      local_pf(i,j,nzb+1) = surf_lsm_h%rad_sw_out(m)
+                   DO  m = surf_lsm_h(0)%start_index(j,i),                        &
+                           surf_lsm_h(0)%end_index(j,i)
+                      local_pf(i,j,nzb+1) = surf_lsm_h(0)%rad_sw_out(m)
                    ENDDO
 !
 !--                Urban-type surfaces
-                   DO  m = surf_usm_h%start_index(j,i),                        &
-                           surf_usm_h%end_index(j,i)
-                      local_pf(i,j,nzb+1) = surf_usm_h%rad_sw_out(m)
+                   DO  m = surf_usm_h(0)%start_index(j,i),                        &
+                           surf_usm_h(0)%end_index(j,i)
+                      local_pf(i,j,nzb+1) = surf_usm_h(0)%rad_sw_out(m)
                    ENDDO
                 ENDDO
