Index: /palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90
===================================================================
--- /palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90	(revision 4565)
+++ /palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90	(revision 4566)
@@ -25,4 +25,10 @@
 ! -----------------
 ! $Id$
+! - revise parametrization for reynolds-stress components, turbulent length- and time scales
+! - revise computation of velocity disturbances to be consistent to Xie and Castro (2008)
+! - change default value of time interval to adjust turbulence parametrization
+! - bugfix in computation of amplitude-tensor (vertical flux of horizontal momentum)
+! 
+! 4562 2020-06-12 08:38:47Z raasch
 ! Parts of r4559 re-formatted
 ! 
@@ -151,8 +157,6 @@
 !> parametrized within the boundary layer.
 !>
-!> @todo test restart
-!>       Enable cyclic_fill
+!> @todo Enable cyclic_fill
 !>       Implement turbulence generation for e and pt
-!> @todo Input of height-constant length scales via namelist
 !> @note <Enter notes on the module>
 !> @bug  Height information from input file is not used. Profiles from input must match with current
@@ -287,5 +291,10 @@
     INTEGER(iwp) ::  id_stg_right       !< right lateral boundary core id in case of turbulence generator
     INTEGER(iwp) ::  id_stg_south       !< south lateral boundary core id in case of turbulence generator
-    INTEGER(iwp) ::  mergp              !< maximum length scale (in gp)
+    INTEGER(iwp) ::  k_zi               !< vertical index of boundary-layer top
+    INTEGER(iwp) ::  mergp_limit = 1000 !< maximum length scale (in gp)
+    INTEGER(iwp) ::  mergp_x            !< maximum length scale in x (in gp)
+    INTEGER(iwp) ::  mergp_xy           !< maximum horizontal length scale (in gp)
+    INTEGER(iwp) ::  mergp_y            !< maximum length scale in y (in gp)
+    INTEGER(iwp) ::  mergp_z            !< maximum length scale in z (in gp)
     INTEGER(iwp) ::  nzb_x_stg          !< lower bound of z coordinate (required for transposing z on PEs along x)
     INTEGER(iwp) ::  nzt_x_stg          !< upper bound of z coordinate (required for transposing z on PEs along x)
@@ -326,4 +335,6 @@
 
 
+    LOGICAL ::  adjustment_step               = .FALSE.  !< control flag indicating that time and lenght scales have been updated and
+                                                         !< no time correlation to the timestep before should be considered
     LOGICAL ::  compute_velocity_seeds_local  = .TRUE.   !< switch to decide whether velocity seeds are computed locally
                                                          !< or if computation is distributed over several processes
@@ -336,10 +347,10 @@
     REAL(wp) ::  blend                     !< value to create gradually and smooth blending of Reynolds stress and length
                                            !< scales above the boundary layer
-    REAL(wp) ::  blend_coeff = -2.3_wp     !< coefficient used to ensure that blending functions decreases to 1/10 after
+    REAL(wp) ::  blend_coeff = -9.3_wp     !< coefficient used to ensure that blending functions decreases to 1/10 after
                                            !< one length scale above ABL top
     REAL(wp) ::  d_l                       !< blend_coeff/length_scale
-    REAL(wp) ::  dt_stg_adjust = 300.0_wp  !< time interval for adjusting turbulence statistics
+    REAL(wp) ::  d_nxy                     !< inverse of the total number of xy-grid points
+    REAL(wp) ::  dt_stg_adjust = 1800.0_wp !< time interval for adjusting turbulence statistics
     REAL(wp) ::  dt_stg_call = 0.0_wp      !< time interval for calling synthetic turbulence generator
-    REAL(wp) ::  length_scale              !< length scale, default is 8 x minimum grid spacing
     REAL(wp) ::  scale_l                   !< scaling parameter used for turbulence parametrization - Obukhov length
     REAL(wp) ::  scale_us                  !< scaling parameter used for turbulence parametrization - friction velocity
@@ -390,6 +401,6 @@
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fwo_yz  !< velocity seed for w at yz plane with new random number
 
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist_xz  !< imposed disturbances at north/south boundary
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist_yz  !< imposed disturbances at north/south boundary
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist_xz !< disturbances for parallel/crosswind/vertical component at north/south boundary
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist_yz !< disturbances for parallel/crosswind/vertical component  at north/south boundary
 
 !
@@ -756,4 +767,12 @@
        ENDIF
     ENDIF
+!
+!-- Assign initial profiles. Note, this is only required if turbulent inflow from the left is 
+!-- desired, not in case of any of the nesting (offline or self nesting) approaches.
+    IF ( .NOT. nesting_offline  .AND.  .NOT.  child_domain )  THEN 
+       u_init = mean_inflow_profiles(:,1)
+       v_init = mean_inflow_profiles(:,2)
+       e_init = MAXVAL( mean_inflow_profiles(:,5) )
+    ENDIF
 
     ALLOCATE ( a11(nzb:nzt+1), a21(nzb:nzt+1), a22(nzb:nzt+1),                                     &
@@ -766,4 +785,14 @@
                tu(nzb:nzt+1),  tv(nzb:nzt+1),  tw(nzb:nzt+1) )
 
+    r11 = 0.0_wp
+    r21 = 0.0_wp
+    r22 = 0.0_wp
+    r31 = 0.0_wp
+    r32 = 0.0_wp
+    r33 = 0.0_wp
+    tu  = 0.0_wp
+    tv  = 0.0_wp
+    tw  = 0.0_wp
+
     ALLOCATE ( dist_xz(nzb:nzt+1,nxl:nxr,3) )
     ALLOCATE ( dist_yz(nzb:nzt+1,nys:nyn,3) )
@@ -831,11 +860,6 @@
     ELSE
 !
-!--    Define length scale for the imposed turbulence, which is defined as 8 times the minimum grid
-!--    spacing
-       length_scale = 8.0_wp * MIN( dx, dy, MINVAL( dzw ) )
-!
-!--    Define constant to gradually decreasing length scales and Reynolds stress above ABL top.
-!--    Assure that no zero length scales are used.
-       d_l = blend_coeff / MAX( length_scale, dx, dy, MINVAL( dzw ) )
+!--    Precalulate the inversion number of xy-grid points - later used for mass conservation
+       d_nxy = 1.0_wp / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
 !
 !--    Set flag indicating that turbulence is parametrized
@@ -851,11 +875,8 @@
        CALL calc_scaling_variables
 !
-!--    Initialize lenght and time scales, which in turn are used to initialize the filter functions.
-       CALL calc_length_and_time_scale
-!
 !--    Parametrize Reynolds-stress tensor, diagonal elements as well as r21 (v'u'), r31 (w'u'),
 !--    r32 (w'v'). Parametrization follows Rotach et al. (1996) and is based on boundary-layer depth,
 !--    friction velocity and velocity scale.
