Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 3058)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 3065)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! dz was replaced by dz(1), error message revised
+! 
+! 3049 2018-05-29 13:52:36Z Giersch
 ! add variable description
 ! 
@@ -4178,5 +4181,5 @@
     IF ( ( constant_flux_layer .OR.  &
            INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )       &
-         .AND. roughness_length >= 0.5 * dz )  THEN
+         .AND. roughness_length >= 0.5 * dz(1) )  THEN
        message_string = 'roughness_length must be smaller than dz/2'
        CALL message( 'check_parameters', 'PA0424', 1, 2, 0, 6, 0 )
@@ -4198,5 +4201,6 @@
 !-- Check if vertical grid stretching is switched off in case of complex 
 !-- terrain simulations
-    IF ( complex_terrain  .AND.  dz_stretch_level < 100000.0_wp )  THEN
+    IF ( complex_terrain  .AND.                                                &
+         ANY( dz_stretch_level_start /= -9999999.9_wp ) )  THEN
        message_string = 'Vertical grid stretching is not allowed for ' //      &
                         'complex_terrain = .T.'
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 3058)
+++ palm/trunk/SOURCE/header.f90	(revision 3065)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Header output concerning stretching revised
+! 
+! 3045 2018-05-28 07:55:41Z Giersch
 ! Error messages revised
 ! 
@@ -800,17 +803,67 @@
 !-- Computational grid
     IF ( .NOT. ocean )  THEN
-       WRITE ( io, 250 )  dx, dy, dz, (nx+1)*dx, (ny+1)*dy, zu(nzt+1)
-       IF ( dz_stretch_level_index < nzt+1 )  THEN
-          WRITE ( io, 252 )  dz_stretch_level, dz_stretch_level_index, &
-                             dz_stretch_factor, dz_max
-       ENDIF
+       WRITE ( io, 250 )  dx, dy
+       
+       DO i = 1, number_stretch_level_start+1
+          WRITE ( io, 253 )  i, dz(i)
+       ENDDO
+       
+       WRITE( io, 251 ) (nx+1)*dx, (ny+1)*dy, zu(nzt+1)
+       
+       IF ( ANY( dz_stretch_level_start_index < nzt+1 ) )  THEN
+          WRITE( io, '(A)', advance='no') ' Vertical stretching starts at height:'
+          DO i = 1, number_stretch_level_start
+             WRITE ( io, '(F10.1,A3)', advance='no' )  dz_stretch_level_start(i), ' m,'
+          ENDDO
+          WRITE( io, '(/,A)', advance='no') ' Vertical stretching starts at index: '
+          DO i = 1, number_stretch_level_start
+             WRITE ( io, '(I12,A1)', advance='no' )  dz_stretch_level_start_index(i), ','
+          ENDDO
+          WRITE( io, '(/,A)', advance='no') ' Vertical stretching ends at height:  '
+          DO i = 1, number_stretch_level_start
+             WRITE ( io, '(F10.1,A3)', advance='no' )  dz_stretch_level_end(i), ' m,'
+          ENDDO
+          WRITE( io, '(/,A)', advance='no') ' Vertical stretching ends at index:   '
+          DO i = 1, number_stretch_level_start
+             WRITE ( io, '(I12,A1)', advance='no' )  dz_stretch_level_end_index(i), ','
+          ENDDO
+          WRITE( io, '(/,A)', advance='no') ' Factor used for stretching:          '
+          DO i = 1, number_stretch_level_start
+             WRITE ( io, '(F12.3,A1)', advance='no' )  dz_stretch_factor_array(i), ','
+          ENDDO
+       ENDIF
+       
     ELSE
-       WRITE ( io, 250 )  dx, dy, dz, (nx+1)*dx, (ny+1)*dy, zu(0)
-       IF ( dz_stretch_level_index > 0 )  THEN
-          WRITE ( io, 252 )  dz_stretch_level, dz_stretch_level_index, &
-                             dz_stretch_factor, dz_max
-       ENDIF
-    ENDIF
-    WRITE ( io, 254 )  nx, ny, nzt+1, MIN( nnx, nx+1 ), MIN( nny, ny+1 ), &
+       WRITE ( io, 250 )  dx, dy
+       DO i = 1, number_stretch_level_start+1
+          WRITE ( io, 253 )  i, dz(i)
+       ENDDO
+       
+       WRITE ( io, 251 ) (nx+1)*dx, (ny+1)*dy, zu(0)
+       
+       IF ( ANY( dz_stretch_level_start_index > 0 ) )  THEN
+          WRITE( io, '(A)', advance='no') ' Vertical stretching starts at height:'
+          DO i = 1, number_stretch_level_start
+             WRITE ( io, '(F10.1,A3)', advance='no' )  dz_stretch_level_start(i), ' m,'
+          ENDDO
+          WRITE( io, '(/,A)', advance='no') ' Vertical stretching starts at index: '
+          DO i = 1, number_stretch_level_start
+             WRITE ( io, '(I12,A1)', advance='no' )  dz_stretch_level_start_index(i), ','
+          ENDDO
+          WRITE( io, '(/,A)', advance='no') ' Vertical stretching ends at height:  '
+          DO i = 1, number_stretch_level_start
+             WRITE ( io, '(F10.1,A3)', advance='no' )  dz_stretch_level_end(i), ' m,'
+          ENDDO
+          WRITE( io, '(/,A)', advance='no') ' Vertical stretching ends at index:   '
+          DO i = 1, number_stretch_level_start
+             WRITE ( io, '(I12,A1)', advance='no' )  dz_stretch_level_end_index(i), ','
+          ENDDO
+          WRITE( io, '(/,A)', advance='no') ' Factor used for stretching:          '
+          DO i = 1, number_stretch_level_start
+             WRITE ( io, '(F12.3,A1)', advance='no' )  dz_stretch_factor_array(i), ','
+          ENDDO
+       ENDIF
+    ENDIF
+    WRITE ( io, 254 )  nx, ny, nzt+1, MIN( nnx, nx+1 ), MIN( nny, ny+1 ),      &
                        MIN( nnz+2, nzt+2 )
     IF ( sloping_surface )  WRITE ( io, 260 )  alpha_surface
@@ -2083,12 +2136,9 @@
 250 FORMAT (//' Computational grid and domain size:'/ &
               ' ----------------------------------'// &
-              ' Grid length:      dx =    ',F7.3,' m    dy =    ',F7.3, &
-              ' m    dz =    ',F7.3,' m'/ &
-              ' Domain size:       x = ',F10.3,' m     y = ',F10.3, &
+              ' Grid length:      dx =    ',F8.3,' m    dy =    ',F8.3, ' m')
+251 FORMAT (  /' Domain size:       x = ',F10.3,' m     y = ',F10.3, &
               ' m  z(u) = ',F10.3,' m'/)
