Index: palm/trunk/SOURCE/init_grid.f90
===================================================================
--- palm/trunk/SOURCE/init_grid.f90	(revision 4360)
+++ palm/trunk/SOURCE/init_grid.f90	(revision 4386)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Allocation statements, comments, naming of variables revised and _wp added to 
+! real type values
+! 
+! 4360 2020-01-07 11:25:50Z suehring
 ! Revise error messages for generic tunnel setup.
 ! 
@@ -199,5 +203,7 @@
 !
 !-- Allocate grid arrays
-    ALLOCATE( x(0:nx), xu(0:nx) )
+    ALLOCATE( x(0:nx) )
+    ALLOCATE( xu(0:nx) )
+    
     DO i = 0, nx
        xu(i) = i * dx
@@ -205,5 +211,7 @@
     ENDDO
 
-    ALLOCATE( y(0:ny), yv(0:ny) )
+    ALLOCATE( y(0:ny) )
+    ALLOCATE( yv(0:ny) )
+    
     DO j = 0, ny
        yv(j) = j * dy
@@ -211,11 +219,15 @@
     ENDDO
 
-!
-!-- Allocate grid arrays
-    ALLOCATE( ddzu(1:nzt+1), ddzw(1:nzt+1), dd2zu(1:nzt), dzu(1:nzt+1),        &
-              dzw(1:nzt+1), zu(nzb:nzt+1), zw(nzb:nzt+1) )
-
-!
-!-- Compute height of u-levels from constant grid length and dz stretch factors
+    ALLOCATE( ddzu(1:nzt+1) )
+    ALLOCATE( ddzw(1:nzt+1) )
+    ALLOCATE( dd2zu(1:nzt) )
+    ALLOCATE( dzu(1:nzt+1) )
+    ALLOCATE( dzw(1:nzt+1) )
+    ALLOCATE( zu(nzb:nzt+1) )
+    ALLOCATE( zw(nzb:nzt+1) )
+
+!
+!-- For constructing an appropriate grid, the vertical grid spacing dz has to
+!-- be specified with a non-negative value in the parameter file
     IF ( dz(1) == -1.0_wp )  THEN
        message_string = 'missing dz'
@@ -248,11 +260,11 @@
 
 !
-!-- The number of specified end levels +1 has to be the same than the number 
+!-- The number of specified end levels +1 has to be the same as the number 
 !-- of specified dz values
     IF ( number_dz /= number_stretch_level_end + 1 ) THEN
-       WRITE( message_string, * ) 'The number of values for dz = ',         &
-                                   number_dz, 'has to be the same than& ',  &
-                                   'the number of values for ',             &
-                                   'dz_stretch_level_end + 1 = ',           &
+       WRITE( message_string, * ) 'The number of values for dz = ',            &
+                                   number_dz, 'has to be the same as& ',       &
+                                   'the number of values for ',                &
+                                   'dz_stretch_level_end + 1 = ',              &
                                    number_stretch_level_end+1
           CALL message( 'init_grid', 'PA0156', 1, 2, 0, 6, 0 )
@@ -260,25 +272,25 @@
     
 !