-       CALL parametrize_reynolds_stress
+       CALL parametrize_turbulence
 !
 !--    Calculate coefficient matrix from Reynolds stress tensor (Lund rotation)
@@ -863,39 +884,7 @@
 
     ENDIF
-
-!
-!-- Assign initial profiles. Note, this is only required if turbulent inflow from the left is
-!-- desired, not in case of any of the nesting (offline or self nesting) approaches.
-    IF ( .NOT. nesting_offline  .AND.  .NOT.  child_domain )  THEN
-       u_init = mean_inflow_profiles(:,1)
-       v_init = mean_inflow_profiles(:,2)
-      !pt_init = mean_inflow_profiles(:,4)
-       e_init = MAXVAL( mean_inflow_profiles(:,5) )
-    ENDIF
-
-!
-!-- Define the size of the filter functions and allocate them.
-    mergp = 0
-
-    ! arrays must be large enough to cover the largest length scale
-    DO  k = nzb, nzt+1
-       j = MAX( ABS( nux(k) ), ABS( nuy(k) ), ABS( nuz(k) ),                                       &
-                ABS( nvx(k) ), ABS( nvy(k) ), ABS( nvz(k) ),                                       &
-                ABS( nwx(k) ), ABS( nwy(k) ), ABS( nwz(k) ) )
-       IF ( j > mergp )  mergp = j
-    ENDDO
-
-!     mergp  = 2 * mergp
-!     mergp = mergp
-
-    ALLOCATE ( bux(-mergp:mergp,nzb:nzt+1),                                                        &
-               buy(-mergp:mergp,nzb:nzt+1),                                                        &
-               buz(-mergp:mergp,nzb:nzt+1),                                                        &
-               bvx(-mergp:mergp,nzb:nzt+1),                                                        &
-               bvy(-mergp:mergp,nzb:nzt+1),                                                        &
-               bvz(-mergp:mergp,nzb:nzt+1),                                                        &
-               bwx(-mergp:mergp,nzb:nzt+1),                                                        &
-               bwy(-mergp:mergp,nzb:nzt+1),                                                        &
-               bwz(-mergp:mergp,nzb:nzt+1) )
+!
+!-- Initial filter functions
+    CALL stg_setup_filter_function
 
 !
@@ -924,16 +913,4 @@
     fw_yz  = 0.0_wp
     fwo_yz = 0.0_wp
-
-!
-!-- Create filter functions
-    CALL stg_filter_func( nux, bux ) !filter ux
-    CALL stg_filter_func( nuy, buy ) !filter uy
-    CALL stg_filter_func( nuz, buz ) !filter uz
-    CALL stg_filter_func( nvx, bvx ) !filter vx
-    CALL stg_filter_func( nvy, bvy ) !filter vy
-    CALL stg_filter_func( nvz, bvz ) !filter vz
-    CALL stg_filter_func( nwx, bwx ) !filter wx
-    CALL stg_filter_func( nwy, bwy ) !filter wy
-    CALL stg_filter_func( nwz, bwz ) !filter wz
 
 #if defined( __parallel )
@@ -1059,22 +1036,20 @@
 !-- core. In case there is only inflow from the left, it is sufficient to generate only random
 !-- numbers for the yz-layer, else random numbers for the xz-layer are also required.
-    ALLOCATE ( id_rand_yz(-mergp+nys:nyn+mergp) )
-    ALLOCATE ( seq_rand_yz(5,-mergp+nys:nyn+mergp) )
+    ALLOCATE ( id_rand_yz(-mergp_limit+nys:nyn+mergp_limit) )
+    ALLOCATE ( seq_rand_yz(5,-mergp_limit+nys:nyn+mergp_limit) )
     id_rand_yz  = 0
     seq_rand_yz = 0
 
-    CALL init_parallel_random_generator( ny, -mergp+nys, nyn+mergp, id_rand_yz, seq_rand_yz )
+    CALL init_parallel_random_generator( ny, -mergp_limit+nys, nyn+mergp_limit, id_rand_yz, seq_rand_yz )
 
     IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent                             &
          .AND.  .NOT.  rans_mode ) )  THEN
-       ALLOCATE ( id_rand_xz(-mergp+nxl:nxr+mergp) )
-       ALLOCATE ( seq_rand_xz(5,-mergp+nxl:nxr+mergp) )
+       ALLOCATE ( id_rand_xz(-mergp_limit+nxl:nxr+mergp_limit) )
+       ALLOCATE ( seq_rand_xz(5,-mergp_limit+nxl:nxr+mergp_limit) )
        id_rand_xz  = 0
        seq_rand_xz = 0
 
-       CALL init_parallel_random_generator( nx, -mergp+nxl, nxr+mergp, id_rand_xz, seq_rand_xz )
+       CALL init_parallel_random_generator( nx, -mergp_limit+nxl, nxr+mergp_limit, id_rand_xz, seq_rand_xz )
     ENDIF
-
-
 
  END SUBROUTINE stg_init
@@ -1086,9 +1061,10 @@
 !> Calculate filter function bxx from length scale nxx following Eg.9 and 10 (Xie and Castro, 2008)
 !--------------------------------------------------------------------------------------------------!
- SUBROUTINE stg_filter_func( nxx, bxx )
-
-    INTEGER(iwp) ::  k    !< loop index
-    INTEGER(iwp) ::  n_k  !< length scale nXX in height k
-    INTEGER(iwp) ::  nf   !< index for length scales
+ SUBROUTINE stg_filter_func( nxx, bxx, mergp )
+
+    INTEGER(iwp) ::  k     !< loop index
+    INTEGER(iwp) ::  mergp !< passed length scale in grid points
+    INTEGER(iwp) ::  n_k   !< length scale nXX in height k
+    INTEGER(iwp) ::  nf    !< index for length scales
 
     INTEGER(iwp), DIMENSION(nzb:nzt+1) ::  nxx  !< length scale (in gp)