-252 FORMAT (' dz constant up to ',F10.3,' m (k=',I4,'), then stretched by', &
-              ' factor:',F6.3/ &
-            ' maximum dz not to be exceeded is dz_max = ',F10.3,' m'/)
-254 FORMAT (' Number of gridpoints (x,y,z):  (0:',I4,', 0:',I4,', 0:',I4,')'/ &
+253 FORMAT ( '                dz(',I1,') =    ', F8.3, ' m')
+254 FORMAT (//' Number of gridpoints (x,y,z):  (0:',I4,', 0:',I4,', 0:',I4,')'/ &
             ' Subdomain size (x,y,z):        (  ',I4,',   ',I4,',   ',I4,')'/)
 260 FORMAT (/' The model has a slope in x-direction. Inclination angle: ',F6.2,&
Index: palm/trunk/SOURCE/init_grid.f90
===================================================================
--- palm/trunk/SOURCE/init_grid.f90	(revision 3058)
+++ palm/trunk/SOURCE/init_grid.f90	(revision 3065)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! New vertical stretching mechanism introduced
+! 
+! 3051 2018-05-30 17:43:55Z suehring
 ! Minor bugfix concerning mapping 3D buildings on top of terrain
 ! 
@@ -325,15 +328,19 @@
                canyon_height, canyon_wall_left, canyon_wall_south,             &
                canyon_width_x, canyon_width_y, constant_flux_layer,            &
-               dp_level_ind_b, dz, dz_max, dz_stretch_factor,                  &
-               dz_stretch_level, dz_stretch_level_index, grid_level,           &
+               dp_level_ind_b, dz, dz_max, dz_stretch_factor,                  &   
+               dz_stretch_factor_array, dz_stretch_level, dz_stretch_level_end,&
+               dz_stretch_level_end_index, dz_stretch_level_start_index,       &
+               dz_stretch_level_start, grid_level,                             &
                force_bound_l, force_bound_r, force_bound_n, force_bound_s,     &
                ibc_uv_b, inflow_l, inflow_n, inflow_r, inflow_s,               &
                masking_method, maximum_grid_level, message_string,             &
-               momentum_advec, nest_domain, nest_bound_l, nest_bound_n,        &
-               nest_bound_r, nest_bound_s, ocean, outflow_l, outflow_n,        &
-               outflow_r, outflow_s, psolver, scalar_advec, topography,        &
-               topography_grid_convention, tunnel_height, tunnel_length,       &
-               tunnel_width_x, tunnel_width_y, tunnel_wall_depth,              &
-               use_surface_fluxes, use_top_fluxes, wall_adjustment_factor
+               momentum_advec, nest_domain, nest_bound_l,                      &
+               nest_bound_n, nest_bound_r, nest_bound_s,                       &
+               number_stretch_level_end, number_stretch_level_start, ocean,    &
+               outflow_l, outflow_n, outflow_r, outflow_s, psolver,            & 
+               scalar_advec, topography, topography_grid_convention,           &
+               tunnel_height, tunnel_length, tunnel_width_x, tunnel_width_y,   &
+               tunnel_wall_depth, use_surface_fluxes, use_top_fluxes,          &
+               wall_adjustment_factor
          
     USE grid_variables,                                                        &
@@ -361,18 +368,22 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  i                !< index variable along x 
-    INTEGER(iwp) ::  j                !< index variable along y
-    INTEGER(iwp) ::  k                !< index variable along z
-    INTEGER(iwp) ::  k_top            !< topography top index on local PE
-    INTEGER(iwp) ::  l                !< loop variable 
-    INTEGER(iwp) ::  nzb_local_max    !< vertical grid index of maximum topography height
-    INTEGER(iwp) ::  nzb_local_min    !< vertical grid index of minimum topography height
+    INTEGER(iwp) ::  i                           !< index variable along x 
+    INTEGER(iwp) ::  j                           !< index variable along y
+    INTEGER(iwp) ::  k                           !< index variable along z
+    INTEGER(iwp) ::  k_top                       !< topography top index on local PE
+    INTEGER(iwp) ::  n                           !< loop variable for stretching
+    INTEGER(iwp) ::  number_dz                   !< number of user-specified dz values       
+    INTEGER(iwp) ::  nzb_local_max               !< vertical grid index of maximum topography height
+    INTEGER(iwp) ::  nzb_local_min               !< vertical grid index of minimum topography height
                                       
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_local      !< index for topography top at cell-center
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_tmp        !< dummy to calculate topography indices on u- and v-grid
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_local  !< index for topography top at cell-center
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_tmp    !< dummy to calculate topography indices on u- and v-grid
 
     INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  topo !< input array for 3D topography and dummy array for setting "outer"-flags
 
+    REAL(wp) ::  dz_level_end  !< distance between calculated height level for u/v-grid and user-specified end level for stretching
     REAL(wp) ::  dz_stretched  !< stretched vertical grid spacing
+    
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  min_dz_stretch_level_end !< Array that contains all minimum heights where the stretching can end
 
 
@@ -391,37 +402,198 @@
 !
 !-- Compute height of u-levels from constant grid length and dz stretch factors
-    IF ( dz == -1.0_wp )  THEN
+    IF ( dz(1) == -1.0_wp )  THEN
        message_string = 'missing dz'
        CALL message( 'init_grid', 'PA0200', 1, 2, 0, 6, 0 ) 
-    ELSEIF ( dz <= 0.0_wp )  THEN
-       WRITE( message_string, * ) 'dz=',dz,' <= 0.0'
+    ELSEIF ( dz(1) <= 0.0_wp )  THEN
+       WRITE( message_string, * ) 'dz=',dz(1),' <= 0.0'
        CALL message( 'init_grid', 'PA0201', 1, 2, 0, 6, 0 )
     ENDIF
 
 !
+!-- Initialize dz_stretch_level_start with the value of dz_stretch_level
+!-- if it was set by the user
+    IF ( dz_stretch_level /= -9999999.9_wp ) THEN
+       dz_stretch_level_start(1) = dz_stretch_level
+    ENDIF
+       
+!
+!-- Determine number of dz values and stretching levels specified by the
+!-- user to allow right controlling of the stretching mechanism and to
+!-- perform error checks
+    number_dz = COUNT( dz /= -1.0_wp )
+    number_stretch_level_start = COUNT( dz_stretch_level_start /=              &
+                                       -9999999.9_wp )
+    number_stretch_level_end = COUNT( dz_stretch_level_end /=                  &
+                                      9999999.9_wp )
+
+!
+!-- The number of specified end levels +1 has to be the same than 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 = ',           &
+                                   number_stretch_level_end+1
+          CALL message( 'init_grid', 'PA0156', 1, 2, 0, 6, 0 )
+    ENDIF
+    
+!
+!--    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 = ',             &
+                                   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.   &
+         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 = ',    &
+                                   number_stretch_level_end
+          CALL message( 'init_grid', 'PA0216', 1, 2, 0, 6, 0 )
+    ENDIF
+
+!
+!-- Initialize dz for the free atmosphere with the value of dz_max
+    IF ( dz(number_stretch_level_start+1) == -1.0_wp .AND.                     &
+         number_stretch_level_start /= 0 ) THEN 
+       dz(number_stretch_level_start+1) = dz_max
+    ENDIF
+       
+!
+!-- Initialize the stretching factor if (infinitely) stretching in the free 
+!-- atmosphere is desired (dz_stretch_level_end was not specified for the 
+!-- free atmosphere)
+    IF ( number_stretch_level_start == number_stretch_level_end + 1 ) THEN 
+       dz_stretch_factor_array(number_stretch_level_start) =                   &
+       dz_stretch_factor
+    ENDIF
+    
+!
+!-- Allocation of arrays for stretching
+    ALLOCATE( min_dz_stretch_level_end(number_stretch_level_start) )
+    
+!
 !-- Define the vertical grid levels
     IF ( .NOT. ocean )  THEN
+    
+!
+!--    The stretching region has to be large enough to allow for a smooth
+!--    transition between two different grid spacings
+       DO n = 1, number_stretch_level_start
+          min_dz_stretch_level_end(n) = dz_stretch_level_start(n) +            &
+                                        4 * MAX( dz(n),dz(n+1) )
+       ENDDO
+
+       IF ( ANY( min_dz_stretch_level_end > dz_stretch_level_end ) ) THEN
+             message_string= 'Eeach dz_stretch_level_end has to be larger ' // &
+                             'than its corresponding value for &' //           &
+                             'dz_stretch_level_start + 4*MAX(dz(n),dz(n+1)) '//&
+                             'to allow for smooth grid stretching'
+             CALL message( 'init_grid', 'PA0224', 1, 2, 0, 6, 0 )
+       ENDIF
+       
+!
+!--    Stretching must not be applied within the prandtl_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 ',&
+                                     'larger than ', dz(1) * 1.5
+             CALL message( 'init_grid', 'PA0226', 1, 2, 0, 6, 0 )
+       ENDIF
+
+!
+!--    The stretching has to start and end on a grid level. Therefore 
+!--    user-specified values have to ''interpolate'' to the next lowest level
+       IF ( number_stretch_level_start /= 0 ) THEN
+          dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) -        &
+                                            dz(1)/2.0) / dz(1) )               &
+                                      * dz(1) + dz(1)/2.0
+       ENDIF
+       
+       IF ( number_stretch_level_start > 1 ) THEN
+          DO n = 2, number_stretch_level_start
+             dz_stretch_level_start(n) = INT( dz_stretch_level_start(n) /      &
+                                              dz(n) ) * dz(n)
+          ENDDO
+       ENDIF
+       
+       IF ( number_stretch_level_end /= 0 ) THEN
+          DO n = 1, number_stretch_level_end
+             dz_stretch_level_end(n) = INT( dz_stretch_level_end(n) /          &
+                                            dz(n+1) ) * dz(n+1)
+          ENDDO
+       ENDIF
+ 
+!
+!--    Determine stretching factor if necessary
+       IF ( number_stretch_level_end >= 1 ) THEN 
+          CALL calculate_stretching_factor( number_stretch_level_end )
+       ENDIF
+
 !
 !--    Grid for atmosphere with surface at z=0 (k=0, w-grid).
+!--    First compute the u- and v-levels. In case of dirichlet bc for u and v
+!--    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.
-
        IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 ) THEN
           zu(0) = 0.0_wp
-      !    zu(0) = - dz * 0.5_wp
        ELSE
-          zu(0) = - dz * 0.5_wp
-       ENDIF
-       zu(1) =   dz * 0.5_wp
-
-       dz_stretch_level_index = nzt+1
-       dz_stretched = dz
+          zu(0) = - dz(1) * 0.5_wp
+       ENDIF
+          
+       zu(1) =   dz(1) * 0.5_wp
+       
+!
+!--    Determine u and v height levels considering the possibility of grid
+!--    stretching in several heights.
+       n = 1
+       dz_stretch_level_start_index = nzt+1
+       dz_stretch_level_end_index = nzt+1
+       dz_stretched = dz(1)
+
+!--    The default value of dz_stretch_level_start is negative, thus the first
+!--    condition is always true. Hence, the second condition is necessary.
        DO  k = 2, nzt+1
-          IF ( dz_stretch_level <= zu(k-1)  .AND.  dz_stretched < dz_max )  THEN
-             dz_stretched = dz_stretched * dz_stretch_factor
-             dz_stretched = MIN( dz_stretched, dz_max )
-             IF ( dz_stretch_level_index == nzt+1 ) dz_stretch_level_index = k-1
-          ENDIF
+          IF ( dz_stretch_level_start(n) <= zu(k-1) .AND.                      &
+               dz_stretch_level_start(n) /= -9999999.9_wp ) THEN
+             dz_stretched = dz_stretched * dz_stretch_factor_array(n)
+             
+             IF ( dz(n) > dz(n+1) ) THEN
+                dz_stretched = MAX( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (higher) dz
+             ELSE
+                dz_stretched = MIN( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (lower) dz
+             ENDIF
+             
+             IF ( dz_stretch_level_start_index(n) == nzt+1 )                         &
+             dz_stretch_level_start_index(n) = k-1
+             
+          ENDIF
+          
           zu(k) = zu(k-1) + dz_stretched
+          
+!
+!--       Make sure that the stretching ends exactly at dz_stretch_level_end 
+          dz_level_end = ABS( zu(k) - dz_stretch_level_end(n) ) 
+          
+          IF ( dz_level_end  < dz(n+1)/3.0 ) THEN
+             zu(k) = dz_stretch_level_end(n)
+             dz_stretched = dz(n+1)
+             dz_stretch_level_end_index(n) = k
+             n = n + 1             
+          ENDIF
        ENDDO
 
@@ -438,4 +610,60 @@
 
     ELSE
+
+!
+!--    The stretching region has to be large enough to allow for a smooth
+!--    transition between two different grid spacings
+       DO n = 1, number_stretch_level_start
+          min_dz_stretch_level_end(n) = dz_stretch_level_start(n) -            &
+                                        4 * MAX( dz(n),dz(n+1) )
+       ENDDO
+       
+       IF ( ANY( min_dz_stretch_level_end < dz_stretch_level_end ) ) THEN
+             message_string= 'Eeach dz_stretch_level_end has to be less ' //   &
+                             'than its corresponding value for &' //           &
+                             'dz_stretch_level_start - 4*MAX(dz(n),dz(n+1)) '//&
+                             'to allow for smooth grid stretching'
+             CALL message( 'init_grid', 'PA0224', 1, 2, 0, 6, 0 )
+       ENDIF
+       
+!
+!--    Stretching must not be applied within the prandtl_layer 
+!--    (first two grid points). For the default case dz_stretch_level_start 
+!--    is negative. Therefore the absolut value is checked here.
+       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
+             CALL message( 'init_grid', 'PA0226', 1, 2, 0, 6, 0 )
+       ENDIF
+
+!
+!--    The stretching has to start and end on a grid level. Therefore 
+!--    user-specified values have to ''interpolate'' to the next highest level
+       IF ( number_stretch_level_start /= 0 ) THEN
+          dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) +        &
+                                            dz(1)/2.0) / dz(1) )               &
+                                      * dz(1) - dz(1)/2.0
+       ENDIF
+       
+       IF ( number_stretch_level_start > 1 ) THEN
+          DO n = 2, number_stretch_level_start
+             dz_stretch_level_start(n) = INT( dz_stretch_level_start(n) /      &
+                                              dz(n) ) * dz(n)
+          ENDDO
+       ENDIF
+       
+       IF ( number_stretch_level_end /= 0 ) THEN
+          DO n = 1, number_stretch_level_end
+             dz_stretch_level_end(n) = INT( dz_stretch_level_end(n) /          &
+                                            dz(n+1) ) * dz(n+1)
+          ENDDO
+       ENDIF
+       
+!
+!--    Determine stretching factor if necessary
+       IF ( number_stretch_level_end >= 1 ) THEN 
+          CALL calculate_stretching_factor( number_stretch_level_end )
+       ENDIF
+
 !
 !--    Grid for ocean with free water surface is at k=nzt (w-grid).