-!--    The number of specified start levels has to be the same or one less than 
-!--    the number of specified dz values
-    IF ( number_dz /= number_stretch_level_start + 1 .AND.                  &
+!-- The number of specified start levels has to be the same or one less than 
+!-- the number of specified dz values
+    IF ( number_dz /= number_stretch_level_start + 1 .AND.                     &
          number_dz /= number_stretch_level_start ) THEN
-       WRITE( message_string, * ) 'The number of values for dz = ',         &
-                                   number_dz, 'has to be the same or one ', &
-                                   'more than& the number of values for ',  &
-                                   'dz_stretch_level_start = ',             &
+       WRITE( message_string, * ) 'The number of values for dz = ',            &
+                                   number_dz, 'has to be the same as or one ', &
+                                   'more than& the number of values for ',     &
+                                   'dz_stretch_level_start = ',                &
                                    number_stretch_level_start
           CALL message( 'init_grid', 'PA0211', 1, 2, 0, 6, 0 )
     ENDIF
     
-!--    The number of specified start levels has to be the same or one more than 
-!--    the number of specified end levels
-    IF ( number_stretch_level_start /= number_stretch_level_end + 1 .AND.   &
+!-- The number of specified start levels has to be the same or one more than 
+!-- the number of specified end levels
+    IF ( number_stretch_level_start /= number_stretch_level_end + 1 .AND.      &
          number_stretch_level_start /= number_stretch_level_end ) THEN
-       WRITE( message_string, * ) 'The number of values for ',              &
-                                  'dz_stretch_level_start = ',              &
-                                   dz_stretch_level_start, 'has to be the ',&
-                                   'same or one more than& the number of ', &
-                                   'values for dz_stretch_level_end = ',    &
+       WRITE( message_string, * ) 'The number of values for ',                 &
+                                  'dz_stretch_level_start = ',                 &
+                                   dz_stretch_level_start, 'has to be the ',   &
+                                   'same or one more than& the number of ',    &
+                                   'values for dz_stretch_level_end = ',       &
                                    number_stretch_level_end
           CALL message( 'init_grid', 'PA0216', 1, 2, 0, 6, 0 )
@@ -306,10 +318,11 @@
 
 !
-!-- Define the vertical grid levels
+!-- Define the vertical grid levels. Start with atmosphere branch
     IF ( .NOT. ocean_mode )  THEN
     
 !
 !--    The stretching region has to be large enough to allow for a smooth
-!--    transition between two different grid spacings
+!--    transition between two different grid spacings. The number 4 is an
+!--    empirical value
        DO n = 1, number_stretch_level_start
           min_dz_stretch_level_end(n) = dz_stretch_level_start(n) +            &
@@ -319,5 +332,5 @@
        IF ( ANY( min_dz_stretch_level_end(1:number_stretch_level_start) >      &
                  dz_stretch_level_end(1:number_stretch_level_start) ) ) THEN
-             message_string= 'Eeach dz_stretch_level_end has to be larger ' // &
+             message_string= 'Each dz_stretch_level_end has to be larger ' //  &
                              'than its corresponding value for &' //           &
                              'dz_stretch_level_start + 4*MAX(dz(n),dz(n+1)) '//&
@@ -327,9 +340,9 @@
        
 !
-!--    Stretching must not be applied within the prandtl_layer 
+!--    Stretching must not be applied within the surface layer 
 !--    (first two grid points). For the default case dz_stretch_level_start 
 !--    is negative. Therefore the absolut value is checked here.
-       IF ( ANY( ABS( dz_stretch_level_start ) < dz(1) * 1.5_wp ) ) THEN
-          WRITE( message_string, * ) 'Eeach dz_stretch_level_start has to be ',&
+       IF ( ANY( ABS( dz_stretch_level_start ) <= dz(1) * 1.5_wp ) ) THEN
+          WRITE( message_string, * ) 'Each dz_stretch_level_start has to be ',&
                                      'larger than ', dz(1) * 1.5
              CALL message( 'init_grid', 'PA0226', 1, 2, 0, 6, 0 )
@@ -338,5 +351,9 @@
 !
 !--    The stretching has to start and end on a grid level. Therefore 
-!--    user-specified values have to ''interpolate'' to the next lowest level
+!--    user-specified values are mapped to the next lowest level. The  
+!--    calculation of the first level is realized differently just because of 
+!--    historical reasons (the advanced/new stretching mechanism was realized  
+!--    in a way that results don't change if the old parameters 
+!--    dz_stretch_level, dz_stretch_factor and dz_max are used)
        IF ( number_stretch_level_start /= 0 ) THEN
           dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) -        &
@@ -358,5 +375,5 @@
           ENDDO
        ENDIF
- 
+
 !
 !--    Determine stretching factor if necessary
@@ -370,5 +387,5 @@
 !--    the first u/v- and w-level (k=0) are defined at same height (z=0). 
 !--    The second u-level (k=1) corresponds to the top of the 
-!--    Prandtl-layer. In case of symmetric boundaries (closed channel flow), 
+!--    surface layer. In case of symmetric boundaries (closed channel flow), 
 !--    the first grid point is always at z=0.
        IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 .OR.                              & 
@@ -390,5 +407,6 @@
 
 !--    The default value of dz_stretch_level_start is negative, thus the first
-!--    condition is always true. Hence, the second condition is necessary.
+!--    condition is true even if no stretching shall be applied. Hence, the 
+!--    second condition is also necessary.
        DO  k = 2, nzt+1-symmetry_flag
           IF ( dz_stretch_level_start(n) <= zu(k-1) .AND.                      &
@@ -452,9 +470,10 @@
        ENDIF
 
-    ELSE
+    ELSE !ocean branch
 
 !
 !--    The stretching region has to be large enough to allow for a smooth
-!--    transition between two different grid spacings
+!--    transition between two different grid spacings. The number 4 is an
+!--    empirical value
        DO n = 1, number_stretch_level_start
           min_dz_stretch_level_end(n) = dz_stretch_level_start(n) -            &
@@ -464,5 +483,5 @@
        IF ( ANY( min_dz_stretch_level_end (1:number_stretch_level_start) <     &
                  dz_stretch_level_end(1:number_stretch_level_start) ) ) THEN
-             message_string= 'Eeach dz_stretch_level_end has to be less ' //   &
+             message_string= 'Each dz_stretch_level_end has to be less ' //   &
                              'than its corresponding value for &' //           &
                              'dz_stretch_level_start - 4*MAX(dz(n),dz(n+1)) '//&
@@ -474,7 +493,7 @@
 !--    Stretching must not be applied close to the surface (last two grid 
 !--    points). For the default case dz_stretch_level_start is negative. 
-       IF ( ANY( dz_stretch_level_start > - dz(1) * 1.5_wp ) ) THEN
-          WRITE( message_string, * ) 'Eeach dz_stretch_level_start has to be ',&
-                                     'less than ', dz(1) * 1.5
+       IF ( ANY( dz_stretch_level_start >= - dz(1) * 1.5_wp ) ) THEN
+          WRITE( message_string, * ) 'Each dz_stretch_level_start has to be ',&
+                                     'less than ', -dz(1) * 1.5
              CALL message( 'init_grid', 'PA0226', 1, 2, 0, 6, 0 )
        ENDIF
@@ -482,5 +501,9 @@
 !
 !--    The stretching has to start and end on a grid level. Therefore 
-!--    user-specified values have to ''interpolate'' to the next highest level
+!--    user-specified values are mapped to the next highest level. The  
+!--    calculation of the first level is realized differently just because of 
+!--    historical reasons (the advanced/new stretching mechanism was realized  
+!--    in a way that results don't change if the old parameters 
+!--    dz_stretch_level, dz_stretch_factor and dz_max are used)
        IF ( number_stretch_level_start /= 0 ) THEN
           dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) +        &
@@ -514,5 +537,5 @@
 !--    below the first w-level (k=0). In case of dirichlet bc the first u- and
 !--    w-level are defined at same height, but staggered from the second level.
-!--    The second u-level (k=1) corresponds to the top of the Prandtl-layer.
+!--    The second u-level (k=1) corresponds to the top of the surface layer.
 !--    z values are negative starting from z=0 (surface)
        zu(nzt+1) =   dz(1) * 0.5_wp
@@ -576,5 +599,5 @@
        ENDIF
 
-    ENDIF
+    ENDIF !End of defining the vertical grid levels
 
 !
@@ -617,5 +640,5 @@
     topo = 0
 !
-!-- Initialize topography by generic topography or read from topography from file.  
+!-- Initialize topography by generic topography or read topography from file.  
     CALL init_topo( topo )
 !
@@ -691,15 +714,15 @@
 !--    Allocate and set the arrays containing the topography height
        IF ( nxr == nx  .AND.  nyn /= ny )  THEN
-          ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn),                             &
-                    zw_w_inner(nxl:nxr+1,nys:nyn) )
+          ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn) )                             
+          ALLOCATE( zw_w_inner(nxl:nxr+1,nys:nyn) )
        ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
