Index: /palm/trunk/SOURCE/init_grid.f90
===================================================================
--- /palm/trunk/SOURCE/init_grid.f90	(revision 4652)
+++ /palm/trunk/SOURCE/init_grid.f90	(revision 4653)
@@ -2411,5 +2411,5 @@
 !-- Pre-calculate topography top indices (former get_topography_top_index
 !-- function)
-    ALLOCATE( topo_top_ind(nysg:nyng,nxlg:nxrg,0:4) )
+    ALLOCATE( topo_top_ind(nysg:nyng,nxlg:nxrg,0:5) )
 !
 !-- Uppermost topography index on scalar grid
@@ -2437,4 +2437,10 @@
     topo_top_ind(:,:,4) = MAXLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,:,:), ibit ) ), DIM=1 ) &
                           - 1
+!
+!-- Uppermost topography index including full-3D geometry
+    ibit = 12
+    DO  k = nzb, nzt+1
+       WHERE( BTEST( wall_flags_total_0(k,:,:), ibit ) )  topo_top_ind(:,:,5) = k
+    ENDDO
 
  END SUBROUTINE set_topo_flags
Index: /palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4652)
+++ /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4653)
@@ -28,4 +28,13 @@
 ! -----------------
 ! $Id$
+! RTM version 4.0:
+! - Support full-3D geometry:
+!   * Downward faces can only be default type (temporary PALM limitation)
+!   * rad_angular_discretization = .FALSE. is not supported for full-3D cases
+! - Order faces in surf(l) by column
+! - Remove separate direction indices for urban and land faces
+! - Bugfix: skyvft output
+!
+! 4644 2020-08-22 15:05:46Z pavelkrc
 ! Bugfix: remove invalid reordering of vertical surfaces in RTM
 !
@@ -426,5 +435,5 @@
 
     USE surface_mod,                                                           &
-        ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win,                       &
+        ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_def_h,           &
                surf_lsm_h, surf_lsm_v, surf_type, surf_usm_h, surf_usm_v,      &
                vertical_surfaces_exist
@@ -711,44 +720,26 @@
     INTEGER(iwp), PARAMETER                        ::  iy = 3                             !< position of j-index in surfl and surf
     INTEGER(iwp), PARAMETER                        ::  ix = 4                             !< position of i-index in surfl and surf
-    INTEGER(iwp), PARAMETER                        ::  im = 5                             !< position of surface m-index in surfl and surf
-    INTEGER(iwp), PARAMETER                        ::  nidx_surf = 5                      !< number of indices in surfl and surf
-
-    INTEGER(iwp), PARAMETER                        ::  nsurf_type = 10                    !< number of surf types incl. phys.(land+urban) & (atm.,sky,boundary) surfaces - 1
-
-    INTEGER(iwp), PARAMETER                        ::  iup_u    = 0                       !< 0 - index of urban upward surface (ground or roof)
-    INTEGER(iwp), PARAMETER                        ::  idown_u  = 1                       !< 1 - index of urban downward surface (overhanging)
-    INTEGER(iwp), PARAMETER                        ::  inorth_u = 2                       !< 2 - index of urban northward facing wall
-    INTEGER(iwp), PARAMETER                        ::  isouth_u = 3                       !< 3 - index of urban southward facing wall
-    INTEGER(iwp), PARAMETER                        ::  ieast_u  = 4                       !< 4 - index of urban eastward facing wall
-    INTEGER(iwp), PARAMETER                        ::  iwest_u  = 5                       !< 5 - index of urban westward facing wall
-
-    INTEGER(iwp), PARAMETER                        ::  iup_l    = 6                       !< 6 - index of land upward surface (ground or roof)
-    INTEGER(iwp), PARAMETER                        ::  inorth_l = 7                       !< 7 - index of land northward facing wall
-    INTEGER(iwp), PARAMETER                        ::  isouth_l = 8                       !< 8 - index of land southward facing wall
-    INTEGER(iwp), PARAMETER                        ::  ieast_l  = 9                       !< 9 - index of land eastward facing wall
-    INTEGER(iwp), PARAMETER                        ::  iwest_l  = 10                      !< 10- index of land westward facing wall
-
-    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  idir = (/0, 0,0, 0,1,-1,0,0, 0,1,-1/)   !< surface normal direction x indices
-    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  jdir = (/0, 0,1,-1,0, 0,0,1,-1,0, 0/)   !< surface normal direction y indices
-    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  kdir = (/1,-1,0, 0,0, 0,1,0, 0,0, 0/)   !< surface normal direction z indices
+    INTEGER(iwp), PARAMETER                        ::  nidx_surf = 4                      !< number of indices in surfl and surf
+
+    INTEGER(iwp), PARAMETER                        ::  nsurf_type = 5                     !< number of surf types = surface directions
+
+    INTEGER(iwp), PARAMETER                        ::  iup    = 0                         !< 0 - index of upward surface (ground or roof)
+    INTEGER(iwp), PARAMETER                        ::  idown  = 1                         !< 1 - index of downward surface (overhanging)
+    INTEGER(iwp), PARAMETER                        ::  inorth = 2                         !< 2 - index of northward facing wall
+    INTEGER(iwp), PARAMETER                        ::  isouth = 3                         !< 3 - index of southward facing wall
+    INTEGER(iwp), PARAMETER                        ::  ieast  = 4                         !< 4 - index of eastward facing wall
+    INTEGER(iwp), PARAMETER                        ::  iwest  = 5                         !< 5 - index of westward facing wall
+
+    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  idir = (/0, 0, 0, 0, 1,-1/)      !< surface normal direction x indices
+    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  jdir = (/0, 0, 1,-1, 0, 0/)      !< surface normal direction y indices
+    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  kdir = (/1,-1, 0, 0, 0, 0/)      !< surface normal direction z indices
     REAL(wp),     DIMENSION(0:nsurf_type)          :: facearea                            !< area of single face in respective
                                                                                           !< direction (will be calc'd)
 
-
-!-- indices and sizes of urban and land surface models
-    INTEGER(iwp)                                   ::  startland        !< start index of block of land and roof surfaces
-    INTEGER(iwp)                                   ::  endland          !< end index of block of land and roof surfaces
-    INTEGER(iwp)                                   ::  nlands           !< number of land and roof surfaces in local processor
-    INTEGER(iwp)                                   ::  startwall        !< start index of block of wall surfaces
-    INTEGER(iwp)                                   ::  endwall          !< end index of block of wall surfaces
-    INTEGER(iwp)                                   ::  nwalls           !< number of wall surfaces in local processor
 
 !-- indices needed for RTM netcdf output subroutines
     INTEGER(iwp), PARAMETER                        :: nd = 5
     CHARACTER(LEN=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
-    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     :: dirint_u = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)
-    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     :: dirint_l = (/ iup_l, isouth_l, inorth_l, iwest_l, ieast_l /)
-    INTEGER(iwp), DIMENSION(0:nd-1)                :: dirstart
-    INTEGER(iwp), DIMENSION(0:nd-1)                :: dirend
+    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     :: dirint = (/ iup, isouth, inorth, iwest, ieast /)
 
 !-- indices and sizes of urban and land surface models
@@ -762,4 +753,6 @@
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfstart        !< starts of blocks of surfaces for individual processors in array surf (indexed from 1)
                                                                         !< respective block for particular processor is surfstart[iproc+1]+1 : surfstart[iproc+1]+nsurfs[iproc+1]
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfl_col_start  !< start of surfaces in surfl for each x,y column (local surfaces)
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfg_col_start  !< start of surfaces in surfl for each x,y column (all surfaces)
 
 !-- block variables needed for calculation of the plant canopy model inside the urban surface model
@@ -786,5 +779,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. 3.0'         !< identification of version of binary svf and restart files
+    CHARACTER(rad_version_len), PARAMETER          ::  rad_version = 'RAD v. 4.0'         !< 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)
@@ -879,5 +872,7 @@
 #endif
     REAL(wp)                                       ::  prototype_lad    !< prototype leaf area density for computing effective optical depth
-    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  nzterr, plantt   !< temporary global arrays for raytracing
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  nzterrt          !< temporary global arrays for raytracing
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  nzterrb          !< temporary global arrays for raytracing
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  plantt           !< temporary global arrays for raytracing
     INTEGER(iwp)                                   ::  plantt_max
 
@@ -1093,11 +1088,9 @@
            cos_zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon,   &
            idir, jdir, kdir, id, iz, iy, ix,                                   &
-           iup_u, inorth_u, isouth_u, ieast_u, iwest_u,                        &
-           iup_l, inorth_l, isouth_l, ieast_l, iwest_l,                        &
-           nsurf_type, nz_urban_b, nz_urban_t, nz_urban, pch, nsurf,                 &
+           iup, inorth, isouth, ieast, iwest,                                  &
+           nsurf_type, nz_urban_b, nz_urban_t, nz_urban, pch, nsurf,           &
            idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct,                             &
-           radiation_interactions, startwall, startland, endland, endwall,     &
-           skyvf, skyvft, radiation_interactions_on, average_radiation,        &
-           rad_sw_in_diff, rad_sw_in_dir
+           radiation_interactions, skyvf, skyvft, radiation_interactions_on,   &
+           average_radiation, rad_sw_in_diff, rad_sw_in_dir
 
 
@@ -5966,9 +5959,9 @@
 !
 !--  Set up thermal radiation from surfaces
-!--  Horizontal walls
      mm = 1
      DO  i = nxl, nxr
         DO  j = nys, nyn
 !
+!--        Horizontal walls
 !--        urban
            DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
@@ -5996,10 +5989,18 @@
               mm = mm + 1
            ENDDO
-        ENDDO
-     ENDDO
-!
-!--  Vertical walls
-     DO  i = nxl, nxr
-        DO  j = nys, nyn
+!
+!--        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
+           ENDDO
+!
+!--        Vertical walls
            DO  l = 0, 3
 !
@@ -6400,186 +6401,106 @@
 !--  Transfer radiation arrays required for energy balance to the respective data types
 !    and claculate relevant radiation model-RTM coupling terms
+     mm = 1
+     DO  i = nxl, nxr
+        DO  j = nys, nyn
+!
+!--        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) -        &
+                                            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
+           ENDDO
+!
+!--        Vertical walls
+           DO  l = 0, 3
+!
+!--           urban
+              DO  m = surf_usm_v(l)%start_index(j,i),                       &
+                      surf_usm_v(l)%end_index(j,i)
+                 surf_usm_v(l)%rad_sw_in(m)  = surfinsw(mm)
+                 surf_usm_v(l)%rad_sw_out(m) = surfoutsw(mm)
+                 surf_usm_v(l)%rad_sw_dir(m) = surfinswdir(mm)
+                 surf_usm_v(l)%rad_sw_dif(m) = surfinswdif(mm)
+                 surf_usm_v(l)%rad_sw_ref(m) = surfinsw(mm) - surfinswdir(mm) -     &
+                                               surfinswdif(mm)
+                 surf_usm_v(l)%rad_sw_res(m) = surfins(mm)
+                 surf_usm_v(l)%rad_lw_in(m)  = surfinlw(mm)
+                 surf_usm_v(l)%rad_lw_out(m) = surfoutlw(mm)
+                 surf_usm_v(l)%rad_net(m)    = surfinsw(mm) - surfoutsw(mm) +       &
+                                               surfinlw(mm) - surfoutlw(mm)
+                 surf_usm_v(l)%rad_net_l(m)  = surf_usm_v(l)%rad_net(m)
+                 surf_usm_v(l)%rad_lw_dif(m) = surfinlwdif(mm)
+                 surf_usm_v(l)%rad_lw_ref(m) = surfinlw(mm) - surfinlwdif(mm)
+                 surf_usm_v(l)%rad_lw_res(m) = surfinl(mm)
+                 mm = mm + 1
+              ENDDO
+!
+!--           land
+              DO  m = surf_lsm_v(l)%start_index(j,i),                       &
+                      surf_lsm_v(l)%end_index(j,i)
+                 surf_lsm_v(l)%rad_sw_in(m)  = surfinsw(mm)
+                 surf_lsm_v(l)%rad_sw_out(m) = surfoutsw(mm)
+                 surf_lsm_v(l)%rad_sw_dir(m) = surfinswdir(mm)
+                 surf_lsm_v(l)%rad_sw_dif(m) = surfinswdif(mm)
+                 surf_lsm_v(l)%rad_sw_ref(m) = surfinsw(mm) - surfinswdir(mm) -     &
+                                               surfinswdif(mm)
+                 surf_lsm_v(l)%rad_sw_res(m) = surfins(mm)
+                 surf_lsm_v(l)%rad_lw_in(m)  = surfinlw(mm)
+                 surf_lsm_v(l)%rad_lw_out(m) = surfoutlw(mm)
+                 surf_lsm_v(l)%rad_net(m)    = surfinsw(mm) - surfoutsw(mm) +       &
+                                               surfinlw(mm) - surfoutlw(mm)
+                 surf_lsm_v(l)%rad_lw_dif(m) = surfinlwdif(mm)
+                 surf_lsm_v(l)%rad_lw_ref(m) = surfinlw(mm) - surfinlwdif(mm)
+                 surf_lsm_v(l)%rad_lw_res(m) = surfinl(mm)
+                 mm = mm + 1
+              ENDDO
+           ENDDO
+        ENDDO
+     ENDDO
 
      DO  i = 1, nsurfl