@@ -444,22 +672,46 @@
 !--    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.
-       zu(nzt+1) =   dz * 0.5_wp
-       zu(nzt)   = - dz * 0.5_wp
-
-       dz_stretch_level_index = 0
-       dz_stretched = dz
+!--    z values are negative starting from z=0 (surface)
+       zu(nzt+1) =   dz(1) * 0.5_wp
+       zu(nzt)   = - dz(1) * 0.5_wp
+
+!
+!--    Determine u and v height levels considering the possibility of grid
+!--    stretching in several heights.
+       n = 1
+       dz_stretch_level_start_index = 0
+       dz_stretch_level_end_index = 0
+       dz_stretched = dz(1)
+
        DO  k = nzt-1, 0, -1
-!
-!--       The default value of dz_stretch_level is positive, thus the first 
-!--       condition is always true. Hence, the second condition is necessary.
-          IF ( dz_stretch_level >= zu(k+1)  .AND.  dz_stretch_level <= 0.0  &
-               .AND.  dz_stretched < dz_max )  THEN
-             dz_stretched = dz_stretched * dz_stretch_factor
-             dz_stretched = MIN( dz_stretched, dz_max )
-             IF ( dz_stretch_level_index == 0 ) dz_stretch_level_index = k+1
-          ENDIF
+          
+          IF ( dz_stretch_level_start(n) >= zu(k+1) ) THEN
+             dz_stretched = dz_stretched * dz_stretch_factor_array(n)
+
+             IF ( dz(n) > dz(n+1) ) THEN
+                dz_stretched = MAX( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (higher) dz
+             ELSE
+                dz_stretched = MIN( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (lower) dz
+             ENDIF
+             
+             IF ( dz_stretch_level_start_index(n) == 0 )                             &
+             dz_stretch_level_start_index(n) = k+1
+             
+          ENDIF
+          
           zu(k) = zu(k+1) - dz_stretched
+          
+!
+!--       Make sure that the stretching ends exactly at dz_stretch_level_end 
+          dz_level_end = ABS( zu(k) - dz_stretch_level_end(n) ) 
+          
+          IF ( dz_level_end  < dz(n+1)/3.0 ) THEN
+             zu(k) = dz_stretch_level_end(n)
+             dz_stretched = dz(n+1)
+             dz_stretch_level_end_index(n) = k
+             n = n + 1             
+          ENDIF
        ENDDO
-
+       
 !
 !--    Compute the w-levels. They are always staggered half-way between the 
@@ -468,5 +720,5 @@
 !--    same height. The top w-level (nzt+1) is not used but set for 
 !--    consistency, since w and all scalar variables are defined up tp nzt+1.
-       zw(nzt+1) = dz
+       zw(nzt+1) = dz(1)
        zw(nzt)   = 0.0_wp
        DO  k = 0, nzt
@@ -855,4 +1107,134 @@
  END SUBROUTINE init_grid
 
+
+! Description:
+! -----------------------------------------------------------------------------!
+!> Calculation of the stretching factor through an iterative method. Ideas were 
+!> taken from the paper "Regional stretched grid generation and its application
+!> to the NCAR RegCM (1999)". Normally, no analytic solution exists because the
+!> system of equations has two variables (r,l) but four requirements 
+!> (l=integer, r=[0,88;1,2], Eq(6), Eq(5) starting from index j=1) which
+!> results into an overdetermined system. 
+!------------------------------------------------------------------------------!
+ SUBROUTINE calculate_stretching_factor( number_end )
+ 
+    USE control_parameters,                                                    &
+        ONLY:  dz, dz_stretch_factor, dz_stretch_factor_array,                 &   
+               dz_stretch_level_end, dz_stretch_level_start, message_string
+  
+    USE kinds
+    
+    IMPLICIT NONE
+    
+    INTEGER(iwp) ::  iterations  !< number of iterations until stretch_factor_lower/upper_limit is reached  
+    INTEGER(iwp) ::  l_rounded   !< after l_rounded grid levels dz(n) is strechted to dz(n+1) with stretch_factor_2 
+    INTEGER(iwp) ::  n           !< loop variable for stretching
+    
+    INTEGER(iwp), INTENT(IN) ::  number_end !< number of user-specified end levels for stretching
+        
+    REAL(wp) ::  delta_l               !< absolute difference between l and l_rounded
+    REAL(wp) ::  delta_stretch_factor  !< absolute difference between stretch_factor_1 and stretch_factor_2
+    REAL(wp) ::  delta_total_new       !< sum of delta_l and delta_stretch_factor for the next iteration (should be as small as possible) 
+    REAL(wp) ::  delta_total_old       !< sum of delta_l and delta_stretch_factor for the last iteration 
+    REAL(wp) ::  distance              !< distance between dz_stretch_level_start and dz_stretch_level_end (stretching region)
+    REAL(wp) ::  l                     !< value that fulfil Eq. (5) in the paper mentioned above together with stretch_factor_1 exactly
+    REAL(wp) ::  numerator             !< numerator of the quotient
+    REAL(wp) ::  stretch_factor_1      !< stretching factor that fulfil Eq. (5) togehter with l exactly
+    REAL(wp) ::  stretch_factor_2      !< stretching factor that fulfil Eq. (6) togehter with l_rounded exactly
+    
+    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
+ 
+ 
+       l = 0
+       DO  n = 1, number_end
+       
+          iterations = 1
+          stretch_factor_1 = 1.0 
+          stretch_factor_2 = 1.0
+          delta_total_old = 1.0
+          
+          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 - distance/dz(n)
+                
+                IF ( numerator > 0.0 ) THEN
+                   l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0
+                   l_rounded = NINT( l )
+                   delta_l = ABS( l_rounded - l ) / l
+                ENDIF
+                
+                stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) )
+                
+                delta_stretch_factor = ABS( stretch_factor_1 -              &
+                                            stretch_factor_2 ) /            &
+                                       stretch_factor_2
+                
+                delta_total_new = delta_l + delta_stretch_factor
+
+!
+!--                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
+                   delta_total_old = delta_total_new
+                ENDIF
+                
+                iterations = iterations + 1
+               
+             ENDDO
+                
+          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) -                   &
+                           dz_stretch_level_start(n) ) 
+                numerator = distance*stretch_factor_1/dz(n) +               &
+                            stretch_factor_1 - distance/dz(n)
+                
+                l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0
+                l_rounded = NINT( l )
+                delta_l = ABS( l_rounded - l ) / l
+                
+                stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) )
+
+                delta_stretch_factor = ABS( stretch_factor_1 -              &
+                                           stretch_factor_2 ) /             &
+                                           stretch_factor_2
+                
+                delta_total_new = delta_l + delta_stretch_factor
+                
+!
+!--                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
+                   delta_total_old = delta_total_new
+                ENDIF
+                
+                iterations = iterations + 1
+             ENDDO
+             
+          ELSE
+             message_string= 'Two adjacent values of dz must be different'
+             CALL message( 'init_grid', 'PA0228', 1, 2, 0, 6, 0 )
+             
+          ENDIF
+       ENDDO
+        
+ END SUBROUTINE calculate_stretching_factor
+ 
+ 
 ! Description:
 ! -----------------------------------------------------------------------------!
@@ -1703,5 +2085,5 @@
 !--       Tunnel-wall depth
           IF ( tunnel_wall_depth == 9999999.9_wp )  THEN  
-             td = MAX ( dx, dy, dz )
+             td = MAX ( dx, dy, dz(1) )
           ELSE
              td = tunnel_wall_depth
Index: palm/trunk/SOURCE/init_masks.f90
===================================================================
--- palm/trunk/SOURCE/init_masks.f90	(revision 3058)
+++ palm/trunk/SOURCE/init_masks.f90	(revision 3065)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! dz_stretch_level was replaced by dz_stretch_level_start
+! 
+! 3049 2018-05-29 13:52:36Z Giersch
 ! Error messages revised
 ! 
@@ -133,6 +136,6 @@
         ONLY:  constant_diffusion, cloud_droplets, cloud_physics,              &
                data_output_masks, data_output_masks_user,                      &
-               doav, doav_n, domask, domask_no, dz, dz_stretch_level, humidity,&
-               mask, masks, mask_scale, mask_i,                                &
+               doav, doav_n, domask, domask_no, dz, dz_stretch_level_start,    &
+               humidity, mask, masks, mask_scale, mask_i,                      &
                mask_i_global, mask_j, mask_j_global, mask_k, mask_k_global,    &
                mask_loop, mask_size, mask_size_l, mask_start_l, mask_x,        &
@@ -141,4 +144,5 @@
                microphysics_morrison, microphysics_seifert, passive_scalar,    &
                ocean, varnamelength
+               
 
     USE grid_variables,                                                        &
@@ -490,5 +494,5 @@
        CALL set_mask_locations( 1, dx, 'dx', nx, 'nx', nxl, nxr )
        CALL set_mask_locations( 2, dy, 'dy', ny, 'ny', nys, nyn )