-          ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn+1),                             &
-                    zw_w_inner(nxl:nxr,nys:nyn+1) )
+          ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn+1) )                             
+          ALLOCATE( zw_w_inner(nxl:nxr,nys:nyn+1) )
        ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
-          ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn+1),                           &
-                    zw_w_inner(nxl:nxr+1,nys:nyn+1) )
+          ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn+1) ) 
+          ALLOCATE( zw_w_inner(nxl:nxr+1,nys:nyn+1) )
        ELSE
-          ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn),                               &
-                    zw_w_inner(nxl:nxr,nys:nyn) )
+          ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn) )
+          ALLOCATE( zw_w_inner(nxl:nxr,nys:nyn) )
        ENDIF
 
@@ -731,16 +754,16 @@
 !-- Allocate outer and inner index arrays for topography and set
 !-- defaults.                    
-    ALLOCATE( nzb_s_inner(nysg:nyng,nxlg:nxrg),                                &
-              nzb_s_outer(nysg:nyng,nxlg:nxrg),                                &
-              nzb_u_inner(nysg:nyng,nxlg:nxrg),                                &
-              nzb_u_outer(nysg:nyng,nxlg:nxrg),                                &
-              nzb_v_inner(nysg:nyng,nxlg:nxrg),                                &
-              nzb_v_outer(nysg:nyng,nxlg:nxrg),                                &
-              nzb_w_inner(nysg:nyng,nxlg:nxrg),                                &
-              nzb_w_outer(nysg:nyng,nxlg:nxrg),                                &
-              nzb_diff_s_inner(nysg:nyng,nxlg:nxrg),                           &
-              nzb_diff_s_outer(nysg:nyng,nxlg:nxrg),                           &
-              nzb_local(nysg:nyng,nxlg:nxrg),                                  &
-              nzb_tmp(nysg:nyng,nxlg:nxrg) )
+    ALLOCATE( nzb_s_inner(nysg:nyng,nxlg:nxrg) )                               
+    ALLOCATE( nzb_s_outer(nysg:nyng,nxlg:nxrg) )                               
+    ALLOCATE( nzb_u_inner(nysg:nyng,nxlg:nxrg) )                               
+    ALLOCATE( nzb_u_outer(nysg:nyng,nxlg:nxrg) )                               
+    ALLOCATE( nzb_v_inner(nysg:nyng,nxlg:nxrg) )                               
+    ALLOCATE( nzb_v_outer(nysg:nyng,nxlg:nxrg) )                               
+    ALLOCATE( nzb_w_inner(nysg:nyng,nxlg:nxrg) )                               
+    ALLOCATE( nzb_w_outer(nysg:nyng,nxlg:nxrg) )                               
+    ALLOCATE( nzb_diff_s_inner(nysg:nyng,nxlg:nxrg) )                          
+    ALLOCATE( nzb_diff_s_outer(nysg:nyng,nxlg:nxrg) )                          
+    ALLOCATE( nzb_local(nysg:nyng,nxlg:nxrg) )                                 
+    ALLOCATE( nzb_tmp(nysg:nyng,nxlg:nxrg) )
 !
 !-- Initialize 2D-index arrays. Note, these will be removed soon! 