-        m  = surfl(im,i)
         d  = surfl(id, i)
-!
-!--     (1) Urban surfaces
-!--     upward-facing
-        IF ( surfl(1,i) == iup_u )  THEN
-           surf_usm_h%rad_sw_in(m)  = surfinsw(i)
-           surf_usm_h%rad_sw_out(m) = surfoutsw(i)
-           surf_usm_h%rad_sw_dir(m) = surfinswdir(i)
-           surf_usm_h%rad_sw_dif(m) = surfinswdif(i)
-           surf_usm_h%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) -        &
-                                      surfinswdif(i)
-           surf_usm_h%rad_sw_res(m) = surfins(i)
-           surf_usm_h%rad_lw_in(m)  = surfinlw(i)
-           surf_usm_h%rad_lw_out(m) = surfoutlw(i)
-           surf_usm_h%rad_net(m)    = surfinsw(i) - surfoutsw(i) +          &
-                                      surfinlw(i) - surfoutlw(i)
-           surf_usm_h%rad_net_l(m)  = surf_usm_h%rad_net(m)
-           surf_usm_h%rad_lw_dif(m) = surfinlwdif(i)
-           surf_usm_h%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)
-           surf_usm_h%rad_lw_res(m) = surfinl(i)
-!
-!--     northward-facding
-        ELSEIF ( surfl(1,i) == inorth_u )  THEN
-           surf_usm_v(0)%rad_sw_in(m)  = surfinsw(i)
-           surf_usm_v(0)%rad_sw_out(m) = surfoutsw(i)
-           surf_usm_v(0)%rad_sw_dir(m) = surfinswdir(i)
-           surf_usm_v(0)%rad_sw_dif(m) = surfinswdif(i)
-           surf_usm_v(0)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) -     &
-                                         surfinswdif(i)
-           surf_usm_v(0)%rad_sw_res(m) = surfins(i)
-           surf_usm_v(0)%rad_lw_in(m)  = surfinlw(i)
-           surf_usm_v(0)%rad_lw_out(m) = surfoutlw(i)
-           surf_usm_v(0)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
-                                         surfinlw(i) - surfoutlw(i)
-           surf_usm_v(0)%rad_net_l(m)  = surf_usm_v(0)%rad_net(m)
-           surf_usm_v(0)%rad_lw_dif(m) = surfinlwdif(i)
-           surf_usm_v(0)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)
-           surf_usm_v(0)%rad_lw_res(m) = surfinl(i)
-!
-!--     southward-facding
-        ELSEIF ( surfl(1,i) == isouth_u )  THEN
-           surf_usm_v(1)%rad_sw_in(m)  = surfinsw(i)
-           surf_usm_v(1)%rad_sw_out(m) = surfoutsw(i)
-           surf_usm_v(1)%rad_sw_dir(m) = surfinswdir(i)
-           surf_usm_v(1)%rad_sw_dif(m) = surfinswdif(i)
-           surf_usm_v(1)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) -     &
-                                         surfinswdif(i)
-           surf_usm_v(1)%rad_sw_res(m) = surfins(i)
-           surf_usm_v(1)%rad_lw_in(m)  = surfinlw(i)
-           surf_usm_v(1)%rad_lw_out(m) = surfoutlw(i)
-           surf_usm_v(1)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
-                                         surfinlw(i) - surfoutlw(i)
-           surf_usm_v(1)%rad_net_l(m)  = surf_usm_v(1)%rad_net(m)
-           surf_usm_v(1)%rad_lw_dif(m) = surfinlwdif(i)
-           surf_usm_v(1)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)
-           surf_usm_v(1)%rad_lw_res(m) = surfinl(i)
-!
-!--     eastward-facing
-        ELSEIF ( surfl(1,i) == ieast_u )  THEN
-           surf_usm_v(2)%rad_sw_in(m)  = surfinsw(i)
-           surf_usm_v(2)%rad_sw_out(m) = surfoutsw(i)
-           surf_usm_v(2)%rad_sw_dir(m) = surfinswdir(i)
-           surf_usm_v(2)%rad_sw_dif(m) = surfinswdif(i)
-           surf_usm_v(2)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) -     &
-                                         surfinswdif(i)
-           surf_usm_v(2)%rad_sw_res(m) = surfins(i)
-           surf_usm_v(2)%rad_lw_in(m)  = surfinlw(i)
-           surf_usm_v(2)%rad_lw_out(m) = surfoutlw(i)
-           surf_usm_v(2)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
-                                         surfinlw(i) - surfoutlw(i)
-           surf_usm_v(2)%rad_net_l(m)  = surf_usm_v(2)%rad_net(m)
-           surf_usm_v(2)%rad_lw_dif(m) = surfinlwdif(i)
-           surf_usm_v(2)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)
-           surf_usm_v(2)%rad_lw_res(m) = surfinl(i)
-!
-!--     westward-facding
-        ELSEIF ( surfl(1,i) == iwest_u )  THEN
-           surf_usm_v(3)%rad_sw_in(m)  = surfinsw(i)
-           surf_usm_v(3)%rad_sw_out(m) = surfoutsw(i)
-           surf_usm_v(3)%rad_sw_dir(m) = surfinswdir(i)
-           surf_usm_v(3)%rad_sw_dif(m) = surfinswdif(i)
-           surf_usm_v(3)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) -     &
-                                         surfinswdif(i)
-           surf_usm_v(3)%rad_sw_res(m) = surfins(i)
-           surf_usm_v(3)%rad_lw_in(m)  = surfinlw(i)
-           surf_usm_v(3)%rad_lw_out(m) = surfoutlw(i)
-           surf_usm_v(3)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
-                                         surfinlw(i) - surfoutlw(i)
-           surf_usm_v(3)%rad_net_l(m)  = surf_usm_v(3)%rad_net(m)
-           surf_usm_v(3)%rad_lw_dif(m) = surfinlwdif(i)
-           surf_usm_v(3)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)
-           surf_usm_v(3)%rad_lw_res(m) = surfinl(i)
-!
-!--     (2) land surfaces
-!--     upward-facing
-        ELSEIF ( surfl(1,i) == iup_l )  THEN
-           surf_lsm_h%rad_sw_in(m)  = surfinsw(i)
-           surf_lsm_h%rad_sw_out(m) = surfoutsw(i)
-           surf_lsm_h%rad_sw_dir(m) = surfinswdir(i)
-           surf_lsm_h%rad_sw_dif(m) = surfinswdif(i)
-           surf_lsm_h%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) -        &
-                                         surfinswdif(i)
-           surf_lsm_h%rad_sw_res(m) = surfins(i)
-           surf_lsm_h%rad_lw_in(m)  = surfinlw(i)
-           surf_lsm_h%rad_lw_out(m) = surfoutlw(i)
-           surf_lsm_h%rad_net(m)    = surfinsw(i) - surfoutsw(i) +          &
-                                      surfinlw(i) - surfoutlw(i)
-           surf_lsm_h%rad_lw_dif(m) = surfinlwdif(i)
-           surf_lsm_h%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)
-           surf_lsm_h%rad_lw_res(m) = surfinl(i)
-!
-!--     northward-facding
-        ELSEIF ( surfl(1,i) == inorth_l )  THEN
-           surf_lsm_v(0)%rad_sw_in(m)  = surfinsw(i)
-           surf_lsm_v(0)%rad_sw_out(m) = surfoutsw(i)
-           surf_lsm_v(0)%rad_sw_dir(m) = surfinswdir(i)
-           surf_lsm_v(0)%rad_sw_dif(m) = surfinswdif(i)
-           surf_lsm_v(0)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) -     &
-                                         surfinswdif(i)
-           surf_lsm_v(0)%rad_sw_res(m) = surfins(i)
-           surf_lsm_v(0)%rad_lw_in(m)  = surfinlw(i)
-           surf_lsm_v(0)%rad_lw_out(m) = surfoutlw(i)
-           surf_lsm_v(0)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
-                                         surfinlw(i) - surfoutlw(i)
-           surf_lsm_v(0)%rad_lw_dif(m) = surfinlwdif(i)
-           surf_lsm_v(0)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)
-           surf_lsm_v(0)%rad_lw_res(m) = surfinl(i)
-!
-!--     southward-facding
-        ELSEIF ( surfl(1,i) == isouth_l )  THEN
-           surf_lsm_v(1)%rad_sw_in(m)  = surfinsw(i)
-           surf_lsm_v(1)%rad_sw_out(m) = surfoutsw(i)
-           surf_lsm_v(1)%rad_sw_dir(m) = surfinswdir(i)
-           surf_lsm_v(1)%rad_sw_dif(m) = surfinswdif(i)
-           surf_lsm_v(1)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) -     &
-                                         surfinswdif(i)
-           surf_lsm_v(1)%rad_sw_res(m) = surfins(i)
-           surf_lsm_v(1)%rad_lw_in(m)  = surfinlw(i)
-           surf_lsm_v(1)%rad_lw_out(m) = surfoutlw(i)
-           surf_lsm_v(1)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
-                                         surfinlw(i) - surfoutlw(i)
-           surf_lsm_v(1)%rad_lw_dif(m) = surfinlwdif(i)
-           surf_lsm_v(1)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)
-           surf_lsm_v(1)%rad_lw_res(m) = surfinl(i)
-!
-!--     eastward-facing
-        ELSEIF ( surfl(1,i) == ieast_l )  THEN
-           surf_lsm_v(2)%rad_sw_in(m)  = surfinsw(i)
-           surf_lsm_v(2)%rad_sw_out(m) = surfoutsw(i)
-           surf_lsm_v(2)%rad_sw_dir(m) = surfinswdir(i)
-           surf_lsm_v(2)%rad_sw_dif(m) = surfinswdif(i)
-           surf_lsm_v(2)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) -     &
-                                         surfinswdif(i)
-           surf_lsm_v(2)%rad_sw_res(m) = surfins(i)
-           surf_lsm_v(2)%rad_lw_in(m)  = surfinlw(i)
-           surf_lsm_v(2)%rad_lw_out(m) = surfoutlw(i)
-           surf_lsm_v(2)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
-                                         surfinlw(i) - surfoutlw(i)
-           surf_lsm_v(2)%rad_lw_dif(m) = surfinlwdif(i)
-           surf_lsm_v(2)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)
-           surf_lsm_v(2)%rad_lw_res(m) = surfinl(i)
-!
-!--     westward-facing
-        ELSEIF ( surfl(1,i) == iwest_l )  THEN
-           surf_lsm_v(3)%rad_sw_in(m)  = surfinsw(i)
-           surf_lsm_v(3)%rad_sw_out(m) = surfoutsw(i)
-           surf_lsm_v(3)%rad_sw_dir(m) = surfinswdir(i)
-           surf_lsm_v(3)%rad_sw_dif(m) = surfinswdif(i)
-           surf_lsm_v(3)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) -     &
-                                         surfinswdif(i)
-           surf_lsm_v(3)%rad_sw_res(m) = surfins(i)
-           surf_lsm_v(3)%rad_lw_in(m)  = surfinlw(i)
-           surf_lsm_v(3)%rad_lw_out(m) = surfoutlw(i)
-           surf_lsm_v(3)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
-                                         surfinlw(i) - surfoutlw(i)
-           surf_lsm_v(3)%rad_lw_dif(m) = surfinlwdif(i)
-           surf_lsm_v(3)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)
-           surf_lsm_v(3)%rad_lw_res(m) = surfinl(i)
-        ENDIF
 !
 !--     RTM coupling terms