-       CALL set_mask_locations( 3, dz, 'dz', nz, 'nz', nzb, nzt )
+       CALL set_mask_locations( 3, dz(1), 'dz', nz, 'nz', nzb, nzt )
 !
 !--    Set global masks along all three dimensions (required by
@@ -727,12 +731,14 @@
 !--          The following line assumes a constant vertical grid spacing within
 !--          the vertical mask range; it fails for vertical grid stretching.
-!--          Maybe revise later. Issue warning but continue execution.
+!--          Maybe revise later. Issue warning but continue execution. ABS(...)
+!--          within the IF statement is necessary because the default value of
+!--          dz_stretch_level_start is -9999999.9_wp.
              loop_stride = NINT( mask_loop(mid,dim,3) * mask_scale(dim) * ddxyz )
 
-             IF ( mask_loop(mid,dim,2) * mask_scale(dim) > dz_stretch_level )  &
-                  THEN
+             IF ( mask_loop(mid,dim,2) * mask_scale(dim) >                     &
+                  ABS( dz_stretch_level_start(1) ) )  THEN
                 WRITE ( message_string, '(A,I3,A,I1,A,F9.3,A,F8.2,3A)' )       &
                      'mask_loop(',mid,',',dim,',2)=', mask_loop(mid,dim,2),    &
-                     ' exceeds dz_stretch_level=',dz_stretch_level,            &
+                     ' exceeds dz_stretch_level=',dz_stretch_level_start(1),   &
                      '.&Vertical mask locations will not ',                    &
                      'match the desired heights within the stretching ',       &
Index: palm/trunk/SOURCE/lpm_advec.f90
===================================================================
--- palm/trunk/SOURCE/lpm_advec.f90	(revision 3058)
+++ palm/trunk/SOURCE/lpm_advec.f90	(revision 3065)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! dz values were replaced by dzw or dz(1) to allow for right vertical stretching
+! 
+! 2969 2018-04-13 11:55:09Z thiele
 ! Bugfix in Interpolation indices.
 ! 
@@ -158,5 +161,5 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  de_dx, de_dy, de_dz, diss, e, km, u, v, w, zu, zw
+        ONLY:  de_dx, de_dy, de_dz, diss, dzw, e, km, u, v, w, zu, zw
 
     USE cpulog
@@ -409,5 +412,5 @@
                             + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) *           &
                             u(k+1,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
-                u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dz *                  &
+                u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dzw(k) *               &
                            ( u_int_u - u_int_l )
              ENDIF
@@ -503,5 +506,5 @@
                           + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
                           ) / ( 3.0_wp * gg ) - v_gtrans
-                v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dz * &
+                v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dzw(k) *               &
                                   ( v_int_u - v_int_l )
              ENDIF
@@ -540,5 +543,5 @@
                             ( gg-dd ) * w(k+1,j+1,i+1) &
                           ) / ( 3.0_wp * gg )
-                w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dz * &
+                w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dzw(k) *               &
                            ( w_int_u - w_int_l )
              ENDIF
@@ -611,5 +614,5 @@
                                ( gg - dd ) * e(k+1,j+1,i+1) &
                             ) / ( 3.0_wp * gg )
-                   e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
+                   e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dzw(k) *            &
                                      ( e_int_u - e_int_l )
                 ENDIF
@@ -637,5 +640,5 @@
                                    ( gg - dd ) * de_dx(k+1,j+1,i+1) &
                                   ) / ( 3.0_wp * gg )
-                   de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / dz * &
+                   de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / dzw(k) *    &
                                               ( de_dx_int_u - de_dx_int_l )
                 ENDIF
@@ -655,5 +658,5 @@
                                    ( gg - dd ) * de_dy(k+1,j+1,i+1) &
                                   ) / ( 3.0_wp * gg )
-                      de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / dz * &
+                      de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / dzw(k) * &
                                                  ( de_dy_int_u - de_dy_int_l )
                 ENDIF
@@ -661,5 +664,5 @@
 !
 !--             Interpolate the TKE gradient along z
-                IF ( zv(n) < 0.5_wp * dz )  THEN
+                IF ( zv(n) < 0.5_wp * dz(1) )  THEN
                    de_dz_int(n) = 0.0_wp
                 ELSE
@@ -678,5 +681,5 @@
                                       ( gg - dd ) * de_dz(k+1,j+1,i+1) &
                                      ) / ( 3.0_wp * gg )
-                      de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) / dz * &
+                      de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) / dzw(k) * &
                                                  ( de_dz_int_u - de_dz_int_l )
                    ENDIF
@@ -699,5 +702,5 @@
                                   ( gg - dd ) * diss(k+1,j+1,i+1) &
                                  ) / ( 3.0_wp * gg )
-                   diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz * &
+                   diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dzw(k) *      &
                                             ( diss_int_u - diss_int_l )
                 ENDIF
Index: palm/trunk/SOURCE/lpm_exchange_horiz.f90
===================================================================
--- palm/trunk/SOURCE/lpm_exchange_horiz.f90	(revision 3058)
+++ palm/trunk/SOURCE/lpm_exchange_horiz.f90	(revision 3065)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! dz was replaced by dz(1) to allow for right vertical stretching
+! 
+! 3049 2018-05-29 13:52:36Z Giersch
 ! Error messages revised
 ! 
@@ -910,5 +913,5 @@
        ip = particle_array(n)%x * ddx
        jp = particle_array(n)%y * ddy
-       kp = particle_array(n)%z / dz + 1 + offset_ocean_nzt
+       kp = particle_array(n)%z / dz(1) + 1 + offset_ocean_nzt
 !
 !--    In case of grid stretching a particle might be above or below the  
Index: palm/trunk/SOURCE/lpm_init.f90
===================================================================
--- palm/trunk/SOURCE/lpm_init.f90	(revision 3058)
+++ palm/trunk/SOURCE/lpm_init.f90	(revision 3065)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! dz was replaced by dzw or dz(1) to allow for right vertical stretching
+! 
+! 3049 2018-05-29 13:52:36Z Giersch
 ! Error messages revised
 ! 
@@ -794,5 +797,5 @@
                                ip = tmp_particle%x * ddx
                                jp = tmp_particle%y * ddy
-                               kp = tmp_particle%z / dz + 1 + offset_ocean_nzt                               
+                               kp = tmp_particle%z / dz(1) + 1 + offset_ocean_nzt                               
                                DO WHILE( zw(kp) < tmp_particle%z ) 
                                   kp = kp + 1
@@ -969,6 +972,6 @@
                                      pdz(particles(n)%group)
                       particles(n)%z = particles(n)%z +                        &
-                              MERGE( rand_contr, SIGN( dz, rand_contr ),       &
-                                     ABS( rand_contr ) < dz                    &
+                              MERGE( rand_contr, SIGN( dzw(kp), rand_contr ),  &
+                                     ABS( rand_contr ) < dzw(kp)               &
                                    )
                    ENDIF
@@ -987,5 +990,5 @@
                    i = particles(n)%x * ddx
                    j = particles(n)%y * ddy
-                   k = particles(n)%z / dz + 1 + offset_ocean_nzt
+                   k = particles(n)%z / dz(1) + 1 + offset_ocean_nzt
                    DO WHILE( zw(k) < particles(n)%z )
                       k = k + 1
Index: palm/trunk/SOURCE/lpm_set_attributes.f90
===================================================================
--- palm/trunk/SOURCE/lpm_set_attributes.f90	(revision 3058)
+++ palm/trunk/SOURCE/lpm_set_attributes.f90	(revision 3065)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! dz was replaced by dzw to allow for right vertical stretching
+! 
+! 2718 2018-01-02 08:49:38Z maronga
 ! Corrected "Former revisions" section
 ! 
@@ -80,8 +83,8 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  pt, u, v, w, zu, zw
+        ONLY:  dzw, pt, u, v, w, zu, zw
 
     USE control_parameters,                                                    &
-        ONLY:  u_gtrans, v_gtrans, dz
+        ONLY:  u_gtrans, v_gtrans
 
     USE cpulog,                                                                &
@@ -207,5 +210,5 @@
                                      ( gg - dd ) * u(k+1,j+1,i+1)                &
                                    ) / ( 3.0_wp * gg ) - u_gtrans
-                         u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dz *         &
+                         u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dzw(k) *      &
                                            ( u_int_u - u_int_l )
                       ENDIF
@@ -240,5 +243,5 @@
                                       ( gg - dd ) * v(k+1,j+1,i+1)                &
                                     ) / ( 3.0_wp * gg ) - v_gtrans
-                         v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dz *         &
+                         v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dzw(k) *      &
                                            ( v_int_u - v_int_l )
                       ENDIF
@@ -344,5 +347,5 @@
                                  ) / ( 3.0_wp * gg ) - sums(k,4)
 
-                      pt_int = pt_int_l + ( zv(n) - zu(k) ) / dz *    &
+                      pt_int = pt_int_l + ( zv(n) - zu(k) ) / dzw(k) *          &
                                           ( pt_int_u - pt_int_l )
 
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 3058)
+++ palm/trunk/SOURCE/modules.f90	(revision 3065)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Variables concerning stretching introduced or revised
+! 
+! 3045 2018-05-28 07:55:41Z Giersch
 ! z_max_do2d removed
 ! 
@@ -1148,5 +1151,4 @@
     INTEGER(iwp) ::  dp_level_ind_b = 0                !< lowest grid index for external pressure gradient forcing
     INTEGER(iwp) ::  dvrp_filecount = 0                !< parameter for dvr visualization software
-    INTEGER(iwp) ::  dz_stretch_level_index            !< vertical grid level index above which the vertical grid spacing is stretched
     INTEGER(iwp) ::  ensemble_member_nr = 0            !< namelist parameter
     INTEGER(iwp) ::  gamma_mg                          !< switch for steering the multigrid cycle: 1: v-cycle, 2: w-cycle
@@ -1190,4 +1192,6 @@
     INTEGER(iwp) ::  num_var_fl                        !< number of sampling/output variables in virtual flight measurements
     INTEGER(iwp) ::  num_var_fl_user=0                 !< number of user-defined sampling/output variables in virtual flight measurements
+    INTEGER(iwp) ::  number_stretch_level_start        !< number of user-specified start levels for stretching
+    INTEGER(iwp) ::  number_stretch_level_end          !< number of user-specified end levels for stretching
     INTEGER(iwp) ::  nz_do3d = -9999                   !< namelist parameter
     INTEGER(iwp) ::  prt_time_count = 0                !< number of output intervals for particle data output
@@ -1199,4 +1203,5 @@
     INTEGER(iwp) ::  timestep_count = 0                !< number of timesteps carried out since the beginning of the initial run
     INTEGER(iwp) ::  y_shift = 0                       !< namelist parameter
+    
     INTEGER(iwp) ::  dist_nxl(0:1)                               !< left boundary of disturbance region
     INTEGER(iwp) ::  dist_nxr(0:1)                               !< right boundary of disturbance region