@@ -996,7 +1019,7 @@
     REAL(wp) ::  dz_stretch_factor_array_2(9) = 1.08_wp  !< Array that contains all stretch_factor_2 that belongs to stretch_factor_1 
     
-    REAL(wp), PARAMETER ::  stretch_factor_interval = 1.0E-06  !< interval for sampling possible stretching factors
-    REAL(wp), PARAMETER ::  stretch_factor_lower_limit = 0.88  !< lowest possible stretching factor
-    REAL(wp), PARAMETER ::  stretch_factor_upper_limit = 1.12  !< highest possible stretching factor
+    REAL(wp), PARAMETER ::  stretch_factor_interval = 1.0E-06_wp  !< interval for sampling possible stretching factors
+    REAL(wp), PARAMETER ::  stretch_factor_lower_limit = 0.88_wp  !< lowest possible stretching factor
+    REAL(wp), PARAMETER ::  stretch_factor_upper_limit = 1.12_wp  !< highest possible stretching factor
  
  
@@ -1005,19 +1028,21 @@
     
        iterations = 1
-       stretch_factor_1 = 1.0 
-       stretch_factor_2 = 1.0
-       delta_total_old = 1.0
+       stretch_factor_1 = 1.0_wp 
+       stretch_factor_2 = 1.0_wp
+       delta_total_old = 1.0_wp
        
+!
+!--    First branch for stretching from rough to fine
        IF ( dz(n) > dz(n+1) ) THEN
           DO WHILE ( stretch_factor_1 >= stretch_factor_lower_limit ) 
              
-             stretch_factor_1 = 1.0 - iterations * stretch_factor_interval
-             distance = ABS( dz_stretch_level_end(n) -                   &
-                        dz_stretch_level_start(n) ) 
-             numerator = distance*stretch_factor_1/dz(n) +               &
+             stretch_factor_1 = 1.0_wp - iterations * stretch_factor_interval
+             distance = ABS( dz_stretch_level_end(n) -                         &
+                        dz_stretch_level_start(n) )   
+             numerator = distance*stretch_factor_1/dz(n) +                     &
                          stretch_factor_1 - distance/dz(n)
              