@@ -1125,4 +1101,68 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculate filter function bxx from length scale nxx following Eg.9 and 10
+!> (Xie and Castro, 2008)
+!------------------------------------------------------------------------------!
+ SUBROUTINE stg_setup_filter_function
+
+    INTEGER(iwp) :: j  !< dummy value to calculate maximum length scale index
+    INTEGER(iwp) :: k  !< loop index along vertical direction
+!
+!-- Define the size of the filter functions and allocate them.
+    mergp_x = 0
+    mergp_y = 0
+    mergp_z = 0
+
+    DO  k = nzb, nzt+1
+       j = MAX( ABS(nux(k)), ABS(nvx(k)), ABS(nwx(k)) )
+       IF ( j > mergp_x )  mergp_x = j
+    ENDDO
+    DO  k = nzb, nzt+1
+       j = MAX( ABS(nuy(k)), ABS(nvy(k)), ABS(nwy(k)) )
+       IF ( j > mergp_y )  mergp_y = j
+    ENDDO
+    DO  k = nzb, nzt+1
+       j = MAX( ABS(nuz(k)), ABS(nvz(k)), ABS(nwz(k)) )
+       IF ( j > mergp_z )  mergp_z = j
+    ENDDO
+
+    mergp_xy = MAX( mergp_x, mergp_y )
+
+    IF ( ALLOCATED( bux ) )  DEALLOCATE( bux )
+    IF ( ALLOCATED( bvx ) )  DEALLOCATE( bvx )
+    IF ( ALLOCATED( bwx ) )  DEALLOCATE( bwx )
+    IF ( ALLOCATED( buy ) )  DEALLOCATE( buy )
+    IF ( ALLOCATED( bvy ) )  DEALLOCATE( bvy )
+    IF ( ALLOCATED( bwy ) )  DEALLOCATE( bwy )
+    IF ( ALLOCATED( buz ) )  DEALLOCATE( buz )
+    IF ( ALLOCATED( bvz ) )  DEALLOCATE( bvz )
+    IF ( ALLOCATED( bwz ) )  DEALLOCATE( bwz )
+
+    ALLOCATE ( bux(-mergp_x:mergp_x,nzb:nzt+1),                                                    &
+               buy(-mergp_y:mergp_y,nzb:nzt+1),                                                    &
+               buz(-mergp_z:mergp_z,nzb:nzt+1),                                                    &
+               bvx(-mergp_x:mergp_x,nzb:nzt+1),                                                    &
+               bvy(-mergp_y:mergp_y,nzb:nzt+1),                                                    &
+               bvz(-mergp_z:mergp_z,nzb:nzt+1),                                                    &
+               bwx(-mergp_x:mergp_x,nzb:nzt+1),                                                    &
+               bwy(-mergp_y:mergp_y,nzb:nzt+1),                                                    &
+               bwz(-mergp_z:mergp_z,nzb:nzt+1) )
+!
+!-- Create filter functions
+    CALL stg_filter_func( nux, bux, mergp_x ) !filter ux
+    CALL stg_filter_func( nuy, buy, mergp_y ) !filter uy
+    CALL stg_filter_func( nuz, buz, mergp_z ) !filter uz
+    CALL stg_filter_func( nvx, bvx, mergp_x ) !filter vx
+    CALL stg_filter_func( nvy, bvy, mergp_y ) !filter vy
+    CALL stg_filter_func( nvz, bvz, mergp_z ) !filter vz
+    CALL stg_filter_func( nwx, bwx, mergp_x ) !filter wx
+    CALL stg_filter_func( nwy, bwy, mergp_y ) !filter wy
+    CALL stg_filter_func( nwz, bwz, mergp_z ) !filter wz
+
+ END SUBROUTINE stg_setup_filter_function
+
 !--------------------------------------------------------------------------------------------------!
 ! Description:
@@ -1262,5 +1302,5 @@
     LOGICAL  ::  stg_call = .FALSE.  !< control flag indicating whether turbulence was updated or only restored from last call
 
-    REAL(wp) ::  dt_stg  !< wheighted subtimestep
+    REAL(wp) ::  dt_stg  !< time interval the STG is called
 
     REAL(wp), DIMENSION(3) ::  mc_factor_l  !< local mass flux correction factor
@@ -1344,5 +1384,5 @@
 !
 !--             Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
-                IF ( tu(k) == 0.0_wp )  THEN
+                IF ( tu(k) == 0.0_wp  .OR.  adjustment_step )  THEN
                    fu_yz(k,j) = fuo_yz(k,j)
                 ELSE
@@ -1351,5 +1391,5 @@
                 ENDIF
 
-                IF ( tv(k) == 0.0_wp )  THEN
+                IF ( tv(k) == 0.0_wp  .OR.  adjustment_step )  THEN
                    fv_yz(k,j) = fvo_yz(k,j)
                 ELSE
@@ -1358,5 +1398,5 @@
                 ENDIF
 
-                IF ( tw(k) == 0.0_wp )  THEN
+                IF ( tw(k) == 0.0_wp  .OR.  adjustment_step )  THEN
                    fw_yz(k,j) = fwo_yz(k,j)
                 ELSE
@@ -1377,7 +1417,7 @@
 !
 !--             Lund rotation following Eq. 17 in Xie and Castro (2008).
-!--             Additional factors are added to improve the variance of v and w
-                dist_yz(k,j,1) = MIN( a11(k) * fu_yz(k,j), 3.0_wp ) *                              &
+                dist_yz(k,j,1) = a11(k) * fu_yz(k,j) *                                             &
                                  MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
+
              ENDDO
           ENDDO
@@ -1391,11 +1431,10 @@
 !--             Additional factors are added to improve the variance of v and w experimental test
 !--             of 1.2
-                dist_yz(k,j,2) = MIN( ( SQRT( a22(k) / MAXVAL( a22 ) ) * 1.2_wp )                  &
-                                 * ( a21(k) * fu_yz(k,j) + a22(k) * fv_yz(k,j) ), 3.0_wp ) *       &
+                dist_yz(k,j,2) = ( a21(k) * fu_yz(k,j) + a22(k) * fv_yz(k,j) ) *                   &
                                  MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
-                dist_yz(k,j,3) = MIN( ( SQRT( a33(k) / MAXVAL( a33 ) ) * 1.3_wp ) *                &
-                                      ( a31(k) * fu_yz(k,j) + a32(k) * fv_yz(k,j) + a33(k)         &
-                                        * fw_yz(k,j) ), 3.0_wp ) *                                 &
-                                      MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
+                dist_yz(k,j,3) = ( a31(k) * fu_yz(k,j) + a32(k) * fv_yz(k,j) +                     &
+                                   a33(k) * fw_yz(k,j) ) *                                         &
+                                 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
+
              ENDDO
           ENDDO
@@ -1510,5 +1549,5 @@
 !
 !--             Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
-                IF ( tu(k) == 0.0_wp )  THEN
+                IF ( tu(k) == 0.0_wp  .OR.  adjustment_step )  THEN
                    fu_xz(k,i) = fuo_xz(k,i)
                 ELSE
@@ -1517,5 +1556,5 @@
                 ENDIF
 
-                IF ( tv(k) == 0.0_wp )  THEN
+                IF ( tv(k) == 0.0_wp  .OR.  adjustment_step )  THEN
                    fv_xz(k,i) = fvo_xz(k,i)
                 ELSE
@@ -1524,5 +1563,5 @@
                 ENDIF
 
-                IF ( tw(k) == 0.0_wp )  THEN
+                IF ( tw(k) == 0.0_wp  .OR.  adjustment_step )  THEN
                    fw_xz(k,i) = fwo_xz(k,i)
                 ELSE
@@ -1546,6 +1585,5 @@
 !--             Additional factors are added to improve the variance of v and w
                 !experimental test of 1.2