@@ -1211,4 +1216,6 @@
     INTEGER(iwp) ::  domask_no(max_masks,0:1) = 0                !< number of masked output quantities
     INTEGER(iwp) ::  domask_time_count(max_masks,0:1)            !< number of output intervals for masked data
+    INTEGER(iwp) ::  dz_stretch_level_end_index(9)               !< vertical grid level index until which the vertical grid spacing is stretched
+    INTEGER(iwp) ::  dz_stretch_level_start_index(9)             !< vertical grid level index above which the vertical grid spacing is stretched
     INTEGER(iwp) ::  mask_size(max_masks,3) = -1                 !< size of mask array per mask and dimension (for netcdf output)
     INTEGER(iwp) ::  mask_size_l(max_masks,3) = -1               !< subdomain size of mask array per mask and dimension (for netcdf output)
@@ -1416,8 +1423,7 @@
     REAL(wp) ::  dt_spinup = 60.0_wp                           !< namelist parameter
     REAL(wp) ::  dt_3d = 1.0_wp                                !< time step
-    REAL(wp) ::  dz = -1.0_wp                                  !< namelist parameter
-    REAL(wp) ::  dz_max = 9999999.9_wp                         !< namelist parameter
+    REAL(wp) ::  dz_max = 1000.0_wp                            !< namelist parameter
     REAL(wp) ::  dz_stretch_factor = 1.08_wp                   !< namelist parameter
-    REAL(wp) ::  dz_stretch_level = 100000.0_wp                !< namelist parameter
+    REAL(wp) ::  dz_stretch_level = -9999999.9_wp              !< namelist parameter
     REAL(wp) ::  e_init = 0.0_wp                               !< namelist parameter
     REAL(wp) ::  e_min = 0.0_wp                                !< namelist parameter
@@ -1529,4 +1535,8 @@
     REAL(wp) ::  dpdxy(1:2) = 0.0_wp                               !< namelist parameter
     REAL(wp) ::  dt_domask(max_masks) = 9999999.9_wp               !< namelist parameter
+    REAL(wp) ::  dz(10) = -1.0_wp                                  !< namelist parameter
+    REAL(wp) ::  dz_stretch_level_start(9) = -9999999.9_wp         !< namelist parameter
+    REAL(wp) ::  dz_stretch_level_end(9) = 9999999.9_wp            !< namelist parameter
+    REAL(wp) ::  dz_stretch_factor_array(9) = 1.08_wp              !< namelist parameter
     REAL(wp) ::  mask_scale(3)                                     !< collective array for mask_scale_x/y/z
     REAL(wp) ::  pt_vertical_gradient(10) = 0.0_wp                 !< namelist parameter
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 3058)
+++ palm/trunk/SOURCE/parin.f90	(revision 3065)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! New initialization parameters added
+! 
+! 3049 2018-05-29 13:52:36Z Giersch
 ! Error messages revised
 ! 
@@ -509,7 +512,7 @@
              dp_external, dp_level_b, dp_smooth, dpdxy, dry_aerosol_radius,    &
              dt, dt_pr_1d, dt_run_control_1d, dt_spinup, dx, dy, dz, dz_max,   &
-             dz_stretch_factor, dz_stretch_level, end_time_1d,                 &
-             ensemble_member_nr, e_init, e_min, fft_method,                    &
-             flux_input_mode, flux_output_mode, forcing,                       &
+             dz_stretch_factor, dz_stretch_level, dz_stretch_level_start,      &
+             dz_stretch_level_end, end_time_1d, ensemble_member_nr, e_init,    &
+             e_min, fft_method, flux_input_mode, flux_output_mode, forcing,    &
              galilei_transformation, humidity,                                 &
              inflow_damping_height, inflow_damping_width,                      &
@@ -579,7 +582,7 @@
              dp_external, dp_level_b, dp_smooth, dpdxy, dry_aerosol_radius,    &
              dt, dt_pr_1d, dt_run_control_1d, dt_spinup, dx, dy, dz, dz_max,   &
-             dz_stretch_factor, dz_stretch_level, end_time_1d,                 &
-             ensemble_member_nr, e_init, e_min, fft_method,                    &
-             flux_input_mode, flux_output_mode, forcing,                       &
+             dz_stretch_factor, dz_stretch_level, dz_stretch_level_start,      &
+             dz_stretch_level_end, end_time_1d, ensemble_member_nr, e_init,    &
+             e_min, fft_method, flux_input_mode, flux_output_mode, forcing,    &
              galilei_transformation, humidity,                                 &
              inflow_damping_height, inflow_damping_width,                      &
Index: palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 3058)
+++ palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 3065)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! dz was replaced by the help of zw to allow for vertical stretching
+! 
+! 3049 2018-05-29 13:52:36Z Giersch
 ! Error messages revised
 ! 
@@ -550,5 +553,5 @@
 
        USE control_parameters,                                                 &
-           ONLY: dz, passive_scalar
+           ONLY: passive_scalar
 
 
@@ -568,5 +571,5 @@
        REAL(wp) ::  canopy_height       !< canopy height (in m)
        
-       canopy_height = pch_index * dz
+       canopy_height = zw(pch_index)
 
        WRITE ( io, 1 )  canopy_mode, canopy_height, pch_index,                 &
@@ -668,5 +671,5 @@
 
        USE control_parameters,                                                 &
-           ONLY: dz, humidity, io_blocks, io_group, message_string, ocean,     &
+           ONLY: humidity, io_blocks, io_group, message_string, ocean,         &
                  passive_scalar, urban_surface
 
@@ -761,5 +764,5 @@
 !--       Use beta function for lad-profile construction
           int_bpdf = 0.0_wp
-          canopy_height = pch_index * dz
+          canopy_height = zw(pch_index)
 
           DO k = 0, pch_index
Index: palm/trunk/SOURCE/pmc_interface_mod.f90
===================================================================
--- palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 3058)
+++ palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 3065)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! dz was replaced by dz(1)
+! 
+! 3049 2018-05-29 13:52:36Z Giersch
 ! Error messages revised
 ! 
@@ -913,5 +916,5 @@
           parent_grid_info_real(5) = lower_left_coord_x + ( nx + 1 ) * dx
           parent_grid_info_real(6) = lower_left_coord_y + ( ny + 1 ) * dy
-          parent_grid_info_real(7) = dz
+          parent_grid_info_real(7) = dz(1)
 
           parent_grid_info_int(1)  = nx
@@ -1318,5 +1321,5 @@
        fval(3) = dx
        fval(4) = dy
-       fval(5) = dz
+       fval(5) = dz(1)
 
        IF ( myid == 0 )  THEN
Index: palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model_mod.f90	(revision 3058)
+++ palm/trunk/SOURCE/radiation_model_mod.f90	(revision 3065)
@@ -28,4 +28,8 @@
 ! -----------------
 ! $Id$
+! dz was replaced by dz(1), error message concerning vertical stretching was 
+! added  
+! 
+! 3049 2018-05-29 13:52:36Z Giersch
 ! Error messages revised
 ! 
@@ -330,4 +334,5 @@
 ! ------------
 !> Radiation models and interfaces
+!> @todo Replace dz(1) appropriatly to account for grid stretching
 !> @todo move variable definitions used in radiation_init only to the subroutine
 !>       as they are no longer required after initialization.
@@ -4371,5 +4376,5 @@
      IF ( plant_canopy )  THEN
          pchf_prep(:) = r_d * (hyp(nzub:nzut) / 100000.0_wp)**0.286_wp &
-                     / (cp * hyp(nzub:nzut) * dx*dy*dz) !< equals to 1 / (rho * c_p * Vbox * T)
+                     / (cp * hyp(nzub:nzut) * dx*dy*dz(1)) !< equals to 1 / (rho * c_p * Vbox * T)
          pctf_prep(:) = r_d * (hyp(nzub:nzut) / 100000.0_wp)**0.286_wp &
                      / (l_v * hyp(nzub:nzut) * dx*dy*dz)
@@ -4396,5 +4401,5 @@
 !--         each dimension, so that this new direction vector will allow us
 !--         to traverse the ray path within grid coordinates directly
-         sunorig_grid = (/ sunorig(1)/dz, sunorig(2)/dy, sunorig(3)/dx /)
+         sunorig_grid = (/ sunorig(1)/dz(1), sunorig(2)/dy, sunorig(3)/dx /)
 !--         sunorig_grid = sunorig_grid / norm2(sunorig_grid)
          sunorig_grid = sunorig_grid / SQRT(SUM(sunorig_grid**2))
@@ -4403,5 +4408,5 @@
 !--            precompute effective box depth with prototype Leaf Area Density
             pc_box_dimshift = MAXLOC(ABS(sunorig), 1) - 1