-             IF ( numerator > 0.0 ) THEN
-                l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0
+             IF ( numerator > 0.0_wp ) THEN
+                l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0_wp
                 l_rounded = NINT( l )
                 delta_l = ABS( l_rounded - l ) / l
@@ -1026,6 +1051,6 @@
              stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) )
              
-             delta_stretch_factor = ABS( stretch_factor_1 -              &
-                                         stretch_factor_2 ) /            &
+             delta_stretch_factor = ABS( stretch_factor_1 -                    &
+                                         stretch_factor_2 ) /                  &
                                     stretch_factor_2
              
@@ -1033,8 +1058,8 @@
 
 !
-!--                stretch_factor_1 is taken to guarantee that the stretching
-!--                procedure ends as close as possible to dz_stretch_level_end.
-!--                stretch_factor_2 would guarantee that the stretched dz(n) is
-!--                equal to dz(n+1) after l_rounded grid levels.
+!--          stretch_factor_1 is taken to guarantee that the stretching
+!--          procedure ends as close as possible to dz_stretch_level_end.
+!--          stretch_factor_2 would guarantee that the stretched dz(n) is
+!--          equal to dz(n+1) after l_rounded grid levels.
              IF (delta_total_new < delta_total_old) THEN
                 dz_stretch_factor_array(n) = stretch_factor_1
@@ -1046,15 +1071,17 @@
             
           ENDDO
-             
+
+!
+!--    Second branch for stretching from fine to rough
        ELSEIF ( dz(n) < dz(n+1) ) THEN 
           DO WHILE ( stretch_factor_1 <= stretch_factor_upper_limit )
                      
-             stretch_factor_1 = 1.0 + iterations * stretch_factor_interval
-             distance = ABS( dz_stretch_level_end(n) -                      &
+             stretch_factor_1 = 1.0_wp + iterations * stretch_factor_interval
+             distance = ABS( dz_stretch_level_end(n) -                         &
                         dz_stretch_level_start(n) ) 
-             numerator = distance*stretch_factor_1/dz(n) +                  &
+             numerator = distance*stretch_factor_1/dz(n) +                     &
                          stretch_factor_1 - distance/dz(n)
              
-             l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0
+             l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0_wp
              l_rounded = NINT( l )
              delta_l = ABS( l_rounded - l ) / l
@@ -1062,6 +1089,6 @@
              stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) )
 
-             delta_stretch_factor = ABS( stretch_factor_1 -                 &
-                                        stretch_factor_2 ) /                &
+             delta_stretch_factor = ABS( stretch_factor_1 -                    &
+                                        stretch_factor_2 ) /                   &
                                         stretch_factor_2
              
@@ -1069,8 +1096,8 @@
              
 !
-!--                stretch_factor_1 is taken to guarantee that the stretching
-!--                procedure ends as close as possible to dz_stretch_level_end.
-!--                stretch_factor_2 would guarantee that the stretched dz(n) is
-!--                equal to dz(n+1) after l_rounded grid levels.
+!--          stretch_factor_1 is taken to guarantee that the stretching
+!--          procedure ends as close as possible to dz_stretch_level_end.
+!--          stretch_factor_2 would guarantee that the stretched dz(n) is
+!--          equal to dz(n+1) after l_rounded grid levels.
              IF (delta_total_new < delta_total_old) THEN
                 dz_stretch_factor_array(n) = stretch_factor_1