@@ -6942,4 +6863,5 @@
        INTEGER(iwp)        ::  k_topo               !< vertical index indicating
                                                     !< topography top for given (j,i)
+       INTEGER(iwp)        ::  icol                 !< flat column number (in (y,x) fortran order)
        INTEGER(iwp)        ::  isurf, ipcgb, imrt
        INTEGER(iwp)        ::  nzptl, nzubl, nzutl
@@ -6965,5 +6887,5 @@
 !--    for any upward-facing wall (see bit 12).
        nzubl = MINVAL( topo_top_ind(nys:nyn,nxl:nxr,0) )
-       nzutl = MAXVAL( topo_top_ind(nys:nyn,nxl:nxr,0) )
+       nzutl = MAXVAL( topo_top_ind(nys:nyn,nxl:nxr,5) )
 
        nzubl = MAX( nzubl, nzb )
@@ -7077,20 +6999,11 @@
 !--    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.
-       startland = 1
-       nsurfl    = surf_usm_h%ns + surf_lsm_h%ns
-       endland   = nsurfl
-       nlands    = endland - startland + 1
-
+       nsurfl    = surf_usm_h%ns + surf_lsm_h%ns + surf_def_h(1)%ns
 !
 !--    Number of vertical surfaces in both USM and LSM. Note that all vertical surface elements are
 !--    already counted in surface_mod.
-       startwall = nsurfl+1
        DO  i = 0,3
           nsurfl = nsurfl + surf_usm_v(i)%ns + surf_lsm_v(i)%ns
        ENDDO
-       endwall = nsurfl
-       nwalls  = endwall - startwall + 1
-       dirstart = (/ startland, startwall, startwall, startwall, startwall /)
-       dirend = (/ endland, endwall, endwall, endwall, endwall /)
 
 !--    fill gridpcbl and pcbl
@@ -7196,15 +7109,22 @@
        ENDIF
 
-!--    add horizontal surface elements (land and urban surfaces)
-!--    TODO: add urban overhanging surfaces (idown_u)
+       ALLOCATE( surfl_col_start(0:nnx*nny-1), surfg_col_start(0:(nx+1)*(ny+1)) )
+!
+!--    Add horizontal and vertical surface elements (land and urban surfaces)
+!--    ordered by x,y column (y most varying)
+!--    TODO: remove the hard coding of l = 0 to l = idirection
+       icol = 0
        DO i = nxl, nxr
            DO j = nys, nyn
+!
+!--           Save column start
+              surfl_col_start(icol) = isurf + 1
+              icol = icol + 1
+!
+!--           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_u,k,j,i,m/)
-                 IF ( rad_angular_discretization ) THEN
-                    gridsurf(iup_u,k,j,i) = isurf
-                 ENDIF
+                 surfl(:,isurf) = (/iup,k,j,i/)
               ENDDO
 
@@ -7212,33 +7132,24 @@
                  k = surf_lsm_h%k(m)
                  isurf = isurf + 1
-                 surfl(:,isurf) = (/iup_l,k,j,i,m/)
-                 IF ( rad_angular_discretization ) THEN
-                    gridsurf(iup_u,k,j,i) = isurf
-                 ENDIF
+                 surfl(:,isurf) = (/iup,k,j,i/)
               ENDDO
 
-           ENDDO
-       ENDDO
-
-!--    add vertical surface elements (land and urban surfaces)
-!--    TODO: remove the hard coding of l = 0 to l = idirection
-       DO i = nxl, nxr
-           DO j = nys, nyn
+              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_u,k,j,i,m/)
-                 IF ( rad_angular_discretization ) THEN
-                    gridsurf(inorth_u,k,j,i) = isurf
-                 ENDIF
+                 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_l,k,j,i,m/)
-                 IF ( rad_angular_discretization ) THEN
-                    gridsurf(inorth_u,k,j,i) = isurf
-                 ENDIF
+                 surfl(:,isurf) = (/inorth,k,j,i/)
               ENDDO
 
@@ -7247,16 +7158,10 @@
                  k = surf_usm_v(l)%k(m)
                  isurf = isurf + 1
-                 surfl(:,isurf) = (/isouth_u,k,j,i,m/)
-                 IF ( rad_angular_discretization ) THEN
-                    gridsurf(isouth_u,k,j,i) = isurf
-                 ENDIF
+                 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_l,k,j,i,m/)
-                 IF ( rad_angular_discretization ) THEN
-                    gridsurf(isouth_u,k,j,i) = isurf
-                 ENDIF
+                 surfl(:,isurf) = (/isouth,k,j,i/)
               ENDDO
 
@@ -7265,16 +7170,10 @@
                  k = surf_usm_v(l)%k(m)
                  isurf = isurf + 1
-                 surfl(:,isurf) = (/ieast_u,k,j,i,m/)
-                 IF ( rad_angular_discretization ) THEN
-                    gridsurf(ieast_u,k,j,i) = isurf
-                 ENDIF
+                 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_l,k,j,i,m/)
-                 IF ( rad_angular_discretization ) THEN
-                    gridsurf(ieast_u,k,j,i) = isurf
-                 ENDIF
+                 surfl(:,isurf) = (/ieast,k,j,i/)
               ENDDO
 
@@ -7283,16 +7182,10 @@
                  k = surf_usm_v(l)%k(m)
                  isurf = isurf + 1
-                 surfl(:,isurf) = (/iwest_u,k,j,i,m/)
-                 IF ( rad_angular_discretization ) THEN
-                    gridsurf(iwest_u,k,j,i) = isurf
-                 ENDIF
+                 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_l,k,j,i,m/)
-                 IF ( rad_angular_discretization ) THEN
-                    gridsurf(iwest_u,k,j,i) = isurf
-                 ENDIF
+                 surfl(:,isurf) = (/iwest,k,j,i/)
               ENDDO
            ENDDO
@@ -7406,6 +7299,32 @@
        surf = surfl
 #endif
-
-!--
+!
+!--    Populate gridsurf with reverse global indices (->surf)
+       IF ( rad_angular_discretization )  THEN
+          DO  isurf = 1, nsurfl
+             gridsurf(surfl(id,isurf),surfl(iz,isurf), &
+                      surfl(iy,isurf),surfl(ix,isurf)) = isurf + surfstart(myid)
+          ENDDO
+       ENDIF
+!
+!--    Gather column start indices
+#if defined( __parallel )
+       CALL MPI_Allgather(surfl_col_start, nnx*nny, MPI_INTEGER,             &
+                          surfg_col_start, nnx*nny, MPI_INTEGER, comm2d, ierr)
+       IF ( ierr /= 0 ) THEN
+          WRITE(9,*) 'Error MPI_AllGather1b:', ierr, surfl_col_start, surfg_col_start
+          FLUSH(9)
+       ENDIF
+!
+!--    Convert local indices (->surfl) to global (->surf)
+       DO  i = 0, numprocs-1
+          surfg_col_start(i*nnx*nny:(i+1)*nnx*nny-1) =                      &
+                surfg_col_start(i*nnx*nny:(i+1)*nnx*nny-1) + surfstart(i)
+       ENDDO
+#else
+       surfg_col_start(0:(nx+1)*(ny+1)-1) = surfl_col_start(0:(nx+1)*(ny+1)-1)
+#endif
+       surfg_col_start((nx+1)*(ny+1)) = nsurf+1
+!
 !--    allocation of the arrays for direct and diffusion radiation
        IF ( debug_output )  CALL debug_message( 'allocation of radiation arrays', 'info' )
@@ -7483,5 +7402,4 @@
         REAL(wp), DIMENSION(:), ALLOCATABLE           :: vffrac0       !< dtto (original values)
         REAL(wp), DIMENSION(:), ALLOCATABLE           :: ztransp       !< array of transparency in z steps
-        INTEGER(iwp)                                  :: lowest_free_ray !< index into zdirs
         INTEGER(iwp), DIMENSION(:), ALLOCATABLE       :: itarget       !< face indices of detected obstacles
         INTEGER(iwp)                                  :: itarg0, itarg1
@@ -7502,6 +7420,8 @@
         INTEGER(iwp)                                  :: max_track_len !< maximum 2d track length
 #if defined( __parallel )
-        INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET:: nzterrl_l
-        INTEGER(iwp), DIMENSION(:,:), POINTER         :: nzterrl
+        INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET:: nzterrtl_l
+        INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET:: nzterrbl_l
+        INTEGER(iwp), DIMENSION(:,:), POINTER         :: nzterrtl
+        INTEGER(iwp), DIMENSION(:,:), POINTER         :: nzterrbl
         INTEGER(iwp)                                  :: minfo
         REAL(wp), DIMENSION(:), POINTER, SAVE         :: lad_s_rma       !< fortran 1D pointer
@@ -7540,20 +7460,32 @@
 
 !--     initialize temporary terrain and plant canopy height arrays (global 2D array!)
-        ALLOCATE( nzterr(0:(nx+1)*(ny+1)-1) )
+        ALLOCATE( nzterrt(0:(nx+1)*(ny+1)-1) )
+        ALLOCATE( nzterrb(0:(nx+1)*(ny+1)-1) )
 #if defined( __parallel )
-        !ALLOCATE( nzterrl(nys:nyn,nxl:nxr) )
-        ALLOCATE( nzterrl_l((nyn-nys+1)*(nxr-nxl+1)) )
-        nzterrl(nys:nyn,nxl:nxr) => nzterrl_l(1:(nyn-nys+1)*(nxr-nxl+1))
-        nzterrl = topo_top_ind(nys:nyn,nxl:nxr,0)
-        CALL MPI_AllGather( nzterrl_l, nnx*nny, MPI_INTEGER, &
-                            nzterr, nnx*nny, MPI_INTEGER, comm2d, ierr )
+        ALLOCATE( nzterrtl_l((nyn-nys+1)*(nxr-nxl+1)) )
+        ALLOCATE( nzterrbl_l((nyn-nys+1)*(nxr-nxl+1)) )
+        nzterrtl(nys:nyn,nxl:nxr) => nzterrtl_l(1:(nyn-nys+1)*(nxr-nxl+1))
+        nzterrbl(nys:nyn,nxl:nxr) => nzterrbl_l(1:(nyn-nys+1)*(nxr-nxl+1))
+        nzterrtl = topo_top_ind(nys:nyn,nxl:nxr,5)
+        nzterrbl = topo_top_ind(nys:nyn,nxl:nxr,0)
+        CALL MPI_AllGather( nzterrtl_l, nnx*nny, MPI_INTEGER, &
+                            nzterrt, nnx*nny, MPI_INTEGER, comm2d, ierr )
         IF ( ierr /= 0 ) THEN
-            WRITE(9,*) 'Error MPI_AllGather1:', ierr, SIZE(nzterrl_l), nnx*nny, &
-                       SIZE(nzterr), nnx*nny
+            WRITE(9,*) 'Error MPI_AllGather1t:', ierr, SIZE(nzterrtl_l), nnx*nny, &
+                       SIZE(nzterrt), nnx*nny
             FLUSH(9)
         ENDIF
-        DEALLOCATE(nzterrl_l)
+        CALL MPI_AllGather( nzterrbl_l, nnx*nny, MPI_INTEGER, &
+                            nzterrb, nnx*nny, MPI_INTEGER, comm2d, ierr )
+        IF ( ierr /= 0 ) THEN
+            WRITE(9,*) 'Error MPI_AllGather1b:', ierr, SIZE(nzterrbl_l), nnx*nny, &
+                       SIZE(nzterrb), nnx*nny
+            FLUSH(9)
+        ENDIF
+        DEALLOCATE(nzterrtl_l)
+        DEALLOCATE(nzterrbl_l)
 #else