-                dist_xz(k,i,2) = MIN( ( SQRT( a22(k) / MAXVAL( a22 ) ) * 1.2_wp )                  &
-                                 * ( a21(k) * fu_xz(k,i) + a22(k) * fv_xz(k,i) ), 3.0_wp ) *       &
+                dist_xz(k,i,2) = ( a21(k) * fu_xz(k,i) + a22(k) * fv_xz(k,i) ) *                   &
                                  MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
              ENDDO
@@ -1559,10 +1597,9 @@
 !--             Lund rotation following Eq. 17 in Xie and Castro (2008).
 !--             Additional factors are added to improve the variance of v and w
-                dist_xz(k,i,1) = MIN( a11(k) * fu_xz(k,i), 3.0_wp ) *                              &
+                dist_xz(k,i,1) = a11(k) * fu_xz(k,i) *                                             &
                                  MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
 
-                dist_xz(k,i,3) = MIN( ( SQRT(a33(k) / MAXVAL( a33 ) ) * 1.3_wp )                   &
-                                 * ( a31(k) * fu_xz(k,i) + a32(k) * fv_xz(k,i) + a33(k)            &
-                                 * fw_xz(k,i) ), 3.0_wp ) *                                        &
+                dist_xz(k,i,3) = ( a31(k) * fu_xz(k,i) + a32(k) * fv_xz(k,i) +                     &
+                                   a33(k) * fw_xz(k,i) ) *                                         &
                                  MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
              ENDDO
@@ -1654,4 +1691,8 @@
 !-- Finally, set time counter for calling STG to zero
     IF ( stg_call )  time_stg_call = 0.0_wp
+!
+!-- Set adjustment step to False to indicate that time correlation can be
+!-- switched-on again.
+    adjustment_step = .FALSE.
 
     IF ( debug_output_timestep )  CALL debug_message( 'stg_main', 'end' )
@@ -1684,6 +1725,6 @@
     REAL(wp) ::  rand_sigma_inv  !< inverse of stdev of random number
 
-    REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1) ::  b_y  !< filter function in y-direction
-    REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1) ::  b_z  !< filter function in z-direction
+    REAL(wp), DIMENSION(-mergp_y:mergp_y,nzb:nzt+1) :: b_y     !< filter function in y-direction
+    REAL(wp), DIMENSION(-mergp_z:mergp_z,nzb:nzt+1) :: b_z     !< filter function in z-direction
 
     REAL(wp), DIMENSION(nzb_x_stg:nzt_x_stg+1,nys:nyn) ::  f_n_l  !<  local velocity seed
@@ -1696,21 +1737,21 @@
 !-- must be defined globally in order to compute the correlation matrices. However, the calculation
 !-- and normalization of random numbers is done locally, while the result is later distributed to
-!-- all processes. Further, please note, a set of random numbers is only calculated for the left
-!-- boundary, while the right boundary uses the same random numbers and thus also computes the same
+!-- all processes. Further, please note, a set of random numbers is only calculated for the west
+!-- boundary, while the east boundary uses the same random numbers and thus also computes the same
 !-- correlation matrix.
-    ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp+nys:nyn+mergp) )
+    ALLOCATE( rand_it(nzb-mergp_z:nzt+1+mergp_z,-mergp_y+nys:nyn+mergp_y) )
     rand_it = 0.0_wp
 
     rand_av        = 0.0_wp
     rand_sigma_inv = 0.0_wp
-    nyz_inv        = 1.0_wp / REAL( ( nzt + 1 + mergp - ( nzb - mergp ) + 1 )                      &
-                     * ( ny + mergp - ( 0 - mergp ) + 1 ), KIND = wp )
+    nyz_inv        = 1.0_wp / REAL( ( nzt + 1 + mergp_z - ( nzb - mergp_z ) + 1 )                  &
+                     * ( ny + mergp_y - ( 0 - mergp_y ) + 1 ), KIND = wp )
 !
 !-- Compute and normalize random numbers.
-    DO  j = nys - mergp, nyn + mergp
+    DO  j = nys - mergp_y, nyn + mergp_y
 !
 !--    Put the random seeds at grid point j
        CALL random_seed_parallel( put=seq_rand_yz(:,j) )
-       DO  k = nzb - mergp, nzt + 1 + mergp
+       DO  k = nzb - mergp_z, nzt + 1 + mergp_z
           CALL random_number_parallel( random_dummy )
           rand_it(k,j) = random_dummy
@@ -1724,10 +1765,10 @@
 !-- random numbers, the inner ghost layers mergp must not be summed-up, else the random numbers on
 !-- the ghost layers will be stronger weighted as they also occur on the inner subdomains.
-    DO  j = MERGE( nys, nys - mergp, nys /= 0 ), MERGE( nyn, nyn + mergp, nyn /= ny )
-       DO  k = nzb - mergp, nzt + 1 + mergp
+    DO  j = MERGE( nys, nys - mergp_y, nys /= 0 ), MERGE( nyn, nyn + mergp_y, nyn /= ny )
+       DO  k = nzb - mergp_z, nzt + 1 + mergp_z
           rand_av = rand_av + rand_it(k,j)
        ENDDO
     ENDDO
-
+    
 #if defined( __parallel )
 !
@@ -1741,6 +1782,6 @@
 !
 !-- Now, compute the variance
-    DO  j = MERGE( nys, nys - mergp, nys /= 0 ), MERGE( nyn, nyn + mergp, nyn /= ny )
-       DO  k = nzb - mergp, nzt + 1 + mergp
+    DO  j = MERGE( nys, nys - mergp_y, nys /= 0 ), MERGE( nyn, nyn + mergp_y, nyn /= ny )
+       DO  k = nzb - mergp_z, nzt + 1 + mergp_z
           rand_sigma_inv = rand_sigma_inv + rand_it(k,j)**2
        ENDDO
@@ -1842,26 +1883,26 @@
  SUBROUTINE stg_generate_seed_xz( n_x, n_z, b_x, b_z, f_n, id_south, id_north )
 