@@ -1794,33 +1821,33 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  bh            !< temporary vertical index of building height
-    INTEGER(iwp) ::  blx           !< grid point number of building size along x
-    INTEGER(iwp) ::  bly           !< grid point number of building size along y
-    INTEGER(iwp) ::  bxl           !< index for left building wall
-    INTEGER(iwp) ::  bxr           !< index for right building wall
-    INTEGER(iwp) ::  byn           !< index for north building wall
-    INTEGER(iwp) ::  bys           !< index for south building wall
-    INTEGER(iwp) ::  ch            !< temporary vertical index for canyon height
-    INTEGER(iwp) ::  cwx           !< grid point number of canyon size along x
-    INTEGER(iwp) ::  cwy           !< grid point number of canyon size along y
-    INTEGER(iwp) ::  cxl           !< index for left canyon wall
-    INTEGER(iwp) ::  cxr           !< index for right canyon wall
-    INTEGER(iwp) ::  cyn           !< index for north canyon wall
-    INTEGER(iwp) ::  cys           !< index for south canyon wall
-    INTEGER(iwp) ::  i             !< index variable along x
-    INTEGER(iwp) ::  j             !< index variable along y
-    INTEGER(iwp) ::  k             !< index variable along z
-    INTEGER(iwp) ::  hv_in         !< heavyside function to model inner tunnel surface 
-    INTEGER(iwp) ::  hv_out        !< heavyside function to model outer tunnel surface 
-    INTEGER(iwp) ::  txe_out       !< end position of outer tunnel wall in x
-    INTEGER(iwp) ::  txs_out       !< start position of outer tunnel wall in x
-    INTEGER(iwp) ::  tye_out       !< end position of outer tunnel wall in y
-    INTEGER(iwp) ::  tys_out       !< start position of outer tunnel wall in y
-    INTEGER(iwp) ::  txe_in        !< end position of inner tunnel wall in x
-    INTEGER(iwp) ::  txs_in        !< start position of inner tunnel wall in x
-    INTEGER(iwp) ::  tye_in        !< end position of inner tunnel wall in y
-    INTEGER(iwp) ::  tys_in        !< start position of inner tunnel wall in y
-    INTEGER(iwp) ::  td            !< tunnel wall depth
-    INTEGER(iwp) ::  th            !< height of outer tunnel wall
+    INTEGER(iwp) ::  bh                !< temporary vertical index of building height
+    INTEGER(iwp) ::  ngp_bx            !< grid point number of building size along x
+    INTEGER(iwp) ::  ngp_by            !< grid point number of building size along y
+    INTEGER(iwp) ::  index_left_bwall  !< index for left building wall
+    INTEGER(iwp) ::  index_right_bwall !< index for right building wall
+    INTEGER(iwp) ::  index_north_bwall !< index for north building wall
+    INTEGER(iwp) ::  index_south_bwall !< index for south building wall
+    INTEGER(iwp) ::  ch                !< temporary vertical index for canyon height
+    INTEGER(iwp) ::  ngp_cx            !< grid point number of canyon size along x
+    INTEGER(iwp) ::  ngp_cy            !< grid point number of canyon size along y
+    INTEGER(iwp) ::  index_left_cwall  !< index for left canyon wall
+    INTEGER(iwp) ::  index_right_cwall !< index for right canyon wall
+    INTEGER(iwp) ::  index_north_cwall !< index for north canyon wall
+    INTEGER(iwp) ::  index_south_cwall !< index for south canyon wall
+    INTEGER(iwp) ::  i                 !< index variable along x
+    INTEGER(iwp) ::  j                 !< index variable along y
+    INTEGER(iwp) ::  k                 !< index variable along z
+    INTEGER(iwp) ::  hv_in             !< heavyside function to model inner tunnel surface 
+    INTEGER(iwp) ::  hv_out            !< heavyside function to model outer tunnel surface 
+    INTEGER(iwp) ::  txe_out           !< end position of outer tunnel wall in x
+    INTEGER(iwp) ::  txs_out           !< start position of outer tunnel wall in x
+    INTEGER(iwp) ::  tye_out           !< end position of outer tunnel wall in y
+    INTEGER(iwp) ::  tys_out           !< start position of outer tunnel wall in y
+    INTEGER(iwp) ::  txe_in            !< end position of inner tunnel wall in x
+    INTEGER(iwp) ::  txs_in            !< start position of inner tunnel wall in x
+    INTEGER(iwp) ::  tye_in            !< end position of inner tunnel wall in y
+    INTEGER(iwp) ::  tys_in            !< start position of inner tunnel wall in y
+    INTEGER(iwp) ::  td                !< tunnel wall depth
+    INTEGER(iwp) ::  th                !< height of outer tunnel wall
 
     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_local         !< index for topography top at cell-center
@@ -1829,5 +1856,5 @@
 !-- Check for correct setting of the namelist parameter topography. If 
 !-- topography information is read from file but topography = 'flat', 
-!-- initialization does not properly.
+!-- initialization does not work properly.
     IF ( ( buildings_f%from_file  .OR.  terrain_height_f%from_file )  .AND.    &
            TRIM( topography ) /= 'read_from_file' )  THEN
@@ -1860,28 +1887,33 @@
 !--       Single rectangular building, by default centered in the middle of the
 !--       total domain