-        nzterr = RESHAPE( topo_top_ind(nys:nyn,nxl:nxr,0), (/(nx+1)*(ny+1)/) )
+        nzterrt = RESHAPE( topo_top_ind(nys:nyn,nxl:nxr,5), (/(nx+1)*(ny+1)/) )
+        nzterrb = RESHAPE( topo_top_ind(nys:nyn,nxl:nxr,0), (/(nx+1)*(ny+1)/) )
 #endif
         IF ( plant_canopy )  THEN
@@ -7701,5 +7633,5 @@
             !--Select a proper half-sphere for 2D raytracing
             SELECT CASE ( td )
-               CASE ( iup_u, iup_l )
+               CASE ( iup )
                   az0 = 0._wp
                   naz = raytrace_discrete_azims
@@ -7708,5 +7640,12 @@
                   nzn = raytrace_discrete_elevs / 2
                   zns = pi / 2._wp / REAL(nzn, wp)
-               CASE ( isouth_u, isouth_l )
+               CASE ( idown )
+                  az0 = 0._wp
+                  naz = raytrace_discrete_azims
+                  azs = 2._wp * pi / REAL(naz, wp)
+                  zn0 = pi / 2._wp
+                  nzn = raytrace_discrete_elevs / 2
+                  zns = pi / 2._wp / REAL(nzn, wp)
+               CASE ( isouth )
                   az0 = pi / 2._wp
                   naz = raytrace_discrete_azims / 2
@@ -7715,5 +7654,5 @@
                   nzn = raytrace_discrete_elevs
                   zns = pi / REAL(nzn, wp)
-               CASE ( inorth_u, inorth_l )
+               CASE ( inorth )
                   az0 = - pi / 2._wp
                   naz = raytrace_discrete_azims / 2
@@ -7722,5 +7661,5 @@
                   nzn = raytrace_discrete_elevs
                   zns = pi / REAL(nzn, wp)
-               CASE ( iwest_u, iwest_l )
+               CASE ( iwest )
                   az0 = pi
                   naz = raytrace_discrete_azims / 2
@@ -7729,5 +7668,5 @@
                   nzn = raytrace_discrete_elevs
                   zns = pi / REAL(nzn, wp)
-               CASE ( ieast_u, ieast_l )
+               CASE ( ieast )
                   az0 = 0._wp
                   naz = raytrace_discrete_azims / 2
@@ -7744,6 +7683,5 @@
 
             ALLOCATE ( zdirs(1:nzn), zcent(1:nzn), zbdry(0:nzn), vffrac(1:nzn*naz), &
-                       ztransp(1:nzn*naz), itarget(1:nzn*naz) )   !FIXME allocate itarget only
-                                                                  !in case of rad_angular_discretization
+                       ztransp(1:nzn*naz), itarget(1:nzn*naz) )
 
             itarg0 = 1
@@ -7751,5 +7689,5 @@
             zcent(:) = (/( zn0+(REAL(izn,wp)-.5_wp)*zns, izn=1, nzn )/)
             zbdry(:) = (/( zn0+REAL(izn,wp)*zns, izn=0, nzn )/)
-            IF ( td == iup_u  .OR.  td == iup_l )  THEN
+            IF ( td == iup )  THEN
                vffrac(1:nzn) = (COS(2 * zbdry(0:nzn-1)) - COS(2 * zbdry(1:nzn))) / 2._wp / REAL(naz, wp)
 !
@@ -7759,4 +7697,10 @@
                ENDDO
 !--            sum of whole vffrac equals 1, verified
+            ELSEIF ( td == idown )  THEN
+               vffrac(1:nzn) = -(COS(2 * zbdry(0:nzn-1)) - COS(2 * zbdry(1:nzn))) / 2._wp / REAL(naz, wp)
+               DO iaz = 1, naz-1
+                  vffrac(iaz*nzn+1:(iaz+1)*nzn) = vffrac(1:nzn)
+               ENDDO
+!--            sum of whole vffrac equals 1, verified
             ENDIF
 !
@@ -7764,5 +7708,5 @@
             DO iaz = 1, naz
                azmid = az0 + (REAL(iaz, wp) - .5_wp) * azs
-               IF ( td /= iup_u  .AND.  td /= iup_l )  THEN
+               IF ( td /= iup  .AND.  td /= idown )  THEN
                   az2 = REAL(iaz, wp) * azs - pi/2._wp
                   az1 = az2 - azs
@@ -7783,25 +7727,27 @@
                                     surfstart(myid) + isurflt, facearea(td),  &
                                     vffrac(itarg0:itarg1), .TRUE., .TRUE.,    &
-                                    .FALSE., lowest_free_ray,                 &
-                                    ztransp(itarg0:itarg1),                   &
+                                    .FALSE., ztransp(itarg0:itarg1),          &
                                     itarget(itarg0:itarg1))
 
-               skyvf(isurflt) = skyvf(isurflt) + &
-                                SUM(vffrac(itarg0:itarg0+lowest_free_ray-1))
-               skyvft(isurflt) = skyvft(isurflt) + &
-                                 SUM(ztransp(itarg0:itarg0+lowest_free_ray-1) &
-                                             * vffrac(itarg0:itarg0+lowest_free_ray-1))
-
+               skyvf(isurflt) = skyvf(isurflt) +                    &
+                                SUM(vffrac(itarg0:itarg1),          &
+                                    MASK=(itarget(itarg0:itarg1)<0))
+               skyvft(isurflt) = skyvft(isurflt) +                                   &
+                                 SUM(ztransp(itarg0:itarg1) * vffrac(itarg0:itarg1), &
+                                     MASK=(itarget(itarg0:itarg1)<0))
+!
 !--            Save direct solar transparency
                j = MODULO(NINT(azmid/                                          &
                                (2._wp*pi)*raytrace_discrete_azims-.5_wp, iwp), &
                           raytrace_discrete_azims)
-
-               DO k = 1, raytrace_discrete_elevs/2
-                  i = dsidir_rev(k-1, j)
-                  IF ( i /= -1  .AND.  k <= lowest_free_ray )  &
-                     dsitrans(isurflt, i) = ztransp(itarg0+k-1)
-               ENDDO
-
+!
+!--            For down direction there is no direct irradiance, otherwise az0=0
+               IF ( td /= idown )  THEN
+                  DO  k = 1, raytrace_discrete_elevs/2
+                     i = dsidir_rev(k-1, j)
+                     IF ( i /= -1  .AND.  itarget(itarg0+k-1) < 0 )  &
+                        dsitrans(isurflt, i) = ztransp(itarg0+k-1)
+                  ENDDO
+               ENDIF
 !
 !--            Advance itarget indices
@@ -7820,5 +7766,5 @@
                itarg0 = 1
                DO WHILE ( itarg0 <= nzn*naz )
-                  IF ( itarget(itarg0) /= -1 )  EXIT
+                  IF ( itarget(itarg0) >= 0 )  EXIT
                   itarg0 = itarg0 + 1
                ENDDO
@@ -7983,5 +7929,5 @@
               CALL raytrace_2d(ta, yxdir, nzn, zdirs,                                &
                                    -999, -999._wp, vffrac, .FALSE., .FALSE., .TRUE., &
-                                   lowest_free_ray, ztransp, itarget)
+                                   ztransp, itarget)
 
 !--           Save direct solar transparency
@@ -7991,5 +7937,5 @@
               DO  k = 1, raytrace_discrete_elevs/2
                  i = dsidir_rev(k-1, j)
-                 IF ( i /= -1  .AND.  k <= lowest_free_ray ) &
+                 IF ( i /= -1  .AND.  itarget(k) < 0 ) &
                     dsitransc(ipcgb, i) = ztransp(k)
               ENDDO
@@ -8010,6 +7956,5 @@
            zns = pi / REAL(nzn, wp)
            ALLOCATE ( zdirs(1:nzn), zcent(1:nzn), zbdry(0:nzn), vffrac(1:nzn*naz), vffrac0(1:nzn), &
-                      ztransp(1:nzn*naz), itarget(1:nzn*naz) )   !FIXME allocate itarget only
-                                                                 !in case of rad_angular_discretization
+                      ztransp(1:nzn*naz), itarget(1:nzn*naz) )
 
            zcent(:) = (/( zn0+(REAL(izn,wp)-.5_wp)*zns, izn=1, nzn )/)
@@ -8059,14 +8004,14 @@
                  CALL raytrace_2d(ta, yxdir, nzn, zdirs,                         &
                                   -999, -999._wp, vffrac(itarg0:itarg1), .TRUE., &
-                                  .FALSE., .TRUE., lowest_free_ray,              &
-                                  ztransp(itarg0:itarg1),                        &
+                                  .FALSE., .TRUE., ztransp(itarg0:itarg1),       &
                                   itarget(itarg0:itarg1))
 
 !--              Sky view factors for MRT
                  mrtsky(imrt) = mrtsky(imrt) + &
-                                  SUM(vffrac(itarg0:itarg0+lowest_free_ray-1))
+                                  SUM(vffrac(itarg0:itarg1),          &
+                                      MASK=(itarget(itarg0:itarg1)<0))
                  mrtskyt(imrt) = mrtskyt(imrt) + &