-    INTEGER(iwp) ::  i           !< loop index in x-direction
-    INTEGER(iwp) ::  id_north    !< core ids at respective boundaries
-    INTEGER(iwp) ::  id_south    !< core ids at respective boundaries
-    INTEGER(iwp) ::  ii          !< loop index in x-direction
-    INTEGER(iwp) ::  k           !< loop index in z-direction
-    INTEGER(iwp) ::  kk          !< loop index in z-direction
-    INTEGER(iwp) ::  send_count  !< send count for MPI_GATHERV
-
-    INTEGER(iwp), DIMENSION(nzb:nzt+1) ::  n_x  !< length scale in x-direction
-    INTEGER(iwp), DIMENSION(nzb:nzt+1) ::  n_z  !< length scale in z-direction
-
-    REAL(wp) ::  nxz_inv         !< inverse of number of grid points in xz-slice
-    REAL(wp) ::  rand_av         !< average of random number
-    REAL(wp) ::  rand_sigma_inv  !< inverse of stdev of random number
-
-    REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1) ::  b_x  !< filter function in x-direction
-    REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1) ::  b_z  !< filter function in z-direction
-
-    REAL(wp), DIMENSION(nzb_y_stg:nzt_y_stg+1,nxl:nxr) ::  f_n_l  !<  local velocity seed
-    REAL(wp), DIMENSION(nzb:nzt+1,nxl:nxr)             ::  f_n    !<  velocity seed
-
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rand_it  !< global array of random numbers
+    INTEGER(iwp) :: i           !< loop index in x-direction
+    INTEGER(iwp) :: id_north    !< core ids at respective boundaries
+    INTEGER(iwp) :: id_south    !< core ids at respective boundaries
+    INTEGER(iwp) :: ii          !< loop index in x-direction
+    INTEGER(iwp) :: k           !< loop index in z-direction
+    INTEGER(iwp) :: kk          !< loop index in z-direction
+    INTEGER(iwp) :: send_count  !< send count for MPI_GATHERV
+
+    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x    !< length scale in x-direction
+    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z    !< length scale in z-direction
+
+    REAL(wp) :: nxz_inv         !< inverse of number of grid points in xz-slice
+    REAL(wp) :: rand_av         !< average of random number
+    REAL(wp) :: rand_sigma_inv  !< inverse of stdev of random number
+
+    REAL(wp), DIMENSION(-mergp_x:mergp_x,nzb:nzt+1) :: b_x     !< filter function in x-direction
+    REAL(wp), DIMENSION(-mergp_z:mergp_z,nzb:nzt+1) :: b_z     !< filter function in z-direction
+    
+    REAL(wp), DIMENSION(nzb_y_stg:nzt_y_stg+1,nxl:nxr) :: f_n_l   !<  local velocity seed
+    REAL(wp), DIMENSION(nzb:nzt+1,nxl:nxr)             :: f_n     !<  velocity seed
+
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rand_it   !< global array of random numbers
 
 !
@@ -1870,21 +1911,21 @@
 !-- must be defined globally in order to compute the correlation matrices. However, the calculation
 !-- and normalization of random numbers is done locally, while the result is later distributed to
-!-- all processes. Further, please note, a set of random numbers is only calculated for the left
-!-- boundary, while the right boundary uses the same random numbers and thus also computes the same
+!-- all processes. Further, please note, a set of random numbers is only calculated for the south
+!-- boundary, while the north boundary uses the same random numbers and thus also computes the same
 !-- correlation matrix.
-    ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp+nxl:nxr+mergp) )
+    ALLOCATE( rand_it(nzb-mergp_z:nzt+1+mergp_z,-mergp_x+nxl:nxr+mergp_x) )
     rand_it = 0.0_wp
 
     rand_av        = 0.0_wp
     rand_sigma_inv = 0.0_wp
-    nxz_inv        = 1.0_wp / REAL( ( nzt + 1 + mergp - ( nzb - mergp ) + 1 )                      &
-                     * ( nx + mergp - ( 0 - mergp ) +1 ), KIND = wp )
+    nxz_inv        = 1.0_wp / REAL( ( nzt + 1 + mergp_z - ( nzb - mergp_z ) + 1 )                  &
+                     * ( nx + mergp_x - ( 0 - mergp_x ) +1 ), KIND = wp )
 !
 !-- Compute and normalize random numbers.
-    DO  i = nxl - mergp, nxr + mergp
+    DO  i = nxl - mergp_x, nxr + mergp_x
 !
 !--    Put the random seeds at grid point ii
        CALL random_seed_parallel( put=seq_rand_xz(:,i) )
-       DO  k = nzb - mergp, nzt + 1 + mergp
+       DO  k = nzb - mergp_z, nzt + 1 + mergp_z
           CALL random_number_parallel( random_dummy )
           rand_it(k,i) = random_dummy
@@ -1899,10 +1940,10 @@
 !-- summed-up, else the random numbers on the ghost layers will be stronger weighted as they
 !-- also occur on the inner subdomains.
-    DO  i = MERGE( nxl, nxl - mergp, nxl /= 0 ), MERGE( nxr, nxr + mergp, nxr /= nx )
-       DO  k = nzb - mergp, nzt + 1 + mergp
+    DO  i = MERGE( nxl, nxl - mergp_x, nxl /= 0 ), MERGE( nxr, nxr + mergp_x, nxr /= nx )
+       DO  k = nzb - mergp_z, nzt + 1 + mergp_z
           rand_av = rand_av + rand_it(k,i)
        ENDDO
     ENDDO
-
+    
 #if defined( __parallel )
 !
@@ -1916,6 +1957,6 @@
 !
 !-- Now, compute the variance
-    DO  i = MERGE( nxl, nxl - mergp, nxl /= 0 ), MERGE( nxr, nxr + mergp, nxr /= nx )
-       DO  k = nzb - mergp, nzt + 1 + mergp
+    DO  i = MERGE( nxl, nxl - mergp_x, nxl /= 0 ), MERGE( nxr, nxr + mergp_x, nxr /= nx )
+       DO  k = nzb - mergp_z, nzt + 1 + mergp_z
           rand_sigma_inv = rand_sigma_inv + rand_it(k,i)**2
        ENDDO
@@ -2004,47 +2045,179 @@
 ! Description:
 ! ------------
