Index: palm/trunk/SOURCE/init_grid.f90
===================================================================
--- palm/trunk/SOURCE/init_grid.f90	(revision 3051)
+++ 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