-                                   SUM(ztransp(itarg0:itarg0+lowest_free_ray-1) &
-                                               * vffrac(itarg0:itarg0+lowest_free_ray-1))
+                                   SUM(ztransp(itarg0:itarg1) * vffrac(itarg0:itarg1), &
+                                       MASK=(itarget(itarg0:itarg1)<0))
 !--              Direct solar transparency for MRT
                  j = MODULO(NINT(azmid/                                         &
@@ -8075,5 +8020,5 @@
                  DO  k = 1, raytrace_discrete_elevs/2
                     i = dsidir_rev(k-1, j)
-                    IF ( i /= -1  .AND.  k <= lowest_free_ray ) &
+                    IF ( i /= -1  .AND.  itarget(itarg0+k-1) < 0 ) &
                        mrtdsit(imrt, i) = ztransp(itarg0+k-1)
                  ENDDO
@@ -8092,5 +8037,5 @@
               itarg0 = 1
               DO WHILE ( itarg0 <= nzn*naz )
-                 IF ( itarget(itarg0) /= -1 )  EXIT
+                 IF ( itarget(itarg0) >= 0 )  EXIT
                  itarg0 = itarg0 + 1
               ENDDO
@@ -8183,5 +8128,5 @@
 
 !--     deallocate temporary global arrays
-        DEALLOCATE(nzterr)
+        DEALLOCATE( nzterrt, nzterrb )
 
         IF ( plant_canopy )  THEN
@@ -8607,5 +8552,5 @@
 !--             in the array nzterr and plantt and id of the coresponding processor
                 CALL radiation_calc_global_offset( box(3), box(2), 0, 1, offs_glob=ig )
-                IF ( box(1) <= nzterr(ig) )  THEN
+                IF ( box(1) <= nzterrb(ig) )  THEN
                     visible = .FALSE.
                     RETURN
@@ -8703,5 +8648,5 @@
    SUBROUTINE raytrace_2d(origin, yxdir, nrays, zdirs, iorig, aorig, vffrac,  &
                               calc_svf, create_csf, skip_1st_pcb,             &
-                              lowest_free_ray, transparency, itarget)
+                              transparency, itarget)
       IMPLICIT NONE
 
@@ -8716,11 +8661,14 @@
       LOGICAL, INTENT(in)                    ::  create_csf    !< whether to create canopy sink factors
       LOGICAL, INTENT(in)                    ::  skip_1st_pcb  !< whether to skip first plant canopy box during raytracing
-      INTEGER(iwp), INTENT(out)              ::  lowest_free_ray !< index into zdirs
       REAL(wp), DIMENSION(nrays), INTENT(OUT) ::  transparency !< transparencies of zdirs paths
-      INTEGER(iwp), DIMENSION(nrays), INTENT(OUT) ::  itarget  !< global indices of target faces for zdirs
-
-      INTEGER(iwp), DIMENSION(nrays)         ::  target_procs
-      REAL(wp)                               ::  horizon       !< highest horizon found after raytracing (z/hdist)
+      INTEGER(iwp), DIMENSION(nrays), INTENT(OUT) ::  itarget  !< global indices of target faces for zdirs or <0 for sky
+
+      LOGICAL                                ::  is_mixed_col  !< whether current column contains full-3D geometry
+      INTEGER(iwp)                           ::  lowest_free_ray  !< index into zdirs
+      INTEGER(iwp)                           ::  lowest_mixed_ray !< index into zdirs
+      REAL(wp)                               ::  full_horizon  !< highest full horizon found after raytracing (z/hdist)
       INTEGER(iwp)                           ::  i, k, l, d
+      INTEGER(iwp)                           ::  iray         !< index into zdirs
+      INTEGER(iwp)                           ::  isurf        !< index into surf(l)
       INTEGER(iwp)                           ::  seldim       !< dimension to be incremented
       REAL(wp), DIMENSION(2)                 ::  yxorigin     !< horizontal copy of origin (y,x)
@@ -8729,6 +8677,8 @@
       REAL(wp)                               ::  nextdist     !< end of current crossing
       REAL(wp)                               ::  crmid        !< midpoint of crossing
-      REAL(wp)                               ::  horz_entry   !< horizon at entry to column
-      REAL(wp)                               ::  horz_exit    !< horizon at exit from column
+      REAL(wp)                               ::  horz_full_entry  !< full horizon at entry to column
+      REAL(wp)                               ::  horz_full_exit   !< full horizon at exit from column
+      REAL(wp)                               ::  horz_mixed_entry !< mixed horizon at entry to column
+      REAL(wp)                               ::  horz_mixed_exit  !< mixed horizon at exit from column
       REAL(wp)                               ::  bdydim       !< boundary for current dimension
       REAL(wp), DIMENSION(2)                 ::  crossdist    !< distances to boundary for dimensions
@@ -8738,5 +8688,7 @@
       INTEGER(iwp), DIMENSION(2)             ::  dimdelta     !< dimension direction = +- 1
       INTEGER(iwp)                           ::  ip           !< number of processor where gridbox reside
-      INTEGER(iwp)                           ::  ig           !< 1D index of gridbox in global 2D array
+      INTEGER(iwp)                           ::  ip_last      !< previous number of processor where gridbox reside
+      INTEGER(iwp)                           ::  ig           !< 1D index of grid column in global 2D array
+      INTEGER(iwp)                           ::  ig_last      !< 1D index of previous column in global 2D array
       INTEGER(iwp)                           ::  maxboxes     !< max no of CSF created
       INTEGER(iwp)                           ::  nly          !< maximum  plant canopy height
@@ -8746,5 +8698,5 @@
       INTEGER(iwp)                           ::  zb1
       INTEGER(iwp)                           ::  nz
-      INTEGER(iwp)                           ::  iz
+      INTEGER(iwp)                           ::  kz
       INTEGER(iwp)                           ::  zsgn
       INTEGER(iwp)                           ::  lastdir      !< wall direction before hitting this column
@@ -8755,4 +8707,7 @@
       INTEGER(iwp)                           ::  wcount       !< RMA window item count
       INTEGER(MPI_ADDRESS_KIND)              ::  wdisp        !< RMA window displacement
+      INTEGER(iwp), DIMENSION(:), ALLOCATABLE::  vert_col     !< vertical faces from previous to current column
+      INTEGER(iwp), DIMENSION(:), ALLOCATABLE::  up_col       !< upward oriented surfaces in current column
+      INTEGER(iwp), DIMENSION(:), ALLOCATABLE::  down_col     !< downward oriented surfaces in current column
 #endif
 
@@ -8769,11 +8724,15 @@
       yxorigin(:) = origin(2:3)
       transparency(:) = 1._wp !-- Pre-set the all rays to transparent before reducing
-      horizon = -HUGE(1._wp)
+      full_horizon = -HUGE(1._wp)
+      lowest_mixed_ray = nrays
       lowest_free_ray = nrays
-      IF ( rad_angular_discretization  .AND.  calc_svf )  THEN
+      IF ( rad_angular_discretization )  THEN
          ALLOCATE(target_surfl(nrays))
          target_surfl(:) = -1
          lastdir = -999
          lastcolumn(:) = -999
+         ALLOCATE( vert_col(nz_urban_b:nz_urban_t) )
+         ALLOCATE( up_col(nz_urban_b:nz_urban_t) )
+         ALLOCATE( down_col(nz_urban_b:nz_urban_t) )
       ENDIF
 
@@ -8806,4 +8765,6 @@
       ENDIF
 
+      ip_last = -1
+      ig_last = -1
       lastdist = 0._wp
 
@@ -8839,39 +8800,178 @@
 
 !--         calculate index of the grid with global indices (column(1),column(2))
-!--         in the array nzterr and plantt and id of the coresponding processor
-            CALL radiation_calc_global_offset( column(2), column(1), 0, 1, offs_glob=ig )
-
-            IF ( lastdist == 0._wp )  THEN
-               horz_entry = -HUGE(1._wp)
+!--         in the array nzterrt/b and plantt and id of the coresponding processor
+            CALL radiation_calc_global_offset( column(2), column(1), 0, 1, &
+                                               iproc=ip, offs_glob=ig )
+
+            IF ( ip_last < 0 )  THEN
+               horz_full_entry = -HUGE(1._wp)
+               horz_mixed_entry = -HUGE(1._wp)
             ELSE
-               horz_entry = (REAL(nzterr(ig), wp) + .5_wp - origin(1)) / lastdist
+               horz_full_entry = (REAL(nzterrb(ig), wp) + .5_wp - origin(1)) / lastdist
+               horz_mixed_entry = (REAL(nzterrt(ig), wp) + .5_wp - origin(1)) / lastdist
             ENDIF
-            horz_exit = (REAL(nzterr(ig), wp) + .5_wp - origin(1)) / nextdist
-
-            IF ( rad_angular_discretization  .AND.  calc_svf )  THEN
-!
-!--            Identify vertical obstacles hit by rays in current column
-               DO WHILE ( lowest_free_ray > 0 )
-                  IF ( zdirs(lowest_free_ray) > horz_entry )  EXIT
+            horz_full_exit = (REAL(nzterrb(ig), wp) + .5_wp - origin(1)) / nextdist
+            horz_mixed_exit = (REAL(nzterrt(ig), wp) + .5_wp - origin(1)) / nextdist
+            is_mixed_col = ( nzterrt(ig) /= nzterrb(ig) )
+
+            IF ( rad_angular_discretization )  THEN
+!
+!--            Identify vertical full obstacles hit by rays in current column,
+!--            mixed rays need to be checked if they are already obstructed.
+               DO WHILE ( lowest_mixed_ray > lowest_free_ray )
+                  IF ( zdirs(lowest_mixed_ray) > horz_full_entry )  EXIT
 !
 !--               This may only happen after 1st column, so lastdir and lastcolumn are valid
-                  CALL request_itarget(lastdir,                                         &
-                        CEILING(-0.5_wp + origin(1) + zdirs(lowest_free_ray)*lastdist), &
-                        lastcolumn(1), lastcolumn(2),                                   &
-                        target_surfl(lowest_free_ray), target_procs(lowest_free_ray))
-                  lowest_free_ray = lowest_free_ray - 1
+                  IF ( target_surfl(lowest_mixed_ray) < 0 )  THEN
+                     CALL request_itarget(lastdir,                                                 &
+                           CEILING(-0.5_wp + origin(1) + zdirs(lowest_mixed_ray)*lastdist),        &
+                           lastcolumn(1), lastcolumn(2),                                           &
+                           target_surfl(lowest_mixed_ray))
+                  ENDIF
+                  lowest_mixed_ray = lowest_mixed_ray - 1
                ENDDO
 !
-!--            Identify horizontal obstacles hit by rays in current column
-               DO WHILE ( lowest_free_ray > 0 )
-                  IF ( zdirs(lowest_free_ray) > horz_exit )  EXIT
-                  CALL request_itarget(iup_u, nzterr(ig)+1, column(1), column(2), &
-                                       target_surfl(lowest_free_ray),           &
-                                       target_procs(lowest_free_ray))
-                  lowest_free_ray = lowest_free_ray - 1
+!--            Identify vertical full obstacles hit by rays in current column,
+!--            free rays need no individual checks.
+               DO WHILE ( lowest_mixed_ray > 0 )
+                  IF ( zdirs(lowest_mixed_ray) > horz_full_entry )  EXIT
+!
+!--               This may only happen after 1st column, so lastdir and lastcolumn are valid
+                  CALL request_itarget(lastdir,                                                    &
+                        CEILING(-0.5_wp + origin(1) + zdirs(lowest_mixed_ray)*lastdist),           &
+                        lastcolumn(1), lastcolumn(2),                                              &
+                        target_surfl(lowest_mixed_ray))
+                  lowest_mixed_ray = lowest_mixed_ray - 1
                ENDDO
-            ENDIF
-
-            horizon = MAX(horizon, horz_entry, horz_exit)
+               IF ( lowest_free_ray > lowest_mixed_ray )  lowest_free_ray = lowest_mixed_ray
+!
+!--            Identify targets for vertical mixed obstacles.
+!--            lowest_mixed_ray now points to bottom of vertical mixed obstacles.
+               IF ( is_mixed_col  .AND.  ip_last >= 0 )  THEN
+!
+!--               Load vertical surfaces belonging to previous column
+                  vert_col(:) = -999
+                  DO  isurf = surfg_col_start(ig_last), surfg_col_start(ig_last+1)-1
+                     IF ( surf(id, isurf) == lastdir )  THEN
+                        vert_col(surf(iz, isurf)) = isurf
+                     ENDIF
+                  ENDDO
+!
+!--               Previously mixed rays need to be checked whether they are obstructed
+                  DO  iray = lowest_mixed_ray, lowest_free_ray-1, -1
+                     IF ( zdirs(iray) > horz_mixed_entry )  EXIT
+                     IF ( target_surfl(iray) >= 0 )  CYCLE
+                     target_surfl(iray) = vert_col(CEILING(-0.5_wp +     &
+                                        origin(1) + zdirs(iray)*lastdist)) ! contains -999 if missing surface
+                  ENDDO
+!
+!--               Previously free rays cannot be obstructed yet
+                  DO  iray = lowest_free_ray, 1, -1
+                     IF ( zdirs(iray) > horz_mixed_entry )  THEN
+!
+!--                     Extend mixed rays by raising the lowest_free ray
+!--                     (remains unchanged if this is hit in 1st iteration)
+                        lowest_free_ray = iray
+                        EXIT
+                     ENDIF
+                     target_surfl(iray) = vert_col(CEILING(-0.5_wp +     &
+                                        origin(1) + zdirs(iray)*lastdist)) ! contains -999 if missing surface
+                  ENDDO
+               ENDIF ! end of mixed horizon
+!
+!--            Identify horizontal full obstacles hit by rays in current column,
+!--            mixed rays need to be checked if they are already obstructed.
+               DO WHILE ( lowest_mixed_ray > lowest_free_ray )
+                  IF ( zdirs(lowest_mixed_ray) > horz_full_exit )  EXIT
+                  IF ( target_surfl(lowest_mixed_ray) < 0 )  THEN
+                     CALL request_itarget(iup, nzterrb(ig)+1, column(1), column(2), &
+                                          target_surfl(lowest_mixed_ray))
+                  ENDIF
+                  lowest_mixed_ray = lowest_mixed_ray - 1
+               ENDDO
+!
+!--            Identify horizontal full obstacles hit by rays in current column,
+!--            free rays need no individual checks.
+               DO WHILE ( lowest_mixed_ray > 0 )
+                  IF ( zdirs(lowest_mixed_ray) > horz_full_exit )  EXIT
+                  CALL request_itarget(iup, nzterrb(ig)+1, column(1), column(2), &
+                                       target_surfl(lowest_mixed_ray))
+                  lowest_mixed_ray = lowest_mixed_ray - 1
+               ENDDO
+               IF ( lowest_free_ray > lowest_mixed_ray )  lowest_free_ray = lowest_mixed_ray
+!
+!--            Identify targets for horizontal mixed obstacles.
+!--            lowest_mixed_ray now points _above_ horizontal full obstacles.
+               IF ( is_mixed_col )  THEN
+!
+!--               Load horizontal surfaces corresponding to current column
+                  up_col(:) = -999
+                  down_col(:) = -999
+                  DO  isurf = surfg_col_start(ig), surfg_col_start(ig+1)-1
+                     SELECT CASE ( surf(id, isurf) )
+                     CASE ( iup )
+                        up_col(surf(iz, isurf)) = isurf
+                     CASE ( idown )
+                        down_col(surf(iz, isurf)) = isurf
+                     ENDSELECT
+                  ENDDO
+!
+!--               Previously mixed rays need to be checked whether they are obstructed
+                  DO  iray = lowest_mixed_ray, lowest_free_ray-1, -1
+                     IF ( zdirs(iray) > MAX(horz_mixed_entry, horz_mixed_exit) )  EXIT
+                     IF ( target_surfl(iray) >= 0 )  CYCLE
+                     IF ( zdirs(iray) <= 0._wp )  THEN
+!
+!--                     Downward pointed ray, cycle k down from entry to exit,
+!--                     search for upward oriented faces
+                        DO  k = FLOOR(0.5_wp + origin(1) + zdirs(iray)*lastdist), &
+                               CEILING(-0.5_wp + origin(1) + zdirs(iray)*nextdist)+1, -1
+                           target_surfl(iray) = up_col(k) ! contains -999 if missing surface
+                           IF ( target_surfl(iray) >= 0 )  EXIT
+                        ENDDO
+                     ELSE
+!
+!--                     Upward pointed ray, cycle k up from entry to exit,
+!--                     search for downward oriented faces
+                        DO  k = CEILING(-0.5_wp + origin(1) + zdirs(iray)*lastdist), &
+                               FLOOR(0.5_wp + origin(1) + zdirs(iray)*nextdist)-1
+                           target_surfl(iray) = down_col(k) ! contains -999 if missing surface
+                           IF ( target_surfl(iray) >= 0 )  EXIT
+                        ENDDO
+                     ENDIF
+                  ENDDO
+!
+!--               Previously free rays cannot be obstructed yet
+                  DO  iray = lowest_free_ray, 1, -1
+                     IF ( zdirs(iray) > MAX(horz_mixed_entry, horz_mixed_exit) )  THEN
+!
+!--                     Extend mixed rays by raising the lowest_free ray
+!--                     (remains unchanged if this is hit in 1st iteration)
+                        lowest_free_ray = iray
+                        EXIT
+                     ENDIF
+                     IF ( zdirs(iray) <= 0._wp )  THEN
+!
+!--                     Downward pointed ray, cycle k down from entry to exit,
+!--                     search upward oriented faces
+                        DO  k = FLOOR(0.5_wp + origin(1) + zdirs(iray)*lastdist), &
+                               CEILING(-0.5_wp + origin(1) + zdirs(iray)*nextdist)+1, -1
+                           target_surfl(iray) = up_col(k) ! contains -999 if missing surface
+                           IF ( target_surfl(iray) >= 0 )  EXIT
+                        ENDDO
+                     ELSE
+!
+!--                     Upward pointed ray, cycle k up from entry to exit,
+!--                     search downward oriented faces
+                        DO  k = CEILING(-0.5_wp + origin(1) + zdirs(iray)*lastdist), &
+                               FLOOR(0.5_wp + origin(1) + zdirs(iray)*nextdist)-1
+                           target_surfl(iray) = down_col(k) ! contains -999 if missing surface
+                           IF ( target_surfl(iray) >= 0 )  EXIT
+                        ENDDO
+                     ENDIF
+                  ENDDO
+               ENDIF ! end of mixed horizon
+            ENDIF ! rad_angular_discretization
+            full_horizon = MAX(full_horizon, horz_full_entry, horz_full_exit)
 
             IF ( plant_canopy )  THEN
@@ -8883,22 +8983,24 @@
          IF ( nextdist + eps >= distance )  EXIT
 
-         IF ( rad_angular_discretization  .AND.  calc_svf )  THEN
+         IF ( rad_angular_discretization )  THEN
 !
 !--         Save wall direction of coming building column (= this air column)
             IF ( seldim == 1 )  THEN
                IF ( dimdelta(seldim) == 1 )  THEN
-                  lastdir = isouth_u
+                  lastdir = isouth
                ELSE
-                  lastdir = inorth_u
+                  lastdir = inorth
                ENDIF
             ELSE
                IF ( dimdelta(seldim) == 1 )  THEN
-                  lastdir = iwest_u
+                  lastdir = iwest
                ELSE
-                  lastdir = ieast_u
+                  lastdir = ieast
                ENDIF
             ENDIF
             lastcolumn = column
          ENDIF
+         ip_last = ip
+         ig_last = ig
          lastdist = nextdist
          dimnext(seldim) = dimnext(seldim) + dimdelta(seldim)
@@ -8921,11 +9023,11 @@
 !--               For fixed view resolution, we need plant canopy even for rays
 !--               to opposing surfaces
-                  lowest_lad = nzterr(ig) + 1
+                  lowest_lad = nzterrb(ig) + 1
                ELSE
 !
-!--               We only need LAD for rays directed above horizon (to sky)
+!--               We only need LAD for rays directed above full horizon (to sky)
                   lowest_lad = CEILING( -0.5_wp + origin(1) +            &
-                                    MIN( horizon * rt2_track_dist(i-1),  & ! entry
-                                         horizon * rt2_track_dist(i)   ) ) ! exit
+                                    MIN( full_horizon * rt2_track_dist(i-1),  & ! entry
+                                         full_horizon * rt2_track_dist(i)   ) ) ! exit
                ENDIF
 !