-!> Parametrization of the Reynolds stress tensor, following the parametrization described in Rotach
-!> et al. (1996), which is applied in state-of-the-art dispserion modelling. Please note, the
-!> parametrization does not distinguish between along-wind and cross-wind turbulence.
-!--------------------------------------------------------------------------------------------------!
- SUBROUTINE parametrize_reynolds_stress
-
-    INTEGER(iwp) ::  k    !< loop index in z-direction
-
-    REAL(wp)     ::  zzi  !< ratio of z/zi
-
-    DO  k = nzb+1, nzt+1
+!> Parametrization of the Reynolds-stress componentes, turbulent length- and time scales. The 
+!> parametrization follows Brost et al. (1982) with modifications described in Rotach et al. (1996),
+!> which is applied in state-of-the-art dispserion modelling.
+!--------------------------------------------------------------------------------------------------!
+  SUBROUTINE parametrize_turbulence
+
+    INTEGER(iwp) :: k                 !< loop index in z-direction
+
+    REAL(wp) ::  corr_term_uh         !< correction term in parametrization of horizontal variances for unstable stratification
+    REAL(wp) ::  d_zi                 !< inverse boundary-layer depth
+    REAL(wp) ::  length_scale_lon     !< length scale in flow direction
+    REAL(wp) ::  length_scale_lon_zi  !< length scale in flow direction at boundary-layer top
+    REAL(wp) ::  length_scale_lat     !< length scale in crosswind direction
+    REAL(wp) ::  length_scale_lat_zi  !< length scale in crosswind direction at boundary-layer top
+    REAL(wp) ::  length_scale_vert    !< length scale in vertical direction
+    REAL(wp) ::  length_scale_vert_zi !< length scale in vertical direction at boundary-layer top
+    REAL(wp) ::  time_scale_zi        !< time scale at boundary-layer top
+    REAL(wp) ::  r11_zi               !< longitudinal variance at boundary-layer top
+    REAL(wp) ::  r22_zi               !< crosswind variance at boundary-layer top
+    REAL(wp) ::  r33_zi               !< vertical variance at boundary-layer top
+    REAL(wp) ::  r31_zi               !< r31 at boundary-layer top
+    REAL(wp) ::  r32_zi               !< r32 at boundary-layer top
+    REAL(wp) ::  rlat1                !< first dummy argument for crosswind compoment of reynolds stress
+    REAL(wp) ::  rlat2                !< second dummy argument for crosswind compoment of reynolds stress
+    REAL(wp) ::  rlon1                !< first dummy argument for longitudinal compoment of reynolds stress
+    REAL(wp) ::  rlon2                !< second dummy argument for longitudinal compoment of reynolds stress
+    REAL(wp) ::  zzi                  !< ratio of z/zi
+
+    REAL(wp), DIMENSION(nzb+1:nzt+1) ::  cos_phi    !< cosine of angle between mean flow direction and x-axis
+    REAL(wp), DIMENSION(nzb+1:nzt+1) ::  phi        !< angle between mean flow direction and x-axis
+    REAL(wp), DIMENSION(nzb+1:nzt+1) ::  sin_phi    !< sine of angle between mean flow direction and x-axis
+
+!
+!-- Calculate the boundary-layer top index. Boundary-layer top is calculated by Richardson-bulk
+!-- criterion according to Heinze et al. (2017).
+    k_zi = MAX( 1, MINLOC( ABS( zu - zi_ribulk ), DIM = 1 ) )
+    IF ( zu(k_zi) > zi_ribulk )  k_zi = k_zi - 1
+    k_zi = MIN( k_zi, nzt )
+!
+!-- Calculate angle between flow direction and x-axis.
+    DO  k = nzb+1, nzt + 1
+       IF ( u_init(k) /= 0.0_wp )  THEN
+          phi(k) = ATAN( v_init(k) / u_init(k) )
+       ELSE
+          phi(k) = 0.5 * pi
+       ENDIF
+!
+!--    Pre-calculate sine and cosine.
+       cos_phi(k) = COS( phi(k) )
+       sin_phi(k) = SIN( phi(k) )
+    ENDDO
+!
+!-- Parametrize Reynolds-stress components. Please note, parametrization is formulated for the
+!-- stream- and spanwise components, while stream- and spanwise direction do not necessarily
+!-- coincide with the grid axis. Hence, components are rotated after computation.
+    d_zi = 1.0_wp / zi_ribulk
+    DO  k = nzb+1, k_zi
+
+       corr_term_uh = MERGE( 0.35_wp * ABS(                                                        &
+                        - zi_ribulk / ( kappa * scale_l - 10E-4_wp )                               &
+                                       )**( 2.0_wp / 3.0_wp ),                                     &
+                          0.0_wp,                                                                  &
+                          scale_l < 0.0_wp )
+!
+!--    Determine normalized height coordinate
+       zzi = zu(k) * d_zi
+!
+!--    Compute longitudinal and crosswind component of reynolds stress. Rotate these components by
+!--    the wind direction to map onto the u- and v-component. Note, since reynolds stress for
+!--    variances cannot become negative, take the absolute values.
+       rlon1 = scale_us**2 * ( corr_term_uh + 5.0_wp - 4.0_wp * zzi )
+       rlat1 = scale_us**2 * ( corr_term_uh + 2.0_wp - zzi )
+
+       r11(k) = ABS(   cos_phi(k) * rlon1 + sin_phi(k) * rlat1 )
+       r22(k) = ABS( - sin_phi(k) * rlon1 + cos_phi(k) * rlat1 )
+!
+!--    For u'v' no parametrization exist so far - ?. For simplicity assume a similar profile as
+!--    for u'w'.
+       r21(k) = 0.5_wp * ( r11(k) + r22(k) )
+!
+!--    w'w'
+       r33(k) = scale_wm**2 * (                                                                    &
+                   1.5_wp * zzi**( 2.0_wp / 3.0_wp ) * EXP( -2.0_wp * zzi )                        &
+                 + ( 1.7_wp - zzi ) * ( scale_us / scale_wm )**2                                   &
+                              )
+!
+!--    u'w' and v'w'. After calculation of the longitudinal and crosswind component
+!--    these are projected along the x- and y-direction. Note, it is assumed that 
+!--    the flux within the boundary points opposite to the vertical gradient.
+       rlon2 = scale_us**2 * ( zzi - 1.0_wp )
+       rlat2 = scale_us**2 * ( 0.4 * zzi * ( 1.0_wp - zzi ) )
+
+       r31(k) = SIGN( ABS(   cos_phi(k) * rlon2 + sin_phi(k) * rlat2 ),                            &
+                   - ( u_init(k+1) - u_init(k-1) ) )
+       r32(k) = SIGN( ABS( - sin_phi(k) * rlon2 + cos_phi(k) * rlat2 ),                            &
+                   - ( v_init(k+1) - v_init(k-1) ) )
+!
+!--    Compute turbulent time scales according to Brost et al. (1982). Note, time scales are
+!--    limited to the adjustment time scales.
+       tu(k)  = MIN( dt_stg_adjust,                                                                &
+                     3.33_wp * zzi * ( 1.0 - 0.67_wp * zzi ) / scale_wm * zi_ribulk )
+
+       tv(k)  = tu(k)
+       tw(k)  = tu(k)
+
+       length_scale_lon  = MIN( SQRT( r11(k) ) * tu(k), zi_ribulk )
+       length_scale_lat  = MIN( SQRT( r22(k) ) * tv(k), zi_ribulk )
+       length_scale_vert = MIN( SQRT( r33(k) ) * tw(k), zi_ribulk )
+!
+!--    Assume isotropic turbulence length scales
+       nux(k) = MAX( INT( length_scale_lon  * ddx     ), 1 )
+       nuy(k) = MAX( INT( length_scale_lat  * ddy     ), 1 )
+       nvx(k) = MAX( INT( length_scale_lon  * ddx     ), 1 )
+       nvy(k) = MAX( INT( length_scale_lat  * ddy     ), 1 )
+       nwx(k) = MAX( INT( length_scale_lon  * ddx     ), 1 )
+       nwy(k) = MAX( INT( length_scale_lat  * ddy     ), 1 )
+       nuz(k) = MAX( INT( length_scale_vert * ddzw(k) ), 1 )
+       nvz(k) = MAX( INT( length_scale_vert * ddzw(k) ), 1 )
+       nwz(k) = MAX( INT( length_scale_vert * ddzw(k) ), 1 )
+
+    ENDDO
+!
+!-- Above boundary-layer top length- and timescales as well as reynolds-stress components are
+!-- reduced to zero by a blending function.
+    length_scale_lon_zi  = SQRT( r11(k_zi) ) * tu(k_zi)
+    length_scale_lat_zi  = SQRT( r22(k_zi) ) * tv(k_zi)
+    length_scale_vert_zi = SQRT( r33(k_zi) ) * tu(k_zi)
+
+    time_scale_zi        = tu(k_zi)
+    r11_zi               = r11(k_zi)
+    r22_zi               = r22(k_zi)
+    r33_zi               = r33(k_zi)
+    r31_zi               = r31(k_zi)
+    r32_zi               = r32(k_zi)
+
+    d_l = blend_coeff / MAX( length_scale_vert_zi, dx, dy, MINVAL( dzw ) )
+
+    DO  k = k_zi+1, nzt+1
 !
 !--    Calculate function to gradually decrease Reynolds stress above ABL top.