-          blx = NINT( building_length_x / dx )
-          bly = NINT( building_length_y / dy )
+          ngp_bx = NINT( building_length_x / dx )
+          ngp_by = NINT( building_length_y / dy )
           bh  = MINLOC( ABS( zw - building_height ), 1 ) - 1
           IF ( ABS( zw(bh)   - building_height ) == &
                ABS( zw(bh+1) - building_height )    )  bh = bh + 1
           IF ( building_wall_left == 9999999.9_wp )  THEN
-             building_wall_left = ( nx + 1 - blx ) / 2 * dx
-          ENDIF
-          bxl = NINT( building_wall_left / dx )
-          bxr = bxl + blx
+             building_wall_left = ( nx + 1 - ngp_bx ) / 2 * dx
+          ENDIF
+          index_left_bwall = NINT( building_wall_left / dx )
+          index_right_bwall = index_left_bwall + ngp_bx
 
           IF ( building_wall_south == 9999999.9_wp )  THEN
-              building_wall_south = ( ny + 1 - bly ) / 2 * dy
-          ENDIF
-          bys = NINT( building_wall_south / dy )
-          byn = bys + bly
+              building_wall_south = ( ny + 1 - ngp_by ) / 2 * dy
+          ENDIF
+          index_south_bwall = NINT( building_wall_south / dy )
+          index_north_bwall = index_south_bwall + ngp_by
 
 !
 !--       Building size has to meet some requirements
-          IF ( ( bxl < 1 ) .OR. ( bxr > nx-1 ) .OR. ( bxr < bxl+3 ) .OR.       &
-               ( bys < 1 ) .OR. ( byn > ny-1 ) .OR. ( byn < bys+3 ) )  THEN
+          IF ( ( index_left_bwall < 1 ) .OR. ( index_right_bwall > nx-1 ) .OR. &
+               ( index_right_bwall < index_left_bwall+3 ) .OR.                 &
+               ( index_south_bwall < 1 ) .OR. ( index_north_bwall > ny-1 ) .OR.&
+               ( index_north_bwall < index_south_bwall+3 ) )  THEN
              WRITE( message_string, * ) 'inconsistent building parameters:',   &
-                                      '&bxl=', bxl, 'bxr=', bxr, 'bys=', bys,  &
-                                      'byn=', byn, 'nx=', nx, 'ny=', ny
+                                      '&index_left_bwall=', index_left_bwall,  &
+                                      'index_right_bwall=', index_right_bwall, &
+                                      'index_south_bwall=', index_south_bwall, &
+                                      'index_north_bwall=', index_north_bwall, &
+                                      'nx=', nx, 'ny=', ny
              CALL message( 'init_grid', 'PA0203', 1, 2, 0, 6, 0 )
           ENDIF
@@ -1891,7 +1923,8 @@
 !
 !--       Define the building. 
-          IF ( bxl <= nxr  .AND.  bxr >= nxl  .AND.                            &
-               bys <= nyn  .AND.  byn >= nys )                                 & 
-             nzb_local(MAX(nys,bys):MIN(nyn,byn),MAX(nxl,bxl):MIN(nxr,bxr)) = bh
+          IF ( index_left_bwall <= nxr  .AND.  index_right_bwall >= nxl  .AND. &
+               index_south_bwall <= nyn  .AND.  index_north_bwall >= nys )     & 
+             nzb_local(MAX(nys,index_south_bwall):MIN(nyn,index_north_bwall),  &
+                       MAX(nxl,index_left_bwall):MIN(nxr,index_right_bwall)) = bh
 !
 !--       Set bit array on basis of nzb_local
@@ -1941,19 +1974,19 @@
 !
 !--          Street canyon in y direction
-             cwx = NINT( canyon_width_x / dx )
+             ngp_cx = NINT( canyon_width_x / dx )
              IF ( canyon_wall_left == 9999999.9_wp )  THEN
-                canyon_wall_left = ( nx + 1 - cwx ) / 2 * dx
-             ENDIF
-             cxl = NINT( canyon_wall_left / dx )
-             cxr = cxl + cwx
+                canyon_wall_left = ( nx + 1 - ngp_cx ) / 2 * dx
+             ENDIF
+             index_left_cwall= NINT( canyon_wall_left / dx )
+             index_right_cwall= index_left_cwall+ ngp_cx
           ELSEIF ( canyon_width_y /= 9999999.9_wp )  THEN
 !
 !--          Street canyon in x direction