@@ -8970,5 +9072,5 @@
       ENDIF ! plant_canopy
 
-      IF ( rad_angular_discretization  .AND.  calc_svf )  THEN
+      IF ( rad_angular_discretization )  THEN
 #if defined( __parallel )
 !--      wait for all gridsurf requests to complete
@@ -8979,18 +9081,13 @@
          ENDIF
 #endif
-!
-!--      recalculate local surf indices into global ones
-         DO i = 1, nrays
-            IF ( target_surfl(i) == -1 )  THEN
-               itarget(i) = -1
-            ELSE
-               itarget(i) = target_surfl(i) + surfstart(target_procs(i))
-            ENDIF
+         itarget(:) = target_surfl(:)
+         DEALLOCATE( target_surfl )
+      ELSE
+         iray = nrays
+         DO WHILE ( iray > 0  .AND.  zdirs(iray) <= full_horizon )
+            itarget(iray) = HUGE(1_iwp)
+            iray = iray - 1
          ENDDO
-
-         DEALLOCATE( target_surfl )
-
-      ELSE
-         itarget(:) = -1
+         itarget(1:iray) = -1
       ENDIF ! rad_angular_discretization
 
@@ -9041,5 +9138,6 @@
 !--               cause less moiree.
                   IF ( .NOT. rad_angular_discretization )  THEN
-                     IF ( zdirs(k) <= horizon )  CYCLE
+                     ! TODO: remove tradiational discretization (incompatible with full 3D)
+                     IF ( zdirs(k) <= full_horizon )  CYCLE
                   ENDIF
 
@@ -9052,5 +9150,5 @@
                      nz = 2
                      rt2_dist(nz) = SQRT(dxxyy)
-                     iz = CEILING(-.5_wp + zorig, iwp)
+                     kz = CEILING(-.5_wp + zorig, iwp)
                   ELSE
                      zexit = MIN(MAX(origin(1) + zdirs(k) * rt2_track_dist(i), zbottom), ztop)
@@ -9062,10 +9160,10 @@
                      qdist = rt2_dist(nz) / (zexit-zorig)
                      rt2_dist(2:nz-1) = (/( ((REAL(l, wp) + .5_wp) * zsgn - zorig) * qdist , l = zb0, zb1 )/)
-                     iz = zb0 * zsgn
+                     kz = zb0 * zsgn
                   ENDIF
 
                   DO  l = 2, nz
-                     IF ( rt2_track_lad(iz, i) > 0._wp )  THEN
-                        curtrans = exp(-ext_coef * rt2_track_lad(iz, i) * (rt2_dist(l)-rt2_dist(l-1)))
+                     IF ( rt2_track_lad(kz, i) > 0._wp )  THEN
+                        curtrans = exp(-ext_coef * rt2_track_lad(kz, i) * (rt2_dist(l)-rt2_dist(l-1)))
 
                         IF ( create_csf )  THEN
@@ -9074,5 +9172,5 @@
                            acsf(ncsfl)%itx = rt2_track(2,i)
                            acsf(ncsfl)%ity = rt2_track(1,i)
-                           acsf(ncsfl)%itz = iz
+                           acsf(ncsfl)%itz = kz
                            acsf(ncsfl)%isurfs = iorig
                            acsf(ncsfl)%rcvf = (1._wp - curtrans)*transparency(k)*vffrac(k)
@@ -9081,10 +9179,11 @@
                         transparency(k) = transparency(k) * curtrans
                      ENDIF
-                     iz = iz + zsgn
+                     kz = kz + zsgn
                   ENDDO ! l = 1, nz - 1
                ENDDO ! k = 1, nrays
             ENDDO ! i = 1, ntrack
-
-            transparency(1:lowest_free_ray) = 1._wp !-- Reset rays above horizon to transparent (see NOTE 6778)
+!
+!--         Reset rays above horizon to transparent (see NOTE 6778)
+            WHERE( itarget < 0 )  transparency = 1.0_wp
          ENDIF
 
@@ -9098,5 +9197,5 @@
 !
 !--            See NOTE 6778 above
-               IF ( zdirs(k) <= horizon )  CYCLE
+               IF ( itarget(k) >= 0 )  CYCLE
 
                zexit = origin(1) + zdirs(k) * rt2_track_dist(i-1)
@@ -9108,5 +9207,5 @@
                   nz = 2
                   rt2_dist(nz) = SQRT(dxxyy)
-                  iz = NINT(zexit, iwp)
+                  kz = NINT(zexit, iwp)
                ELSE
                   zorig = MIN(MAX(origin(1) + zdirs(k) * rt2_track_dist(i), zbottom), ztop)
@@ -9118,10 +9217,10 @@
                   qdist = rt2_dist(nz) / (zexit-zorig)
                   rt2_dist(2:nz-1) = (/( ((REAL(l, wp) + .5_wp) * zsgn - zorig) * qdist , l = zb0, zb1 )/)
-                  iz = zb0 * zsgn
+                  kz = zb0 * zsgn
                ENDIF
 
                DO  l = 2, nz
-                  IF ( rt2_track_lad(iz, i) > 0._wp )  THEN
-                     curtrans = exp(-ext_coef * rt2_track_lad(iz, i) * (rt2_dist(l)-rt2_dist(l-1)))
+                  IF ( rt2_track_lad(kz, i) > 0._wp )  THEN
+                     curtrans = exp(-ext_coef * rt2_track_lad(kz, i) * (rt2_dist(l)-rt2_dist(l-1)))
 
                      IF ( create_csf )  THEN
@@ -9130,6 +9229,5 @@
                         acsf(ncsfl)%itx = rt2_track(2,i)
                         acsf(ncsfl)%ity = rt2_track(1,i)
-                        acsf(ncsfl)%itz = iz
-                        IF ( itarget(k) /= -1 ) STOP 1 !FIXME remove after test
+                        acsf(ncsfl)%itz = kz
                         acsf(ncsfl)%isurfs = -1
                         acsf(ncsfl)%rcvf = (1._wp - curtrans)*transparency(k)*aorig*vffrac(k)
@@ -9138,5 +9236,5 @@
                      transparency(k) = transparency(k) * curtrans
                   ENDIF
-                  iz = iz + zsgn
+                  kz = kz + zsgn
                ENDDO ! l = 1, nz - 1
             ENDDO ! k = 1, nrays
@@ -9144,20 +9242,11 @@
       ENDIF ! plant_canopy
 
-      IF ( .NOT. (rad_angular_discretization  .AND.  calc_svf) )  THEN
-!
-!--      Just update lowest_free_ray according to horizon
-         DO WHILE ( lowest_free_ray > 0 )
-            IF ( zdirs(lowest_free_ray) > horizon )  EXIT
-            lowest_free_ray = lowest_free_ray - 1
-         ENDDO
-      ENDIF
-
    CONTAINS
 