-!--    The function decreases to 1/10 after one length scale above the ABL top.
        blend = MIN( 1.0_wp, EXP( d_l * zu(k) - d_l * zi_ribulk ) )
-!
-!--    Determine normalized height coordinate
-       zzi = zu(k) / zi_ribulk
 !
 !--    u'u' and v'v'. Assume isotropy. Note, add a small negative number to the denominator, else
 !--    the mergpe-function can crash if scale_l is zero.
-       r11(k) = scale_us**2 * ( MERGE( 0.35_wp *                                                   &
-                         ABS( - zi_ribulk / ( kappa * scale_l - 10E-4_wp ) )**( 2.0_wp / 3.0_wp ), &
-                         0.0_wp, scale_l < 0.0_wp ) + 5.0_wp - 4.0_wp * zzi                                 &
-                              ) * blend
-
-       r22(k) = r11(k)
-!
-!--    w'w'
-       r33(k) = scale_wm**2 * ( 1.5_wp * zzi**( 2.0_wp / 3.0_wp ) * EXP( -2.0_wp * zzi )           &
-                + ( 1.7_wp - zzi ) * ( scale_us / scale_wm )**2 )  * blend
-!
-!--    u'w' and v'w'. Assume isotropy.
-       r31(k) = - scale_us**2 * ( 1.01_wp - MIN( zzi, 1.0_wp ) ) * blend
-       r32(k) = r31(k)
-!
-!--    For u'v' no parametrization exist so far - ?. For simplicity assume a similar profile as for
-!--    u'w'.
-       r21(k) = r31(k)
+       r11(k) = r11_zi * blend
+       r22(k) = r22_zi * blend
+       r33(k) = r33_zi * blend
+       r31(k) = r31_zi * blend
+       r32(k) = r32_zi * blend
+       r21(k) = 0.5_wp * ( r11(k) + r22(k) )
+
+!
+!--    Compute turbulent time scales according to Brost et al. (1982).
+!--    Note, time scales are limited to the adjustment time scales.
+       tu(k)  = time_scale_zi * blend
+       tv(k)  = time_scale_zi * blend
+       tw(k)  = time_scale_zi * blend
+
+       length_scale_lon  = length_scale_lon_zi  * blend
+       length_scale_lat  = length_scale_lat_zi  * blend
+       length_scale_vert = length_scale_vert_zi * blend
+!
+!--    Assume isotropic turbulence length scales
+       nux(k) = MAX( INT( length_scale_lon  * ddx ),     1 )
+       nuy(k) = MAX( INT( length_scale_lat  * ddy ),     1 )
+       nvx(k) = MAX( INT( length_scale_lon  * ddx ),     1 )
+       nvy(k) = MAX( INT( length_scale_lat  * ddy ),     1 )
+       nwx(k) = MAX( INT( length_scale_lon  * ddx ),     1 )
+       nwy(k) = MAX( INT( length_scale_lat  * ddy ),     1 )
+       nuz(k) = MAX( INT( length_scale_vert * ddzw(k) ), 1 )
+       nvz(k) = MAX( INT( length_scale_vert * ddzw(k) ), 1 )
+       nwz(k) = MAX( INT( length_scale_vert * ddzw(k) ), 1 )
     ENDDO
-
-!
-!-- Set bottom boundary condition
+!
+!-- Set bottom boundary condition for reynolds stresses
     r11(nzb) = r11(nzb+1)
     r22(nzb) = r22(nzb+1)
@@ -2055,6 +2228,23 @@
     r32(nzb) = r32(nzb+1)
 
-
- END SUBROUTINE parametrize_reynolds_stress
+!
+!-- Set bottom boundary condition for the length and time scales
+    nux(nzb) = nux(nzb+1)
+    nuy(nzb) = nuy(nzb+1)
+    nuz(nzb) = nuz(nzb+1)
+    nvx(nzb) = nvx(nzb+1)
+    nvy(nzb) = nvy(nzb+1)
+    nvz(nzb) = nvz(nzb+1)
+    nwx(nzb) = nwx(nzb+1)
+    nwy(nzb) = nwy(nzb+1)
+    nwz(nzb) = nwz(nzb+1)
+
+    tu(nzb)  = tu(nzb+1)
+    tv(nzb)  = tv(nzb+1)
+    tw(nzb)  = tw(nzb+1)
+
+    adjustment_step = .TRUE.
+
+ END SUBROUTINE parametrize_turbulence
 
 !--------------------------------------------------------------------------------------------------!
@@ -2071,5 +2261,5 @@
     DO  k = nzb+1, nzt+1
        IF ( r11(k) > 10E-6_wp )  THEN
-          a11(k) = SQRT( r11(k) )
+          a11(k) = SQRT( r11(k) ) 
           a21(k) = r21(k) / a11(k)
           a31(k) = r31(k) / a11(k)
@@ -2081,9 +2271,9 @@
     ENDDO
     DO  k = nzb+1, nzt+1
-       a22(k) = r22(k) - a21(k)**2
+       a22(k) = r22(k) - a21(k)**2 
        IF ( a22(k) > 10E-6_wp )  THEN
           a22(k) = SQRT( a22(k) )
-          a32(k) = r32(k) - a21(k) * a31(k) / a22(k)
-       ELSE
+          a32(k) = ( r32(k) - a21(k) * a31(k) ) / a22(k)
+       ELSE 
           a22(k) = 10E-8_wp
           a32(k) = 10E-8_wp
@@ -2128,25 +2318,15 @@
     CALL calc_scaling_variables
 !
-!-- Set length and time scales depending on boundary-layer height
-    CALL calc_length_and_time_scale
-!
-!-- Parametrize Reynolds-stress tensor, diagonal elements as well as r21 (v'u'), r31 (w'u'),
-!-- r32 (w'v'). Parametrization follows Rotach et al. (1996) and is based on boundary-layer depth,
-!-- friction velocity and velocity scale.
-    CALL parametrize_reynolds_stress
-!
-!-- Calculate coefficient matrix from Reynolds stress tensor (Lund rotation)
+!-- Parametrize Reynolds-stress tensor. Parametrization follows Brost et al. (1982) with 
+!-- modifications described in Rotach et al. (1996) and is based on boundary-layer depth, friction
+!-- velocity and velocity scale. 
+    CALL parametrize_turbulence
+!
+!-- Calculate coefficient matrix from Reynolds stress tensor  
+!-- (Lund rotation)
     CALL calc_coeff_matrix
 !