-             cwy = NINT( canyon_width_y / dy )
+             ngp_cy = NINT( canyon_width_y / dy )
              IF ( canyon_wall_south == 9999999.9_wp )  THEN
-                canyon_wall_south = ( ny + 1 - cwy ) / 2 * dy
-             ENDIF
-             cys = NINT( canyon_wall_south / dy )
-             cyn = cys + cwy
+                canyon_wall_south = ( ny + 1 - ngp_cy ) / 2 * dy
+             ENDIF
+             index_south_cwall = NINT( canyon_wall_south / dy )
+             index_north_cwall = index_south_cwall + ngp_cy
       
           ELSE
@@ -1971,18 +2004,21 @@
 !--       Street canyon size has to meet some requirements
           IF ( canyon_width_x /= 9999999.9_wp )  THEN
-             IF ( ( cxl < 1 ) .OR. ( cxr > nx-1 ) .OR. ( cwx < 3 ) .OR.        &
+             IF ( ( index_left_cwall< 1 ) .OR. ( index_right_cwall> nx-1 ) .OR.&
+                  ( ngp_cx < 3 ) .OR. ( ch < 3 ) )  THEN
+                WRITE( message_string, * ) 'inconsistent canyon parameters:',  &
+                                           '&index_left_cwall=', index_left_cwall, &
+                                           ' index_right_cwall=', index_right_cwall, &
+                                           ' ngp_cx=', ngp_cx,                 &
+                                           ' ch=', ch, ' nx=', nx, ' ny=', ny
+                CALL message( 'init_grid', 'PA0205', 1, 2, 0, 6, 0 ) 
+             ENDIF
+          ELSEIF ( canyon_width_y /= 9999999.9_wp )  THEN
+             IF ( ( index_south_cwall < 1 ) .OR.                               &
+                  ( index_north_cwall > ny-1 ) .OR. ( ngp_cy < 3 ) .OR.        &
                   ( ch < 3 ) )  THEN
                 WRITE( message_string, * ) 'inconsistent canyon parameters:',  &
-                                           '&cxl=', cxl, ' cxr=', cxr,         &
-                                           ' cwx=', cwx,                       &
-                                           ' ch=', ch, ' nx=', nx, ' ny=', ny
-                CALL message( 'init_grid', 'PA0205', 1, 2, 0, 6, 0 ) 
-             ENDIF
-          ELSEIF ( canyon_width_y /= 9999999.9_wp )  THEN
-             IF ( ( cys < 1 ) .OR. ( cyn > ny-1 ) .OR. ( cwy < 3 ) .OR.        &
-                  ( ch < 3 ) )  THEN
-                WRITE( message_string, * ) 'inconsistent canyon parameters:',  &
-                                           '&cys=', cys, ' cyn=', cyn,         &
-                                           ' cwy=', cwy,                       &
+                                           '&index_south_cwall=', index_south_cwall, & 
+                                           ' index_north_cwall=', index_north_cwall, &
+                                           ' ngp_cy=', ngp_cy,                 &
                                            ' ch=', ch, ' nx=', nx, ' ny=', ny
                 CALL message( 'init_grid', 'PA0206', 1, 2, 0, 6, 0 ) 
@@ -2000,9 +2036,9 @@
           nzb_local = ch
           IF ( canyon_width_x /= 9999999.9_wp )  THEN
-             IF ( cxl <= nxr  .AND.  cxr >= nxl )                              &
-                nzb_local(:,MAX(nxl,cxl+1):MIN(nxr,cxr-1)) = 0
+             IF ( index_left_cwall<= nxr  .AND.  index_right_cwall>= nxl )     &
+                nzb_local(:,MAX(nxl,index_left_cwall+1):MIN(nxr,index_right_cwall-1)) = 0
           ELSEIF ( canyon_width_y /= 9999999.9_wp )  THEN
-             IF ( cys <= nyn  .AND.  cyn >= nys )                              &          
-                nzb_local(MAX(nys,cys+1):MIN(nyn,cyn-1),:) = 0
+             IF ( index_south_cwall <= nyn  .AND.  index_north_cwall >= nys )  &          
+                nzb_local(MAX(nys,index_south_cwall+1):MIN(nyn,index_north_cwall-1),:) = 0
           ENDIF
 !