-      SUBROUTINE request_itarget( d, z, y, x, isurfl, iproc )
+      SUBROUTINE request_itarget( d, z, y, x, isurfl )
 
          INTEGER(iwp), INTENT(in)            ::  d, z, y, x
          INTEGER(iwp), TARGET, INTENT(out)   ::  isurfl
-         INTEGER(iwp), INTENT(out)           ::  iproc
+         INTEGER(iwp)                        ::  iproc
 
 #if defined( __parallel )
@@ -9181,5 +9270,4 @@
 !--      set index target_surfl(i)
          isurfl = gridsurf(d,z,y,x)
-         iproc  = 0  ! required to avoid compile error about unused variable in serial mode
 #endif
 
@@ -9299,39 +9387,34 @@
 
 !-- first check: are the two surfaces directed in the same direction
-        IF ( (d==iup_u  .OR.  d==iup_l )                             &
-             .AND. (d2==iup_u  .OR. d2==iup_l) ) RETURN
-        IF ( (d==isouth_u  .OR.  d==isouth_l ) &
-             .AND.  (d2==isouth_u  .OR.  d2==isouth_l) ) RETURN
-        IF ( (d==inorth_u  .OR.  d==inorth_l ) &
-             .AND.  (d2==inorth_u  .OR.  d2==inorth_l) ) RETURN
-        IF ( (d==iwest_u  .OR.  d==iwest_l )     &
-             .AND.  (d2==iwest_u  .OR.  d2==iwest_l ) ) RETURN
-        IF ( (d==ieast_u  .OR.  d==ieast_l )     &
-             .AND.  (d2==ieast_u  .OR.  d2==ieast_l ) ) RETURN
+        IF ( d==iup     .AND.  d2==iup    ) RETURN
+        IF ( d==isouth  .AND.  d2==isouth ) RETURN
+        IF ( d==inorth  .AND.  d2==inorth ) RETURN
+        IF ( d==iwest   .AND.  d2==iwest  ) RETURN
+        IF ( d==ieast   .AND.  d2==ieast  ) RETURN
 
 !-- second check: are surfaces facing away from each other
         SELECT CASE (d)
-            CASE (iup_u, iup_l)                     !< upward facing surfaces
+            CASE (iup)                  !< upward facing surfaces
                 IF ( z2 < z ) RETURN
-            CASE (isouth_u, isouth_l)               !< southward facing surfaces
+            CASE (isouth)               !< southward facing surfaces
                 IF ( y2 > y ) RETURN
-            CASE (inorth_u, inorth_l)               !< northward facing surfaces
+            CASE (inorth)               !< northward facing surfaces
                 IF ( y2 < y ) RETURN
-            CASE (iwest_u, iwest_l)                 !< westward facing surfaces
+            CASE (iwest)                !< westward facing surfaces
                 IF ( x2 > x ) RETURN
-            CASE (ieast_u, ieast_l)                 !< eastward facing surfaces
+            CASE (ieast)                !< eastward facing surfaces
                 IF ( x2 < x ) RETURN
         END SELECT
 
         SELECT CASE (d2)
-            CASE (iup_u)                            !< ground, roof
+            CASE (iup)                  !< ground, roof
                 IF ( z < z2 ) RETURN
-            CASE (isouth_u, isouth_l)               !< south facing
+            CASE (isouth)               !< south facing
                 IF ( y > y2 ) RETURN
-            CASE (inorth_u, inorth_l)               !< north facing
+            CASE (inorth)               !< north facing
                 IF ( y < y2 ) RETURN
-            CASE (iwest_u, iwest_l)                 !< west facing
+            CASE (iwest)                !< west facing
                 IF ( x > x2 ) RETURN
-            CASE (ieast_u, ieast_l)                 !< east facing
+            CASE (ieast)                !< east facing
                 IF ( x < x2 ) RETURN
             CASE (-1)
@@ -9938,5 +10021,5 @@
     INTEGER(iwp) ::  l, m !< index of current surface element
     
-    INTEGER(iwp)                                       :: ids, idsint_u, idsint_l, isurf
+    INTEGER(iwp)                                       :: ids, idsint, isurf
     CHARACTER(LEN=varnamelength)                       :: var
 
@@ -9951,6 +10034,5 @@
            IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) )  THEN
                ids = i
-               idsint_u = dirint_u(ids)
-               idsint_l = dirint_l(ids)
+               idsint = dirint(ids)
                var = var(:k-j)
                EXIT
@@ -10399,6 +10481,6 @@
           CASE ( 'rtm_rad_net' )
 !--           array of complete radiation balance
-              DO isurf = dirstart(ids), dirend(ids)
-                 IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                 IF ( surfl(id,isurf) == idsint )  THEN
                     surfradnet_av(isurf) = surfradnet_av(isurf) +               &
                                            surfinsw(isurf) - surfoutsw(isurf) + &
@@ -10409,6 +10491,6 @@
           CASE ( 'rtm_rad_insw' )
 !--           array of sw radiation falling to surface after i-th reflection
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfinsw_av(isurf) = surfinsw_av(isurf) + surfinsw(isurf)
                   ENDIF
@@ -10417,6 +10499,6 @@
           CASE ( 'rtm_rad_inlw' )
 !--           array of lw radiation falling to surface after i-th reflection
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfinlw_av(isurf) = surfinlw_av(isurf) + surfinlw(isurf)
                   ENDIF
@@ -10425,6 +10507,6 @@
           CASE ( 'rtm_rad_inswdir' )
 !--           array of direct sw radiation falling to surface from sun
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfinswdir_av(isurf) = surfinswdir_av(isurf) + surfinswdir(isurf)
                   ENDIF
@@ -10433,6 +10515,6 @@
           CASE ( 'rtm_rad_inswdif' )
 !--           array of difusion sw radiation falling to surface from sky and borders of the domain
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfinswdif_av(isurf) = surfinswdif_av(isurf) + surfinswdif(isurf)
                   ENDIF
@@ -10441,6 +10523,6 @@
           CASE ( 'rtm_rad_inswref' )
 !--           array of sw radiation falling to surface from reflections
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfinswref_av(isurf) = surfinswref_av(isurf) + surfinsw(isurf) - &
                                           surfinswdir(isurf) - surfinswdif(isurf)
@@ -10451,6 +10533,6 @@
           CASE ( 'rtm_rad_inlwdif' )
 !--           array of sw radiation falling to surface after i-th reflection
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfinlwdif_av(isurf) = surfinlwdif_av(isurf) + surfinlwdif(isurf)
                   ENDIF
@@ -10459,6 +10541,6 @@
           CASE ( 'rtm_rad_inlwref' )
 !--           array of lw radiation falling to surface from reflections
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfinlwref_av(isurf) = surfinlwref_av(isurf) + &
                                           surfinlw(isurf) - surfinlwdif(isurf)
@@ -10468,6 +10550,6 @@
           CASE ( 'rtm_rad_outsw' )
 !--           array of sw radiation emitted from surface after i-th reflection
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfoutsw_av(isurf) = surfoutsw_av(isurf) + surfoutsw(isurf)
                   ENDIF
@@ -10476,6 +10558,6 @@
           CASE ( 'rtm_rad_outlw' )
 !--           array of lw radiation emitted from surface after i-th reflection
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfoutlw_av(isurf) = surfoutlw_av(isurf) + surfoutlw(isurf)
                   ENDIF
@@ -10484,6 +10566,6 @@
           CASE ( 'rtm_rad_ressw' )
 !--           array of residua of sw radiation absorbed in surface after last reflection
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfins_av(isurf) = surfins_av(isurf) + surfins(isurf)
                   ENDIF
@@ -10492,6 +10574,6 @@
           CASE ( 'rtm_rad_reslw' )
 !--           array of residua of lw radiation absorbed in surface after last reflection
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfinl_av(isurf) = surfinl_av(isurf) + surfinl(isurf)
                   ENDIF
@@ -10696,6 +10778,6 @@
           CASE ( 'rtm_rad_net' )
 !--           array of complete radiation balance
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfradnet_av(isurf) = surfradnet_av(isurf) / REAL( average_count_3d, kind=wp )
                   ENDIF
@@ -10704,6 +10786,6 @@
           CASE ( 'rtm_rad_insw' )
 !--           array of sw radiation falling to surface after i-th reflection
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfinsw_av(isurf) = surfinsw_av(isurf) / REAL( average_count_3d, kind=wp )
                   ENDIF
@@ -10712,6 +10794,6 @@
           CASE ( 'rtm_rad_inlw' )
 !--           array of lw radiation falling to surface after i-th reflection
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfinlw_av(isurf) = surfinlw_av(isurf) / REAL( average_count_3d, kind=wp )
                   ENDIF
@@ -10720,6 +10802,6 @@
           CASE ( 'rtm_rad_inswdir' )
 !--           array of direct sw radiation falling to surface from sun
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfinswdir_av(isurf) = surfinswdir_av(isurf) / REAL( average_count_3d, kind=wp )
                   ENDIF
@@ -10728,6 +10810,6 @@
           CASE ( 'rtm_rad_inswdif' )
 !--           array of difusion sw radiation falling to surface from sky and borders of the domain
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfinswdif_av(isurf) = surfinswdif_av(isurf) / REAL( average_count_3d, kind=wp )
                   ENDIF
@@ -10736,6 +10818,6 @@
           CASE ( 'rtm_rad_inswref' )
 !--           array of sw radiation falling to surface from reflections
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfinswref_av(isurf) = surfinswref_av(isurf) / REAL( average_count_3d, kind=wp )
                   ENDIF
@@ -10744,6 +10826,6 @@
           CASE ( 'rtm_rad_inlwdif' )
 !--           array of sw radiation falling to surface after i-th reflection
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfinlwdif_av(isurf) = surfinlwdif_av(isurf) / REAL( average_count_3d, kind=wp )
                   ENDIF
@@ -10752,6 +10834,6 @@
           CASE ( 'rtm_rad_inlwref' )
 !--           array of lw radiation falling to surface from reflections
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfinlwref_av(isurf) = surfinlwref_av(isurf) / REAL( average_count_3d, kind=wp )
                   ENDIF
@@ -10760,6 +10842,6 @@
           CASE ( 'rtm_rad_outsw' )
 !--           array of sw radiation emitted from surface after i-th reflection
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfoutsw_av(isurf) = surfoutsw_av(isurf) / REAL( average_count_3d, kind=wp )
                   ENDIF
@@ -10768,6 +10850,6 @@
           CASE ( 'rtm_rad_outlw' )
 !--           array of lw radiation emitted from surface after i-th reflection
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfoutlw_av(isurf) = surfoutlw_av(isurf) / REAL( average_count_3d, kind=wp )
                   ENDIF
@@ -10776,6 +10858,6 @@
           CASE ( 'rtm_rad_ressw' )
 !--           array of residua of sw radiation absorbed in surface after last reflection
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfins_av(isurf) = surfins_av(isurf) / REAL( average_count_3d, kind=wp )
                   ENDIF
@@ -10784,6 +10866,6 @@
           CASE ( 'rtm_rad_reslw' )
 !--           array of residua of lw radiation absorbed in surface after last reflection
-              DO isurf = dirstart(ids), dirend(ids)
-                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+              DO isurf = 1, nsurfl
+                  IF ( surfl(id,isurf) == idsint )  THEN
                       surfinl_av(isurf) = surfinl_av(isurf) / REAL( average_count_3d, kind=wp )
                   ENDIF
@@ -11354,5 +11436,5 @@
 
     CHARACTER (len=varnamelength)                   :: var, surfid
-    INTEGER(iwp)                                    :: ids,idsint_u,idsint_l,isurf,isvf,isurfs,isurflt,ipcgb
+    INTEGER(iwp)                                    :: ids,idsint,isurf,isvf,isurfs,isurflt,ipcgb
     INTEGER(iwp)                                    :: is, js, ks, istat
 
@@ -11378,6 +11460,5 @@
            IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) )  THEN
               ids = i
-              idsint_u = dirint_u(ids)
-              idsint_l = dirint_l(ids)
+              idsint = dirint(ids)
               var = var(:k-j)
               EXIT
@@ -11596,6 +11677,6 @@
       CASE ( 'rtm_rad_net' )
 !--     array of complete radiation balance
-         DO isurf = dirstart(ids), dirend(ids)
-            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+         DO isurf = 1, nsurfl
+            IF ( surfl(id,isurf) == idsint )  THEN
                IF ( av == 0 )  THEN
                   local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = &
@@ -11609,6 +11690,6 @@
       CASE ( 'rtm_rad_insw' )
 !--      array of sw radiation falling to surface after i-th reflection
-         DO isurf = dirstart(ids), dirend(ids)
-            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+         DO isurf = 1, nsurfl
+            IF ( surfl(id,isurf) == idsint )  THEN
                IF ( av == 0 )  THEN
                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinsw(isurf)
@@ -11621,6 +11702,6 @@
       CASE ( 'rtm_rad_inlw' )
 !--      array of lw radiation falling to surface after i-th reflection
-         DO isurf = dirstart(ids), dirend(ids)
-            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+         DO isurf = 1, nsurfl
+            IF ( surfl(id,isurf) == idsint )  THEN
                IF ( av == 0 )  THEN
                   local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinlw(isurf)
@@ -11633,6 +11714,6 @@
       CASE ( 'rtm_rad_inswdir' )
 !--      array of direct sw radiation falling to surface from sun
-         DO isurf = dirstart(ids), dirend(ids)
-            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+         DO isurf = 1, nsurfl
+            IF ( surfl(id,isurf) == idsint )  THEN
                IF ( av == 0 )  THEN
                   local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinswdir(isurf)
@@ -11645,6 +11726,6 @@
       CASE ( 'rtm_rad_inswdif' )
 !--      array of difusion sw radiation falling to surface from sky and borders of the domain
-         DO isurf = dirstart(ids), dirend(ids)
-            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+         DO isurf = 1, nsurfl
+            IF ( surfl(id,isurf) == idsint )  THEN
                IF ( av == 0 )  THEN
                   local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinswdif(isurf)
@@ -11657,6 +11738,6 @@
       CASE ( 'rtm_rad_inswref' )
 !--      array of sw radiation falling to surface from reflections
-         DO isurf = dirstart(ids), dirend(ids)
-            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+         DO isurf = 1, nsurfl
+            IF ( surfl(id,isurf) == idsint )  THEN
                IF ( av == 0 )  THEN
                   local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = &
@@ -11670,6 +11751,6 @@
       CASE ( 'rtm_rad_inlwdif' )
 !--      array of difusion lw radiation falling to surface from sky and borders of the domain
-         DO isurf = dirstart(ids), dirend(ids)
-            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+         DO isurf = 1, nsurfl
+            IF ( surfl(id,isurf) == idsint )  THEN
                IF ( av == 0 )  THEN
                   local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinlwdif(isurf)
@@ -11682,6 +11763,6 @@
       CASE ( 'rtm_rad_inlwref' )
 !--      array of lw radiation falling to surface from reflections
-         DO isurf = dirstart(ids), dirend(ids)
-            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+         DO isurf = 1, nsurfl
+            IF ( surfl(id,isurf) == idsint )  THEN
                IF ( av == 0 )  THEN
                   local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinlw(isurf) - surfinlwdif(isurf)
@@ -11694,6 +11775,6 @@
       CASE ( 'rtm_rad_outsw' )
 !--      array of sw radiation emitted from surface after i-th reflection
-         DO isurf = dirstart(ids), dirend(ids)
-            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+         DO isurf = 1, nsurfl
+            IF ( surfl(id,isurf) == idsint )  THEN
                IF ( av == 0 )  THEN
                   local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfoutsw(isurf)
@@ -11706,6 +11787,6 @@
       CASE ( 'rtm_rad_outlw' )
 !--      array of lw radiation emitted from surface after i-th reflection
-         DO isurf = dirstart(ids), dirend(ids)
-            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+         DO isurf = 1, nsurfl
+            IF ( surfl(id,isurf) == idsint )  THEN
                IF ( av == 0 )  THEN
                   local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfoutlw(isurf)
@@ -11718,6 +11799,6 @@
       CASE ( 'rtm_rad_ressw' )
 !--      average of array of residua of sw radiation absorbed in surface after last reflection
-         DO isurf = dirstart(ids), dirend(ids)
-            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+         DO isurf = 1, nsurfl
+            IF ( surfl(id,isurf) == idsint )  THEN
                IF ( av == 0 )  THEN
                   local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfins(isurf)
@@ -11730,6 +11811,6 @@
       CASE ( 'rtm_rad_reslw' )
 !--      average of array of residua of lw radiation absorbed in surface after last reflection
-         DO isurf = dirstart(ids), dirend(ids)
-            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+         DO isurf = 1, nsurfl
+            IF ( surfl(id,isurf) == idsint )  THEN
                IF ( av == 0 )  THEN
                   local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinl(isurf)
@@ -11838,6 +11919,6 @@
 !
 !--      sky view factor
-         DO isurf = dirstart(ids), dirend(ids)
-            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+         DO isurf = 1, nsurfl
+            IF ( surfl(id,isurf) == idsint )  THEN
                local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = skyvf(isurf)
             ENDIF
@@ -11847,6 +11928,6 @@
 !
 !--      sky view factor
-         DO isurf = dirstart(ids), dirend(ids)
-            IF ( surfl(id,isurf) == ids )  THEN
+         DO isurf = 1, nsurfl
+            IF ( surfl(id,isurf) == idsint )  THEN
                local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = skyvft(isurf)
             ENDIF
@@ -11866,5 +11947,5 @@
 
             IF ( surf(ix,isurfs) == is  .AND.  surf(iy,isurfs) == js  .AND. surf(iz,isurfs) == ks  .AND. &
-                 (surfl(id,isurflt) == idsint_u .OR. surfl(id,isurflt) == idsint_l ) ) THEN
+                 surfl(id,isurflt) == idsint ) THEN
 !
 !--            correct source surface
@@ -11876,6 +11957,6 @@
 !
 !--      surface albedo
-         DO isurf = dirstart(ids), dirend(ids)
-            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+         DO isurf = 1, nsurfl
+            IF ( surfl(id,isurf) == idsint )  THEN
                local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = albedo_surf(isurf)
             ENDIF
@@ -11885,6 +11966,6 @@
 !
 !--      surface emissivity, weighted average
-         DO isurf = dirstart(ids), dirend(ids)
-            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
+         DO isurf = 1, nsurfl
+            IF ( surfl(id,isurf) == idsint )  THEN
                local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = emiss_surf(isurf)
             ENDIF
Index: /palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4652)
+++ /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4653)
@@ -362,14 +362,9 @@
                force_radiation_call,                                                               &
                id,                                                                                 &
-               ieast_l,                                                                            &
-               ieast_u,                                                                            &
-               inorth_l,                                                                           &
-               inorth_u,                                                                           &
-               isouth_l,                                                                           &
-               isouth_u,                                                                           &
-               iup_l,                                                                              &
-               iup_u,                                                                              &
-               iwest_l,                                                                            &
-               iwest_u,                                                                            &
+               ieast,                                                                              &
+               inorth,                                                                             &
+               isouth,                                                                             &
+               iup,                                                                                &
+               iwest,                                                                              &
                nz_urban_b,                                                                         &
                nz_urban_t,                                                                         &
@@ -1343,5 +1338,5 @@
     INTEGER(iwp), PARAMETER                            ::  nd = 5                                  !< number of directions
     CHARACTER(LEN=6), DIMENSION(0:nd-1), PARAMETER     ::  dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
-    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER         ::  dirint = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)
+    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER         ::  dirint = (/ iup, isouth, inorth, iwest, ieast /)
 
 
@@ -2460,5 +2455,5 @@
 
     INTEGER(iwp), PARAMETER                            ::  nd = 5  !< number of directions
-    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER         ::  dirint =  (/    iup_u, isouth_u, inorth_u,  iwest_u,  ieast_u /)  !<
+    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER         ::  dirint =  (/      iup,   isouth,   inorth,    iwest,    ieast /)  !<
     INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER         ::  diridx =  (/       -1,        1,        0,        3,        2 /)
                                                            !< index for surf_*_v: 0:3 = (North, South, East, West)
@@ -2531,5 +2526,5 @@
 !
 !--       Array of surface height (z)
-          IF ( idsint == iup_u )  THEN
+          IF ( idsint == iup )  THEN
              DO  m = 1, surf_usm_h%ns
                 i = surf_usm_h%i(m)
@@ -2551,5 +2546,5 @@
 !
 !--       Surface category
-          IF ( idsint == iup_u )  THEN
+          IF ( idsint == iup )  THEN
              DO  m = 1, surf_usm_h%ns
                 i = surf_usm_h%i(m)
@@ -2571,5 +2566,5 @@
 !
 !--       Transmissivity window tiles
-          IF ( idsint == iup_u )  THEN
+          IF ( idsint == iup )  THEN
              DO  m = 1, surf_usm_h%ns
                 i = surf_usm_h%i(m)
@@ -2592,5 +2587,5 @@
 !--       Array of sensible heat flux from surfaces
           IF ( av == 0 )  THEN
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -2609,5 +2604,5 @@
              ENDIF
           ELSE
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -2632,5 +2627,5 @@
 !--       Array of latent heat flux from surfaces
           IF ( av == 0 )  THEN
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -2649,5 +2644,5 @@
              ENDIF
           ELSE
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -2671,5 +2666,5 @@
 !--       Array of latent heat flux from vegetation surfaces
           IF ( av == 0 )  THEN
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -2688,5 +2683,5 @@
              ENDIF
           ELSE
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -2710,5 +2705,5 @@
 !--       Array of latent heat flux from surfaces with liquid
           IF ( av == 0 )  THEN
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -2727,5 +2722,5 @@
              ENDIF
           ELSE
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -2749,5 +2744,5 @@
 !--       Array of heat flux from ground (land, wall, roof)
           IF ( av == 0 )  THEN
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -2766,5 +2761,5 @@
              ENDIF
           ELSE
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -2788,5 +2783,5 @@
 !--       Array of heat flux from window ground (land, wall, roof)
           IF ( av == 0 )  THEN
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -2805,5 +2800,5 @@
              ENDIF
           ELSE
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -2827,5 +2822,5 @@
 !--       Array of heat flux from green ground (land, wall, roof)
           IF ( av == 0 )  THEN
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -2844,5 +2839,5 @@
              ENDIF
           ELSE
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -2866,5 +2861,5 @@
 !--       Array of heat flux from indoor ground (land, wall, roof)
           IF ( av == 0 )  THEN
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -2883,5 +2878,5 @@
              ENDIF
           ELSE
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -2905,5 +2900,5 @@
 !--       Array of heat flux from indoor window ground (land, wall, roof)
           IF ( av == 0 )  THEN
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -2922,5 +2917,5 @@
              ENDIF
           ELSE
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -2944,5 +2939,5 @@
 !--       Surface temperature for surfaces
           IF ( av == 0 )  THEN
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -2961,5 +2956,5 @@
              ENDIF
           ELSE
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -2983,5 +2978,5 @@
 !--       Surface temperature for window surfaces
           IF ( av == 0 )  THEN
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -3001,5 +2996,5 @@
 
           ELSE
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -3025,5 +3020,5 @@
 !--       Surface temperature for green surfaces
           IF ( av == 0 )  THEN
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -3043,5 +3038,5 @@
 
           ELSE
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -3067,5 +3062,5 @@
 !--       Near surface temperature for whole surfaces
           IF ( av == 0 )  THEN
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -3086,5 +3081,5 @@
 
           ELSE
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -3109,5 +3104,5 @@
 !--       Wall temperature for  iwl layer of walls and land
           IF ( av == 0 )  THEN
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -3126,5 +3121,5 @@
              ENDIF
           ELSE
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -3148,5 +3143,5 @@
 !--       Window temperature for iwl layer of walls and land
           IF ( av == 0 )  THEN
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -3165,5 +3160,5 @@
              ENDIF
           ELSE
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -3187,5 +3182,5 @@
 !--       Green temperature for  iwl layer of walls and land
           IF ( av == 0 )  THEN
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -3204,5 +3199,5 @@
              ENDIF
           ELSE
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -3226,5 +3221,5 @@
 !--       Soil water content for  iwl layer of walls and land
           IF ( av == 0 )  THEN
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
@@ -3237,5 +3232,5 @@
              ENDIF
           ELSE
-             IF ( idsint == iup_u )  THEN
+             IF ( idsint == iup )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