-            CALL box_absorb(CSHIFT((/dz,dy,dx/), pc_box_dimshift),      &
+            CALL box_absorb(CSHIFT((/dz(1),dy,dx/), pc_box_dimshift),      &
                                 60, prototype_lad,                          &
                                 CSHIFT(ABS(sunorig), pc_box_dimshift),      &
@@ -4801,5 +4806,5 @@
         IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx
         IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy
-        IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz
+        IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz(1)
      ENDDO
 !
@@ -5014,4 +5019,7 @@
     SUBROUTINE radiation_interaction_init
 
+       USE control_parameters,                                                 &
+           ONLY:  dz_stretch_level_start
+           
        USE netcdf_data_input_mod,                                              &
            ONLY:  leaf_area_density_f
@@ -5095,4 +5103,16 @@
        nzpt = nzptl
 #endif
+!
+!--    Stretching (non-uniform grid spacing) is not considered in the radiation 
+!--    model. Therefore, vertical stretching has to be applied above the area 
+!--    where the parts of the radiation model which assume constant grid spacing
+!--    are active. ABS (...) is required because the default value of 
+!--    dz_stretch_level_start is -9999999.9_wp (negative). 
+       IF ( ABS( dz_stretch_level_start(1) ) <= zw(nzut) ) THEN
+          WRITE( message_string, * ) 'The lowest level where vertical ',       &
+                                     'stretching is applied have to be ',      &
+                                     'greater than ', zw(nzut)
+          CALL message( 'init_grid', 'PA0496', 1, 2, 0, 6, 0 )
+       ENDIF 
 !
 !--    global number of urban and plant layers
@@ -5390,5 +5410,5 @@
             IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx
             IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy
-            IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz
+            IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz(1)
         ENDDO
 
@@ -5512,5 +5532,5 @@
 
 !--                 unit vector source -> target
-                    uv = (/ (ta(1)-sa(1))*dz, (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)
+                    uv = (/ (ta(1)-sa(1))*dz(1), (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)
                     sqdist = SUM(uv(:)**2)
                     uv = uv / SQRT(sqdist)
@@ -5677,5 +5697,5 @@
 
 !--               unit vector source -> target
-                  uv = (/ (ta(1)-sa(1))*dz, (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)
+                  uv = (/ (ta(1)-sa(1))*dz(1), (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)
                   sqdist = SUM(uv(:)**2)
                   uv = uv / SQRT(sqdist)
@@ -6127,5 +6147,5 @@
         ENDIF
         uvect(:) = delta(:) / distance
-        realdist = SQRT(SUM( (uvect(:)*(/dz,dy,dx/))**2 ))
+        realdist = SQRT(SUM( (uvect(:)*(/dz(1),dy,dx/))**2 ))
 
         lastdist = 0._wp
@@ -6509,5 +6529,5 @@
                      zb1 = CEILING(zexit * zsgn - .5_wp) - 1  ! because it must be smaller than exit
                      nz = MAX(zb1 - zb0 + 3, 2)
-                     rt2_dist(nz) = SQRT(((zexit-zorig)*dz)**2 + dxxyy)
+                     rt2_dist(nz) = SQRT(((zexit-zorig)*dz(1))**2 + dxxyy)
                      qdist = rt2_dist(nz) / (zexit-zorig)
                      rt2_dist(2:nz-1) = (/( ((REAL(l, wp) + .5_wp) * zsgn - zorig) * qdist , l = zb0, zb1 )/)
@@ -6567,5 +6587,5 @@
                   zb1 = CEILING(zexit * zsgn - .5_wp) - 1  ! because it must be smaller than exit
                   nz = MAX(zb1 - zb0 + 3, 2)
-                  rt2_dist(nz) = SQRT(((zexit-zorig)*dz)**2 + dxxyy)
+                  rt2_dist(nz) = SQRT(((zexit-zorig)*dz(1))**2 + dxxyy)
                   qdist = rt2_dist(nz) / (zexit-zorig)
                   rt2_dist(2:nz-1) = (/( ((REAL(l, wp) + .5_wp) * zsgn - zorig) * qdist , l = zb0, zb1 )/)
Index: palm/trunk/SOURCE/read_restart_data_mod.f90
===================================================================
--- palm/trunk/SOURCE/read_restart_data_mod.f90	(revision 3058)
+++ palm/trunk/SOURCE/read_restart_data_mod.f90	(revision 3065)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! New parameters concerning vertical grid stretching have been added
+! 
+! 3056 2018-06-04 07:49:35Z Giersch
 ! found variable has to be set to false inside overlap loop
 ! 
@@ -407,6 +410,12 @@
              CASE ( 'dz_stretch_factor' )
                 READ ( 13 )  dz_stretch_factor
+             CASE ( 'dz_stretch_factor_array' )
+                READ ( 13 )  dz_stretch_factor_array
              CASE ( 'dz_stretch_level' )
                 READ ( 13 )  dz_stretch_level
+             CASE ( 'dz_stretch_level_end' )
+                READ ( 13 )  dz_stretch_level_end
+             CASE ( 'dz_stretch_level_start' )
+                READ ( 13 )  dz_stretch_level_start
              CASE ( 'e_min' )
                 READ ( 13 )  e_min
@@ -1251,5 +1260,5 @@
     THEN
        WRITE( message_string, * ) 'number of PEs or virtual PE-grid changed ', &
-                        'in restart run  PE', myid, ' will read from files ', &
+                        'in restart run & PE', myid, ' will read from files ', &
                          file_list(1:files_to_be_opened)
        CALL message( 'rrd_local', 'PA0285', 0, 0, 0, 6, 0 )
Index: palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90
===================================================================
--- palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90	(revision 3058)
+++ palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90	(revision 3065)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Error message related to vertical stretching has been added, dz was replaced
+! by dz(1)
+! 
+! 3051 2018-05-30 17:43:55Z suehring
 ! Bugfix in calculation of initial Reynolds-stress tensor.
 ! 
@@ -313,5 +317,6 @@
 
     USE control_parameters,                                                    &
-        ONLY:  bc_lr, bc_ns, forcing, nest_domain, rans_mode, turbulent_inflow
+        ONLY:  bc_lr, bc_ns, forcing, nest_domain, number_stretch_level_start, &
+               rans_mode, turbulent_inflow
 
     USE pmc_interface,                                                         &
@@ -365,4 +370,10 @@
                               'is not allowed'
           CALL message( 'stg_check_parameters', 'PA0039', 1, 2, 0, 6, 0 )
+       ENDIF
+       
+       IF ( number_stretch_level_start > 0 )  THEN
+          message_string = 'Using synthetic turbulence generator ' //          &
+                           'in combination with stretching is not allowed'
+          CALL message( 'stg_check_parameters', 'PA0420', 1, 2, 0, 6, 0 )
        ENDIF
 
@@ -607,11 +618,12 @@
 
 !
-!--       Convert length scales from meter to number of grid points
+!--       Convert length scales from meter to number of grid points. Attention:
+!--       Does not work if grid stretching is used
           nuy(k) = INT( luy * ddy )
-          nuz(k) = INT( luz / dz  )
+          nuz(k) = INT( luz / dz(1)  )
           nvy(k) = INT( lvy * ddy )
-          nvz(k) = INT( lvz / dz  )
+          nvz(k) = INT( lvz / dz(1)  )
           nwy(k) = INT( lwy * ddy )
-          nwz(k) = INT( lwz / dz  )
+          nwz(k) = INT( lwz / dz(1)  )
 !
 !--       Workaround, assume isotropic turbulence
Index: palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3058)
+++ palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3065)
@@ -28,4 +28,8 @@
 ! -----------------
 ! $Id$
+! Unused array dxdir was removed, dz was replaced by dzu to consider vertical
+! grid stretching
+! 
+! 3049 2018-05-29 13:52:36Z Giersch
 ! Error messages revised
 ! 
@@ -282,5 +286,5 @@
 #if ! defined( __nopointer )
     USE arrays_3d,                                                             &
-        ONLY:  hyp, zu, pt, pt_1, pt_2, p, u, v, w, hyp, tend
+        ONLY:  dzu, hyp, zu, pt, pt_1, pt_2, p, u, v, w, hyp, tend
 #endif
 
@@ -292,5 +296,5 @@
     
     USE control_parameters,                                                    &
-        ONLY:  coupling_start_time, dz, topography, dt_3d,                     &
+        ONLY:  coupling_start_time, topography, dt_3d,                         &
                intermediate_timestep_count, initializing_actions,              &
                intermediate_timestep_count_max, simulated_time, end_time,      &
@@ -7207,5 +7211,4 @@
         REAL(wp), DIMENSION(nzb:nzt)          :: exn                !< value of the Exner function in layers
         
-        REAL(wp), DIMENSION(0:4)              :: dxdir              !< surface normal direction gridbox length
         REAL(wp)                              :: dtime              !< simulated time of day (in UTC)
         INTEGER(iwp)                          :: dhour              !< simulated hour of day (in UTC)
@@ -7213,5 +7216,4 @@
 
 
-        dxdir = (/dz,dy,dy,dx,dx/)
 #if ! defined( __nopointer )
         exn(nzb:nzt) = (hyp(nzb:nzt) / 100000.0_wp )**0.286_wp          !< Exner function
@@ -7634,5 +7636,5 @@
                          (dtime/3600.0_wp-REAL(dhour,wp))*aheatprof(k,dhour+1)
                  IF ( aheat(k,j,i) > 0.0_wp )  THEN
-                    pt(k,j,i) = pt(k,j,i) + aheat(k,j,i)*acoef*dt_3d/(exn(k)*rho_cp*dz)
+                    pt(k,j,i) = pt(k,j,i) + aheat(k,j,i)*acoef*dt_3d/(exn(k)*rho_cp*dzu(k))
                  ENDIF
               ENDIF
Index: palm/trunk/SOURCE/user_init_grid.f90
===================================================================
--- palm/trunk/SOURCE/user_init_grid.f90	(revision 3058)
+++ palm/trunk/SOURCE/user_init_grid.f90	(revision 3065)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! dz was replaced by dz(1)
+! 
+! 2718 2018-01-02 08:49:38Z maronga
 ! Corrected "Former revisions" section
 ! 
@@ -108,7 +111,7 @@
 !--       topography, bit is 1 for atmospheric grid point. 
 !--       The following example shows how to prescribe sine-like topography 
-!--       along x-direction with amplitude of 10 * dz and wavelength 10 * dy.
+!--       along x-direction with amplitude of 10 * dz(1) and wavelength 10 * dy.
 !           DO  i = nxlg, nxrg
-!              h_topo = 10.0_wp * dz * (SIN(3.14_wp*0.5_wp)*i*dx / ( 5.0_wp * dy ) )**2
+!              h_topo = 10.0_wp * dz(1) * (SIN(3.14_wp*0.5_wp)*i*dx / ( 5.0_wp * dy ) )**2
 ! 
 !              k_topo = MINLOC( ABS( zw - h_topo ), 1 ) - 1
Index: palm/trunk/SOURCE/vertical_nesting_mod.f90
===================================================================
--- palm/trunk/SOURCE/vertical_nesting_mod.f90	(revision 3058)
+++ palm/trunk/SOURCE/vertical_nesting_mod.f90	(revision 3065)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! dz was replaced by dz(1), error messages related to vertical grid stretching 
+! have been added
+! 
+! 3049 2018-05-29 13:52:36Z Giersch
 ! Error messages revised
 ! 
@@ -73,4 +77,5 @@
 !> after spin-up of the CG
 !>
+!> @todo Replace dz(1) appropriatly to account for grid stretching
 !> @todo Ensure that code can be compiled for serial and parallel mode. Please
 !>       check the placement of the directive "__parallel".
@@ -3639,5 +3644,6 @@
 
        USE control_parameters,                                                    &
-           ONLY:  coupling_mode, coupling_mode_remote, coupling_topology, dz
+           ONLY:  coupling_mode, coupling_mode_remote, coupling_topology, dz,     &
+                  dz_stretch_level_start, message_string
 
        USE grid_variables,                                                        &
@@ -3678,40 +3684,40 @@
            dxc = dx
            dyc = dy
-           dzc = dz
+           dzc = dz(1)
            cg_nprocs = numprocs
 
            IF ( myid == 0 )  THEN
 
-               CALL MPI_SEND( nxc, 1, MPI_INTEGER  , numprocs, 1,  comm_inter,  &
+               CALL MPI_SEND( nxc, 1, MPI_INTEGER  , numprocs, 1,  comm_inter, &
                    ierr )
-               CALL MPI_SEND( nyc, 1, MPI_INTEGER  , numprocs, 2,  comm_inter,  &
+               CALL MPI_SEND( nyc, 1, MPI_INTEGER  , numprocs, 2,  comm_inter, &
                    ierr )
-               CALL MPI_SEND( nzc, 1, MPI_INTEGER  , numprocs, 3,  comm_inter,  &
+               CALL MPI_SEND( nzc, 1, MPI_INTEGER  , numprocs, 3,  comm_inter, &
                    ierr )
-               CALL MPI_SEND( dxc, 1, MPI_REAL     , numprocs, 4,  comm_inter,  &
+               CALL MPI_SEND( dxc, 1, MPI_REAL     , numprocs, 4,  comm_inter, &
                    ierr )
-               CALL MPI_SEND( dyc, 1, MPI_REAL     , numprocs, 5,  comm_inter,  &
+               CALL MPI_SEND( dyc, 1, MPI_REAL     , numprocs, 5,  comm_inter, &
                    ierr )
-               CALL MPI_SEND( dzc, 1, MPI_REAL     , numprocs, 6,  comm_inter,  &
+               CALL MPI_SEND( dzc, 1, MPI_REAL     , numprocs, 6,  comm_inter, &
                    ierr )
-               CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 7,  comm_inter,  &
+               CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 7,  comm_inter, &
                    ierr )
                CALL MPI_SEND( cg_nprocs, 1, MPI_INTEGER, numprocs, 8,  comm_inter,  &
                    ierr )
-               CALL MPI_RECV( nxf, 1, MPI_INTEGER,   numprocs, 21, comm_inter,  &
+               CALL MPI_RECV( nxf, 1, MPI_INTEGER,   numprocs, 21, comm_inter, &
                    status, ierr )
-               CALL MPI_RECV( nyf, 1, MPI_INTEGER,   numprocs, 22, comm_inter,  &
+               CALL MPI_RECV( nyf, 1, MPI_INTEGER,   numprocs, 22, comm_inter, &
                    status, ierr )
-               CALL MPI_RECV( nzf, 1, MPI_INTEGER,   numprocs, 23, comm_inter,  &
+               CALL MPI_RECV( nzf, 1, MPI_INTEGER,   numprocs, 23, comm_inter, &
                    status, ierr )
-               CALL MPI_RECV( dxf, 1, MPI_REAL,      numprocs, 24, comm_inter,  &
+               CALL MPI_RECV( dxf, 1, MPI_REAL,      numprocs, 24, comm_inter, &
                    status, ierr )
-               CALL MPI_RECV( dyf, 1, MPI_REAL,      numprocs, 25, comm_inter,  &
+               CALL MPI_RECV( dyf, 1, MPI_REAL,      numprocs, 25, comm_inter, &
                    status, ierr )
-               CALL MPI_RECV( dzf, 1, MPI_REAL,      numprocs, 26, comm_inter,  &
+               CALL MPI_RECV( dzf, 1, MPI_REAL,      numprocs, 26, comm_inter, &
                    status, ierr )
-               CALL MPI_RECV( pdims_partner, 2, MPI_INTEGER,                     &
+               CALL MPI_RECV( pdims_partner, 2, MPI_INTEGER,                   &
                    numprocs, 27, comm_inter, status, ierr )
-               CALL MPI_RECV( fg_nprocs, 1, MPI_INTEGER,                     &
+               CALL MPI_RECV( fg_nprocs, 1, MPI_INTEGER,                       &
                    numprocs, 28, comm_inter, status, ierr )
            ENDIF
@@ -3725,4 +3731,13 @@
            CALL MPI_BCAST( pdims_partner, 2, MPI_INTEGER, 0, comm2d, ierr )
            CALL MPI_BCAST( fg_nprocs,  1, MPI_INTEGER, 0, comm2d, ierr )
+           
+!
+!--        Check if stretching is used within the nested domain. ABS(...) is 
+!--        necessary because of the default value of -9999999.9_wp (negative)
+           IF ( ABS( dz_stretch_level_start(1) ) <= (nzf+1)*dzf )  THEN       
+               message_string = 'Stretching in the parent domain is '//        &
+                                'only allowed above the nested domain'
+               CALL message( 'vertical_nesting_mod', 'PA0497', 1, 2, 0, 6, 0 )
+           ENDIF
 
        ELSEIF ( coupling_mode ==  'vnested_fine' )  THEN
@@ -3733,5 +3748,5 @@
            dxf = dx
            dyf = dy
-           dzf = dz
+           dzf = dz(1)
            fg_nprocs = numprocs
 
@@ -3774,5 +3789,5 @@
 
        ENDIF
-  
+       
        ngp_c = ( nxc+1 + 2 * nbgp ) * ( nyc+1 + 2 * nbgp )
        ngp_f = ( nxf+1 + 2 * nbgp ) * ( nyf+1 + 2 * nbgp )
@@ -3910,11 +3925,11 @@
  
 #if defined( __parallel )
-          USE arrays_3d,                                                             &
+          USE arrays_3d,                                                       &
               ONLY:  zu, zw
               
-          USE control_parameters,                                                    &
-              ONLY:  coupling_mode
+          USE control_parameters,                                              &
+              ONLY:  coupling_mode, message_string, number_stretch_level_start
              
-          USE indices,                                                               &
+          USE indices,                                                         &
               ONLY:  nzt
          
@@ -3925,5 +3940,6 @@
           IMPLICIT NONE
       
-          !-- Allocate and Exchange zuc and zuf, zwc and zwf
+!
+!--       Allocate and Exchange zuc and zuf, zwc and zwf
           IF ( coupling_mode(1:8)  == 'vnested_' )  THEN
           
@@ -3932,7 +3948,8 @@
           
              IF ( coupling_mode ==  'vnested_crse' )  THEN
-          
-                   zuc = zu
-                   zwc = zw
+                
+                zuc = zu
+                zwc = zw
+                
                 IF ( myid == 0 )  THEN
           
@@ -3954,6 +3971,15 @@
              ELSEIF ( coupling_mode ==  'vnested_fine' )  THEN
           
-                   zuf = zu
-                   zwf = zw
+!
+!--             Check if stretching is used within the nested domain
+                IF ( number_stretch_level_start > 0 )  THEN
+                   message_string = 'Stretching in the nested domain is not '//&
+                           'allowed'
+                   CALL message( 'vertical_nesting_mod', 'PA0498', 1, 2, 0, 6, 0 )
+                ENDIF
+                
+                zuf = zu
+                zwf = zw
+                
                 IF ( myid == 0 )  THEN
           
Index: palm/trunk/SOURCE/virtual_flight_mod.f90
===================================================================
--- palm/trunk/SOURCE/virtual_flight_mod.f90	(revision 3058)
+++ palm/trunk/SOURCE/virtual_flight_mod.f90	(revision 3065)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! IF statement revised to consider new vertical grid stretching
+! 
+! 3049 2018-05-29 13:52:36Z Giersch
 ! Error messages revised
 ! 
@@ -734,5 +737,5 @@
 
        USE control_parameters,                                                 &
-           ONLY:  dz, dz_stretch_level   
+           ONLY:  dz, dz_stretch_level_start   
   
       USE grid_variables,                                                     &
@@ -777,6 +780,6 @@
 !--    if flight position is above the vertical grid-stretching level.
 !--    fac=1 if variable is on scalar grid level, fac=0 for w-component.
-       IF ( z_pos(l) < dz_stretch_level )  THEN
-          k = ( z_pos(l) + fac * 0.5_wp * dz ) / dz
+       IF ( z_pos(l) < dz_stretch_level_start(1) )  THEN
+          k = ( z_pos(l) + fac * 0.5_wp * dz(1) ) / dz(1)
        ELSE
 !
Index: palm/trunk/SOURCE/wind_turbine_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/wind_turbine_model_mod.f90	(revision 3058)
+++ palm/trunk/SOURCE/wind_turbine_model_mod.f90	(revision 3065)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! dz was replaced by dz(1), error message concerning grid stretching has been
+! introduced
+! 
+! 3049 2018-05-29 13:52:36Z Giersch
 ! Error messages revised
 ! 
@@ -132,4 +136,5 @@
 !> automatically).
 !>
+!> @todo Replace dz(1) appropriatly to account for grid stretching
 !> @todo Revise code according to PALM Coding Standard
 !> @todo Implement ADM and ALM turbine models
@@ -740,5 +745,5 @@
        delta_r_factor = segment_width
        delta_t_factor = segment_length
-       delta_r_init   = delta_r_factor * MIN( dx, dy, dz)
+       delta_r_init   = delta_r_factor * MIN( dx, dy, dz(1))
 
        DO inot = 1, nturbines
@@ -937,4 +942,7 @@
 
     
+       USE control_parameters,                                                 &
+           ONLY:  dz_stretch_level_start
+    
        IMPLICIT NONE
 
@@ -992,4 +1000,17 @@
                  pi**( 1.0_wp / 6.0_wp ) * eps_kernel
 !
+!--    Stretching (non-uniform grid spacing) is not considered in the wind 
+!--    turbine model. Therefore, vertical stretching has to be applied above
+!--    the area where the wtm is active. ABS (...) is required because the 
+!--    default value of dz_stretch_level_start is -9999999.9_wp (negative). 
+       IF ( ABS( dz_stretch_level_start(1) ) <= MAXVAL(rcz) + MAXVAL(rr) +     &
+       eps_min)  THEN
+          WRITE( message_string, * ) 'The lowest level where vertical ',       &
+                                     'stretching is applied &have to be ',     &
+                                     'greater than ', MAXVAL(rcz) +            &
+                                      MAXVAL(rr) + eps_min
+          CALL message( 'init_grid', 'PA0484', 1, 2, 0, 6, 0 )
+       ENDIF 
+!
 !--    Square of eps_min:
        eps_min2 = eps_min**2
@@ -1022,5 +1043,5 @@
           i_hub(inot) = INT(   rcx(inot)                 / dx )
           j_hub(inot) = INT( ( rcy(inot) + 0.5_wp * dy ) / dy )
-          k_hub(inot) = INT( ( rcz(inot) + 0.5_wp * dz ) / dz )
+          k_hub(inot) = INT( ( rcz(inot) + 0.5_wp * dz(1) ) / dz(1) )
 
 !
@@ -1031,5 +1052,5 @@
           i_smear(inot) = CEILING( ( rr(inot) + eps_min ) / dx )
           j_smear(inot) = CEILING( ( rr(inot) + eps_min ) / dy )
-          k_smear(inot) = CEILING( ( rr(inot) + eps_min ) / dz )
+          k_smear(inot) = CEILING( ( rr(inot) + eps_min ) / dz(1) )
        
        ENDDO
@@ -1121,6 +1142,6 @@
 !--       surface. All points between z=0 and z=dz/s would already be contained
 !--       in grid box 1.
-          index_nacb(inot) = INT( ( rcz(inot) - rnac(inot) ) / dz ) + 1
-          index_nact(inot) = INT( ( rcz(inot) + rnac(inot) ) / dz ) + 1
+          index_nacb(inot) = INT( ( rcz(inot) - rnac(inot) ) / dz(1) ) + 1
+          index_nact(inot) = INT( ( rcz(inot) + rnac(inot) ) / dz(1) ) + 1
 
 !
@@ -1148,5 +1169,5 @@
                             IF ( j == tower_s )  THEN
                                tow_cd_surf(k,j,i) = ( rcz(inot) -              &
-                                    ( k_hub(inot) * dz - 0.5_wp * dz ) )  *    & ! extension in z-direction
+                                    ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) )*& ! extension in z-direction
                                   ( ( tower_s + 1.0_wp + 0.5_wp ) * dy    -    &
                                     ( rcy(inot) - 0.5_wp * dtow(inot) ) ) *    & ! extension in y-direction
@@ -1154,5 +1175,5 @@
                             ELSEIF ( j == tower_n )  THEN
                                tow_cd_surf(k,j,i) = ( rcz(inot)            -   &
-                                    ( k_hub(inot) * dz - 0.5_wp * dz ) )  *    & ! extension in z-direction
+                                    ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) )*& ! extension in z-direction
                                   ( ( rcy(inot) + 0.5_wp * dtow(inot) )   -    &
                                     ( tower_n + 0.5_wp ) * dy )           *    & ! extension in y-direction
@@ -1163,5 +1184,5 @@
                             ELSE
                                tow_cd_surf(k,j,i) = ( rcz(inot) -              &
-                                    ( k_hub(inot) * dz - 0.5_wp * dz ) )  *    &
+                                    ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) )*&
                                     dy * turb_cd_tower(inot)
                             ENDIF
@@ -1169,6 +1190,6 @@
 !--                      tower lies completely within one grid box:
                          ELSE
-                            tow_cd_surf(k,j,i) = ( rcz(inot)                 - &
-                                       ( k_hub(inot) * dz - 0.5_wp * dz ) ) *  &
+                            tow_cd_surf(k,j,i) = ( rcz(inot) - ( k_hub(inot) * &
+                                       dz(1) - 0.5_wp * dz(1) ) ) *            &
                                        dtow(inot) * turb_cd_tower(inot)
                          ENDIF
@@ -1182,10 +1203,10 @@
 !--                         leftmost and rightmost grid box:
                             IF ( j == tower_s )  THEN                         
-                               tow_cd_surf(k,j,i) = dz * (                     &
+                               tow_cd_surf(k,j,i) = dz(1) * (                  &
                                       ( tower_s + 1 + 0.5_wp ) * dy         -  &
                                       ( rcy(inot) - 0.5_wp * dtow(inot) )      &
                                                         ) * turb_cd_tower(inot)
                             ELSEIF ( j == tower_n )  THEN
-                               tow_cd_surf(k,j,i) = dz * (                     &
+                               tow_cd_surf(k,j,i) = dz(1) * (                  &
                                       ( rcy(inot) + 0.5_wp * dtow(inot) )   -  &
                                       ( tower_n + 0.5_wp ) * dy                &
@@ -1195,10 +1216,11 @@
 !--                         (where tow_cd_surf = grid box area):
                             ELSE
-                               tow_cd_surf(k,j,i) = dz * dy * turb_cd_tower(inot)
+                               tow_cd_surf(k,j,i) = dz(1) * dy *               &
+                                                    turb_cd_tower(inot)
                             ENDIF
 !
 !--                         tower lies completely within one grid box:
                          ELSE
-                            tow_cd_surf(k,j,i) = dz * dtow(inot) *             &
+                            tow_cd_surf(k,j,i) = dz(1) * dtow(inot) *          &
                                                 turb_cd_tower(inot)
                          ENDIF ! end if larger than grid box
@@ -1302,6 +1324,6 @@
        ENDDO   ! end of loop over turbines
 
-       tow_cd_surf   = tow_cd_surf   / ( dx * dy * dz )      ! Normalize tower drag
-       nac_cd_surf = nac_cd_surf / ( dx * dy * dz )      ! Normalize nacelle drag
+       tow_cd_surf   = tow_cd_surf   / ( dx * dy * dz(1) )  ! Normalize tower drag
+       nac_cd_surf = nac_cd_surf / ( dx * dy * dz(1) )      ! Normalize nacelle drag
 
        CALL wtm_read_blade_tables
@@ -1746,5 +1768,5 @@
                    ii =   rbx(ring,rseg) * ddx
                    jj = ( rby(ring,rseg) - 0.5_wp * dy ) * ddy
-                   kk = ( rbz(ring,rseg) + 0.5_wp * dz ) / dz
+                   kk = ( rbz(ring,rseg) + 0.5_wp * dz(1) ) / dz(1)
 !
 !--                Interpolate only if all required information is available on
@@ -1776,5 +1798,5 @@
 
                          u_int_1_l(inot,ring,rseg) = u_int_l          +        &
-                                     ( rbz(ring,rseg) - zu(kk) ) / dz *        &
+                                     ( rbz(ring,rseg) - zu(kk) ) / dz(1) *     &
                                      ( u_int_u - u_int_l )
 
@@ -1791,5 +1813,5 @@
                    ii = ( rbx(ring,rseg) - 0.5_wp * dx ) * ddx
                    jj =   rby(ring,rseg)                 * ddy
-                   kk = ( rbz(ring,rseg) + 0.5_wp * dz ) / dz 
+                   kk = ( rbz(ring,rseg) + 0.5_wp * dz(1) ) / dz(1) 
 !
 !--                Interpolate only if all required information is available on
@@ -1821,5 +1843,5 @@
 
                          v_int_1_l(inot,ring,rseg) = v_int_l +                 &
-                                     ( rbz(ring,rseg) - zu(kk) ) / dz *        &
+                                     ( rbz(ring,rseg) - zu(kk) ) / dz(1) *     &
                                      ( v_int_u - v_int_l )
 
@@ -1836,5 +1858,5 @@
                    ii = ( rbx(ring,rseg) - 0.5_wp * dx ) * ddx
                    jj = ( rby(ring,rseg) - 0.5_wp * dy ) * ddy
-                   kk =   rbz(ring,rseg)                 / dz
+                   kk =   rbz(ring,rseg)                 / dz(1)
 !
 !--                Interpolate only if all required information is available on
@@ -1866,5 +1888,5 @@
 
                          w_int_1_l(inot,ring,rseg) = w_int_l +                 &
-                                     ( rbz(ring,rseg) - zw(kk) ) / dz *        &
+                                     ( rbz(ring,rseg) - zw(kk) ) / dz(1) *     &
                                      ( w_int_u - w_int_l )
                       ELSE
@@ -2273,11 +2295,11 @@
                             dist_u_3d = ( i * dx               - rbx(ring,rseg) )**2 + &
                                         ( j * dy + 0.5_wp * dy - rby(ring,rseg) )**2 + &
-                                        ( k * dz - 0.5_wp * dz - rbz(ring,rseg) )**2
+                                        ( k * dz(1) - 0.5_wp * dz(1) - rbz(ring,rseg) )**2
                             dist_v_3d = ( i * dx + 0.5_wp * dx - rbx(ring,rseg) )**2 + &
                                         ( j * dy               - rby(ring,rseg) )**2 + &
-                                        ( k * dz - 0.5_wp * dz - rbz(ring,rseg) )**2
+                                        ( k * dz(1) - 0.5_wp * dz(1) - rbz(ring,rseg) )**2
                             dist_w_3d = ( i * dx + 0.5_wp * dx - rbx(ring,rseg) )**2 + &
                                         ( j * dy + 0.5_wp * dy - rby(ring,rseg) )**2 + &
-                                        ( k * dz               - rbz(ring,rseg) )**2
+                                        ( k * dz(1)               - rbz(ring,rseg) )**2
 
 !
Index: palm/trunk/SOURCE/write_restart_data_mod.f90
===================================================================
--- palm/trunk/SOURCE/write_restart_data_mod.f90	(revision 3058)
+++ palm/trunk/SOURCE/write_restart_data_mod.f90	(revision 3065)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! New parameters concerning vertical grid stretching have been added
+! 
+! 3004 2018-04-27 12:33:25Z Giersch
 ! precipitation_rate_av removed
 ! 
@@ -380,5 +383,5 @@
        CALL wrd_write_string( 'dz' ) 
        WRITE ( 14 )  dz
-
+       
        CALL wrd_write_string( 'dz_max' ) 
        WRITE ( 14 )  dz_max
@@ -386,8 +389,17 @@
        CALL wrd_write_string( 'dz_stretch_factor' ) 
        WRITE ( 14 )  dz_stretch_factor
-
+       
+       CALL wrd_write_string( 'dz_stretch_factor_array' ) 
+       WRITE ( 14 )  dz_stretch_factor_array
+       
        CALL wrd_write_string( 'dz_stretch_level' ) 
        WRITE ( 14 )  dz_stretch_level
 
+       CALL wrd_write_string( 'dz_stretch_level_end' ) 
+       WRITE ( 14 )  dz_stretch_level_end
+       
+       CALL wrd_write_string( 'dz_stretch_level_start' ) 
+       WRITE ( 14 )  dz_stretch_level_start
+       
        CALL wrd_write_string( 'e_min' ) 
        WRITE ( 14 )  e_min