-!-- Determine filter functions on basis of updated length scales
-    CALL stg_filter_func( nux, bux ) !filter ux
-    CALL stg_filter_func( nuy, buy ) !filter uy
-    CALL stg_filter_func( nuz, buz ) !filter uz
-    CALL stg_filter_func( nvx, bvx ) !filter vx
-    CALL stg_filter_func( nvy, bvy ) !filter vy
-    CALL stg_filter_func( nvz, bvz ) !filter vz
-    CALL stg_filter_func( nwx, bwx ) !filter wx
-    CALL stg_filter_func( nwy, bwy ) !filter wy
-    CALL stg_filter_func( nwz, bwz ) !filter wz
+!-- Setup filter functions according to updated length scales.
+    CALL stg_setup_filter_function
 !
 !-- Reset time counter for controlling next adjustment to zero
@@ -2157,82 +2337,4 @@
  END SUBROUTINE stg_adjust
 
-
-!--------------------------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Calculates turbuluent length and time scales if these are not available from measurements.
-!--------------------------------------------------------------------------------------------------!
- SUBROUTINE calc_length_and_time_scale
-
-    INTEGER(iwp) ::  k  !< loop index in z-direction
-
-    REAL(wp) ::  length_scale_dum  !< effectively used length scale
-
-!
-!-- In initial call the boundary-layer depth can be zero. This case, set minimum value for
-!-- boundary-layer depth, to setup length scales correctly.
-    zi_ribulk = MAX( zi_ribulk, zw(nzb+2) )
-!
-!-- Set-up default turbulent length scales. From the numerical point of view the imposed
-!-- perturbations should not be immediately dissipated by the numerics. The numerical dissipation,
-!-- however, acts on scales up to 8 x the grid spacing. For this reason, set the turbulence length
-!-- scale to 8 time the grid spacing. Further, above the boundary layer height, set turbulence
-!-- lenght scales to zero (equivalent to not imposing any perturbations) in order to save
-!-- computational costs. Typical time scales are derived by assuming Taylors's hypothesis, using
-!-- the length scales and the mean profiles of u- and v-component.
-    DO  k = nzb+1, nzt+1
-!
-!--    Determine blending parameter. Within the boundary layer length scales are constant, while
-!--    above lengths scales approach gradully zero, i.e. imposed turbulence is not correlated in
-!--    space and time, just white noise, which saves computations power as the loops for the
-!--    computation of the filter functions depend on the length scales. The value decreases to
-!--    1/10 after one length scale above the ABL top.
-       blend = MIN( 1.0_wp, EXP( d_l * zu(k) - d_l * zi_ribulk ) )
-!
-!--    Assume isotropic turbulence length scales
-       nux(k) = MAX( INT( length_scale * ddx ), 1 ) * blend
-       nuy(k) = MAX( INT( length_scale * ddy ), 1 ) * blend
-       nvx(k) = MAX( INT( length_scale * ddx ), 1 ) * blend
-       nvy(k) = MAX( INT( length_scale * ddy ), 1 ) * blend
-       nwx(k) = MAX( INT( length_scale * ddx ), 1 ) * blend
-       nwy(k) = MAX( INT( length_scale * ddy ), 1 ) * blend
-!
-!--    Along the vertical direction limit the length scale further by the boundary-layer depth to
-!--    assure that no length scales larger than the boundary-layer depth are used
-       length_scale_dum = MIN( length_scale, zi_ribulk )
-
-       nuz(k) = MAX( INT( length_scale_dum * ddzw(k) ), 1 ) * blend
-       nvz(k) = MAX( INT( length_scale_dum * ddzw(k) ), 1 ) * blend
-       nwz(k) = MAX( INT( length_scale_dum * ddzw(k) ), 1 ) * blend
-!
-!--    Limit time scales, else they become very larger for low wind speed, imposing long-living
-!--    inflow perturbations which in turn propagate further into the model domain. Use u_init and
-!--    v_init to calculate the time scales, which will be equal to the inflow profiles, both, in
-!--    offline nesting mode or in dirichlet/radiation mode.
-       tu(k) = MIN( dt_stg_adjust, length_scale / ( ABS( u_init(k) ) + 0.1_wp ) ) * blend
-       tv(k) = MIN( dt_stg_adjust, length_scale / ( ABS( v_init(k) ) + 0.1_wp ) ) * blend
-!
-!--    Time scale of w-component is a mixture from u- and v-component.
-       tw(k)  = SQRT( tu(k)**2 + tv(k)**2 ) * blend
-
-    ENDDO
-!
-!-- Set bottom boundary condition for the length and time scales
-    nux(nzb) = nux(nzb+1)
-    nuy(nzb) = nuy(nzb+1)
-    nuz(nzb) = nuz(nzb+1)
-    nvx(nzb) = nvx(nzb+1)
-    nvy(nzb) = nvy(nzb+1)
-    nvz(nzb) = nvz(nzb+1)
-    nwx(nzb) = nwx(nzb+1)
-    nwy(nzb) = nwy(nzb+1)
-    nwz(nzb) = nwz(nzb+1)
-
-    tu(nzb)  = tu(nzb+1)
-    tv(nzb)  = tv(nzb+1)
-    tw(nzb)  = tw(nzb+1)
-
-
- END SUBROUTINE calc_length_and_time_scale
 
 !--------------------------------------------------------------------------------------------------!
@@ -2305,8 +2407,8 @@
 #endif
 
-    scale_us     = scale_us     / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
-    shf_mean     = shf_mean     / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
-    scale_l      = scale_l      / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
-    pt_surf_mean = pt_surf_mean / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
+    scale_us     = scale_us     * d_nxy
+    shf_mean     = shf_mean     * d_nxy
+    scale_l      = scale_l      * d_nxy
+    pt_surf_mean = pt_surf_mean * d_nxy
 !
 !-- Compute mean convective velocity scale. Note, in case the mean heat flux is negative, set
@@ -2318,4 +2420,9 @@
     ENDIF
 !
+!-- At the initial run the friction velocity is initialized close to zero, leading to almost zero
+!-- disturbances at the boundaries. In order to obtain disturbances nevertheless, set a minium
+!-- value of friction velocity for the reynolds-stress parametrization.
+    IF ( time_since_reference_point <= 0.0_wp )  scale_us = MAX( scale_us, 0.05_wp )
+!
 !-- Finally, in order to consider also neutral or stable stratification, compute momentum velocity
 !-- scale from u* and convective velocity scale, according to Rotach et al. (1996).
