Index: palm/trunk/SOURCE/nesting_offl_mod.f90
===================================================================
--- palm/trunk/SOURCE/nesting_offl_mod.f90	(revision 4020)
+++ palm/trunk/SOURCE/nesting_offl_mod.f90	(revision 4022)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! Detection of boundary-layer depth in stable boundary layer on basis of 
+! boundary data improved 
+! Routine for boundary-layer depth calculation renamed and made public
+! 
+! 3987 2019-05-22 09:52:13Z kanani
 ! Introduce alternative switch for debug output during timestepping
 ! 
@@ -117,6 +122,11 @@
 !
 !-- Public subroutines
-    PUBLIC nesting_offl_bc, nesting_offl_check_parameters, nesting_offl_header,&
-           nesting_offl_init, nesting_offl_mass_conservation, nesting_offl_parin 
+    PUBLIC nesting_offl_bc,                                                    &
+           nesting_offl_calc_zi,                                               &
+           nesting_offl_check_parameters,                                      &
+           nesting_offl_header,                                                &
+           nesting_offl_init,                                                  &
+           nesting_offl_mass_conservation,                                     &
+           nesting_offl_parin 
 !
 !-- Public variables
@@ -126,4 +136,8 @@
        MODULE PROCEDURE nesting_offl_bc
     END INTERFACE nesting_offl_bc
+    
+    INTERFACE nesting_offl_calc_zi
+       MODULE PROCEDURE nesting_offl_calc_zi
+    END INTERFACE nesting_offl_calc_zi
     
     INTERFACE nesting_offl_check_parameters
@@ -806,5 +820,5 @@
 !--    Further, adjust Rayleigh damping height in case of time-changing conditions.
 !--    Therefore, calculate boundary-layer depth first. 
-       CALL calc_zi
+       CALL nesting_offl_calc_zi
        CALL adjust_sponge_layer 
     
@@ -833,5 +847,5 @@
 !> bulk-Richardson criterion. 
 !------------------------------------------------------------------------------!
-    SUBROUTINE calc_zi
+    SUBROUTINE nesting_offl_calc_zi
        
        USE basic_constants_and_equations_mod,                                  &
@@ -845,8 +859,11 @@
        IMPLICIT NONE 
 
-       INTEGER(iwp) :: i            !< loop index in x-direction
-       INTEGER(iwp) :: j            !< loop index in y-direction
-       INTEGER(iwp) :: k            !< loop index in z-direction
-       INTEGER(iwp) :: k_surface    !< topography top index in z-direction
+       INTEGER(iwp) :: i                            !< loop index in x-direction
+       INTEGER(iwp) :: j                            !< loop index in y-direction
+       INTEGER(iwp) :: k                            !< loop index in z-direction
+       INTEGER(iwp) :: k_max_loc                    !< index of maximum wind speed along z-direction
+       INTEGER(iwp) :: k_surface                    !< topography top index in z-direction
+       INTEGER(iwp) :: num_boundary_gp_non_cyclic   !< number of non-cyclic boundaries, used for averaging ABL depth
+       INTEGER(iwp) :: num_boundary_gp_non_cyclic_l !< number of non-cyclic boundaries, used for averaging ABL depth
     
        REAL(wp) ::  ri_bulk                 !< bulk Richardson number
@@ -854,6 +871,4 @@
        REAL(wp) ::  topo_max                !< maximum topography level in model domain
        REAL(wp) ::  topo_max_l              !< maximum topography level in subdomain
-       REAL(wp) ::  u_comp                  !< u-component
-       REAL(wp) ::  v_comp                  !< v-component
        REAL(wp) ::  vpt_surface             !< near-surface virtual potential temperature
        REAL(wp) ::  zi_l                    !< mean boundary-layer depth on subdomain   
@@ -861,4 +876,5 @@
     
        REAL(wp), DIMENSION(nzb:nzt+1) ::  vpt_col !< vertical profile of virtual potential temperature at (j,i)-grid point
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  uv_abs  !< vertical profile of horizontal wind speed at (j,i)-grid point
 
        
@@ -867,5 +883,11 @@
 !--    Start with the left and right boundaries. 
        zi_l      = 0.0_wp
+       num_boundary_gp_non_cyclic_l = 0
        IF ( bc_dirichlet_l  .OR.  bc_dirichlet_r )  THEN
+!
+!--       Sum-up and store number of boundary grid points used for averaging
+!--       ABL depth
+          num_boundary_gp_non_cyclic_l = num_boundary_gp_non_cyclic_l +        &
+                                         nxr - nxl + 1
 !
 !--       Determine index along x. Please note, index indicates boundary
@@ -898,13 +920,22 @@
 !--          the boundary, where velocity inside the building is zero 
 !--          (k_surface is the index of the lowest upward-facing surface).
+             uv_abs(:) = SQRT( MERGE( u(:,j,i+1), u(:,j,i),                   &
+                                      bc_dirichlet_l )**2 +                   &
+                               v(:,j,i)**2 )
+!
+!--          Determine index of the maximum wind speed                               
+             k_max_loc = MAXLOC( uv_abs(:), DIM = 1 ) - 1
+
              zi_local = 0.0_wp
              DO  k = k_surface+1, nzt
-                u_comp = MERGE( u(k,j,i+1), u(k,j,i), bc_dirichlet_l )
-                v_comp = v(k,j,i)
                 ri_bulk = zu(k) * g / vpt_surface *                            &
                           ( vpt_col(k) - vpt_surface ) /                       &
-                          ( u_comp * u_comp + v_comp * v_comp + 1E-5_wp ) 
-                       
-                IF ( zi_local == 0.0_wp  .AND.  ri_bulk > ri_bulk_crit )       &
+                          ( uv_abs(k) + 1E-5_wp ) 
+!
+!--             Check if critical Richardson number is exceeded. Further, check
+!--             if there is a maxium in the wind profile in order to detect also
+!--             ABL heights in the stable boundary layer. 
+                IF ( zi_local == 0.0_wp  .AND.                                 &
+                     ( ri_bulk > ri_bulk_crit  .OR.  k == k_max_loc ) )        &
                    zi_local = zu(k)           
              ENDDO
@@ -920,4 +951,7 @@
 !--    Do the same at the north and south boundaries.
        IF ( bc_dirichlet_s  .OR.  bc_dirichlet_n )  THEN
+       
+          num_boundary_gp_non_cyclic_l = num_boundary_gp_non_cyclic_l +        &
+                                         nxr - nxl + 1
        
           j = MERGE( -1, nyn + 1, bc_dirichlet_s )
@@ -935,13 +969,22 @@
              ENDIF
           
+             uv_abs(:) = SQRT( u(:,j,i)**2 +                                   &
+                               MERGE( v(:,j+1,i), v(:,j,i),                    &
+                               bc_dirichlet_s )**2 )
+!
+!--          Determine index of the maximum wind speed
+             k_max_loc = MAXLOC( uv_abs(:), DIM = 1 ) - 1
+          
              zi_local = 0.0_wp
              DO  k = k_surface+1, nzt                
-                u_comp = u(k,j,i)
-                v_comp = MERGE( v(k,j+1,i), v(k,j,i), bc_dirichlet_s )
                 ri_bulk = zu(k) * g / vpt_surface *                            &
                        ( vpt_col(k) - vpt_surface ) /                          &
-                       ( u_comp * u_comp + v_comp * v_comp + 1E-5_wp ) 
-                       
-                IF ( zi_local == 0.0_wp  .AND.  ri_bulk > 0.25_wp )            &
+                       ( uv_abs(k) + 1E-5_wp ) 
+!
+!--             Check if critical Richardson number is exceeded. Further, check
+!--             if there is a maxium in the wind profile in order to detect also
+!--             ABL heights in the stable boundary layer.                        
+                IF ( zi_local == 0.0_wp  .AND.                                 &
+                     ( ri_bulk > ri_bulk_crit  .OR.  k == k_max_loc ) )        &
                    zi_local = zu(k)    
              ENDDO
@@ -955,8 +998,12 @@
        CALL MPI_ALLREDUCE( zi_l, zi_ribulk, 1, MPI_REAL, MPI_SUM,              &
                            comm2d, ierr )
+       CALL MPI_ALLREDUCE( num_boundary_gp_non_cyclic_l,                       &
+                           num_boundary_gp_non_cyclic,                         &
+                           1, MPI_INTEGER, MPI_SUM, comm2d, ierr )
 #else
        zi_ribulk = zi_l
+       num_boundary_gp_non_cyclic = num_boundary_gp_non_cyclic_l
 #endif
-       zi_ribulk = zi_ribulk / REAL( 2 * nx + 2 * ny, KIND = wp )
+       zi_ribulk = zi_ribulk / REAL( num_boundary_gp_non_cyclic, KIND = wp )
 !
 !--    Finally, check if boundary layer depth is not below the any topography. 
@@ -968,12 +1015,12 @@
        
 #if defined( __parallel )
-       CALL MPI_ALLREDUCE( topo_max_l, topo_max, 1, MPI_REAL, MPI_MAX,            &
+       CALL MPI_ALLREDUCE( topo_max_l, topo_max, 1, MPI_REAL, MPI_MAX,         &
                            comm2d, ierr )
 #else
        topo_max     = topo_max_l
 #endif
-       zi_ribulk = MAX( zi_ribulk, topo_max )
-       
-    END SUBROUTINE calc_zi
+!        zi_ribulk = MAX( zi_ribulk, topo_max )
+       
+    END SUBROUTINE nesting_offl_calc_zi
     
     
@@ -1319,5 +1366,5 @@
 !--    boundary data. This is requiered for initialize the synthetic turbulence
 !--    generator correctly. 
-       CALL calc_zi
+       CALL nesting_offl_calc_zi
        
 !
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 4020)
+++ palm/trunk/SOURCE/parin.f90	(revision 4022)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Change default top boundary condition for pressure to Neumann in offline
+! nesting case
+! 
+! 4017 2019-06-06 12:16:46Z schwenkel
 ! Introduce alternative switch for debug output during timestepping
 ! 
@@ -983,14 +987,7 @@
              bc_pt_t  = 'nesting_offline'
              bc_q_t   = 'nesting_offline'
-           !  bc_s_t   = 'nesting_offline'  ! scalar boundary condition is not clear
+           !  bc_s_t   = 'nesting_offline'  ! scalar boundary condition is not clear yet
            !  bc_cs_t  = 'nesting_offline'  ! same for chemical species
-!
-!--          For the pressure set Dirichlet conditions, in contrast to the 
-!--          self nesting. This gives less oscilations within the 
-!--          free atmosphere since the pressure solver has more degrees of 
-!--          freedom. In constrast to the self nesting, this might be 
-!--          justified since the top boundary is far away from the domain
-!--          of interest. 
-             bc_p_t   = 'dirichlet' !'neumann'
+             bc_p_t   = 'neumann'
           ENDIF
 
Index: palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90
===================================================================
--- palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90	(revision 4020)
+++ palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90	(revision 4022)
@@ -25,4 +25,17 @@
 ! -----------------
 ! $Id$
+! Several bugfixes and improvements
+! - revise bias correction of the imposed perturbations (correction via volume
+!   flow can create instabilities in case the mean volume flow is close to zero)
+! - introduce lower limits in calculation of coefficient matrix, else the 
+!   calculation may become numerically unstable
+! - impose perturbations every timestep, even though no new set of perturbations
+!   is generated in case dt_stg_call /= dt_3d
+! - Implement a gradual decrease of Reynolds stress and length scales above 
+!   ABL height (within 1 length scale above ABL depth to 1/10) rather than a 
+!   discontinuous decrease
+! - Bugfix in non-nested case: use ABL height for parametrized turbulence
+! 
+! 3987 2019-05-22 09:52:13Z kanani
 ! Introduce alternative switch for debug output during timestepping
 ! 
@@ -206,18 +219,68 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  mean_inflow_profiles, q, pt, u, v, w, zu, zw
+        ONLY:  dzw,                                                            &
+               ddzw,                                                           &
+               drho_air,                                                       &
+               mean_inflow_profiles,                                           &
+               q,                                                              &
+               pt,                                                             & 
+               u,                                                              &
+               u_init,                                                         &
+               v,                                                              &
+               v_init,                                                         &
+               w,                                                              & 
+               zu,                                                             &
+               zw
 
     USE basic_constants_and_equations_mod,                                     &
-        ONLY:  g, kappa, pi
+        ONLY:  g,                                                              &
+               kappa,                                                          &
+               pi
 
     USE control_parameters,                                                    &
-        ONLY:  debug_output_timestep,                                          &
+        ONLY:  bc_lr,                                                          &
+               bc_ns,                                                          &
+               child_domain,                                                   &
+               coupling_char,                                                  &
+               debug_output_timestep,                                          &
+               dt_3d,                                                          &
+               e_init,                                                         &
                initializing_actions,                                           &
+               intermediate_timestep_count,                                    &
+               intermediate_timestep_count_max,                                &
+               length,                                                         &
                message_string,                                                 &
+               nesting_offline,                                                &
                num_mean_inflow_profiles,                                       &
-               syn_turb_gen
+               rans_mode,                                                      &
+               restart_string,                                                 &
+               syn_turb_gen,                                                   &
+               time_since_reference_point,                                     &
+               turbulent_inflow
+               
+    USE grid_variables,                                                        &
+        ONLY:  ddx,                                                            &
+               ddy,                                                            & 
+               dx,                                                             &
+               dy
 
     USE indices,                                                               &
-        ONLY:  nbgp, nzb, nzt, nxl, nxlg, nxr, nxrg, nys, nyn, nyng, nysg
+        ONLY:  nbgp,                                                           & 
+               nz,                                                             & 
+               nzb,                                                            &
+               nzt,                                                            &
+               nx,                                                             & 
+               nxl,                                                            & 
+               nxlu,                                                           &
+               nxlg,                                                           & 
+               nxr,                                                            & 
+               nxrg,                                                           &
+               ny,                                                             &
+               nys,                                                            &
+               nysv,                                                           &
+               nyn,                                                            &
+               nyng,                                                           &
+               nysg,                                                           &
+               wall_flags_0
 
     USE kinds
@@ -228,11 +291,22 @@
 
     USE nesting_offl_mod,                                                      &
-        ONLY:  zi_ribulk
+        ONLY:  nesting_offl_calc_zi,                                           &
+               zi_ribulk
 
     USE pegrid,                                                                &
-        ONLY:  comm1dx, comm1dy, comm2d, ierr, myidx, myidy, pdims
+        ONLY:  comm1dx,                                                        &
+               comm1dy,                                                        &
+               comm2d,                                                         &
+               ierr,                                                           &
+               myidx,                                                          &
+               myidy,                                                          &
+               pdims
+        
+    USE pmc_interface,                                                         &
+        ONLY : rans_mode_parent
 
     USE transpose_indices,                                                     &
-        ONLY: nzb_x, nzt_x
+        ONLY: nzb_x,                                                           &
+              nzt_x
 
 
@@ -263,4 +337,9 @@
     INTEGER(iwp) ::  nzt_y_stg          !< upper bound of z coordinate (required for transposing z on PEs along y)
 
+    INTEGER(iwp), DIMENSION(3) ::  nr_non_topo_xz = 0 !< number of non-topography grid points at xz cross-sections,
+                                                      !< required for bias correction of imposed perturbations
+    INTEGER(iwp), DIMENSION(3) ::  nr_non_topo_yz = 0 !< number of non-topography grid points at yz cross-sections,
+                                                      !< required for bias correction of imposed perturbations
+    
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  displs_xz      !< displacement for MPI_GATHERV
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  recv_count_xz  !< receive count for MPI_GATHERV
@@ -278,7 +357,12 @@
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  seed           !< seed of random number for rn-generator
 
+    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 
+                                          !< one length scale above ABL top
+    REAL(wp) ::  d_l                      !< blend_coeff/length_scale
+    REAL(wp) ::  length_scale             !< length scale, default is 8 x minimum grid spacing
     REAL(wp) ::  dt_stg_adjust = 300.0_wp !< time interval for adjusting turbulence statistics
     REAL(wp) ::  dt_stg_call = 5.0_wp     !< time interval for calling synthetic turbulence generator
-    REAL(wp) ::  mc_factor = 1.0_wp       !< mass flux correction factor
     REAL(wp) ::  scale_l                  !< scaling parameter used for turbulence parametrization - Obukhov length
     REAL(wp) ::  scale_us                 !< scaling parameter used for turbulence parametrization - friction velocity
@@ -286,4 +370,6 @@
     REAL(wp) ::  time_stg_adjust = 0.0_wp !< time counter for adjusting turbulence information   
     REAL(wp) ::  time_stg_call = 0.0_wp   !< time counter for calling generator   
+    
+    REAL(wp), DIMENSION(3) ::  mc_factor = 1.0_wp !< correction factor for the u,v,w-components to maintain original mass flux
     
     
@@ -421,13 +507,4 @@
  SUBROUTINE stg_check_parameters
 
-
-    USE control_parameters,                                                    &
-        ONLY:  bc_lr, bc_ns, child_domain, nesting_offline, rans_mode,         &
-               turbulent_inflow
-
-    USE pmc_interface,                                                         &
-        ONLY : rans_mode_parent
-
-
     IMPLICIT NONE
 
@@ -548,21 +625,4 @@
  SUBROUTINE stg_init
 
-
-    USE arrays_3d,                                                             &
-        ONLY:  ddzw, u_init, v_init
-
-    USE control_parameters,                                                    &
-        ONLY:  child_domain, coupling_char, e_init, nesting_offline, rans_mode
-
-    USE grid_variables,                                                        &
-        ONLY:  ddy
-
-    USE indices,                                                               &
-        ONLY:  nz
-
-    USE pmc_interface,                                                         &
-        ONLY : rans_mode_parent
-
-
     IMPLICIT NONE
 
@@ -581,5 +641,7 @@
     INTEGER(iwp) :: nseed                    !< dimension of random number seed
     INTEGER(iwp) :: startseed = 1234567890   !< start seed for random number generator
-
+    
+    INTEGER(iwp), DIMENSION(3) ::  nr_non_topo_xz_l = 0 !< number of non-topography grid points at xz-cross-section on subdomain
+    INTEGER(iwp), DIMENSION(3) ::  nr_non_topo_yz_l = 0 !< number of non-topography grid points at yz-cross-section on subdomain
 !
 !-- Dummy variables used for reading profiles
@@ -782,6 +844,21 @@
     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 decrease 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 ) )
+!
 !--    Set flag indicating that turbulence is parametrized
        parametrize_inflow_turbulence = .TRUE.
+!
+!--    In case of dirichlet inflow boundary conditions only at one lateral 
+!--    boundary, i.e. in the case no offline or self nesting is applied but 
+!--    synthetic turbulence shall be parametrized nevertheless, the 
+!--    boundary-layer depth need to determined first.
+       IF ( .NOT. nesting_offline  .AND.  .NOT.  child_domain )                &
+          CALL nesting_offl_calc_zi
 !
 !--    Determine boundary-layer depth, which is used to initialize lenght 
@@ -897,5 +974,5 @@
              DO  k = nzb, nzt+1
                 IF  ( a11(k) .NE. 0._wp ) THEN
-                   fu_yz(k,j) = ( u(k,j,i) / mc_factor - u_init(k) ) / a11(k)
+                   fu_yz(k,j) = ( u(k,j,i) - u_init(k) ) / a11(k)
                 ELSE
                    fu_yz(k,j) = 0._wp
@@ -903,6 +980,6 @@
 
                 IF  ( a22(k) .NE. 0._wp ) THEN
-                   fv_yz(k,j) = ( v(k,j,i) / mc_factor - a21(k) * fu_yz(k,j) - &
-                               v_init(k) ) / a22(k)
+                   fv_yz(k,j) = ( v(k,j,i) -                                  &
+                                  a21(k) * fu_yz(k,j) - v_init(k) ) / a22(k)
                 ELSE
                    fv_yz(k,j) = 0._wp
@@ -910,6 +987,7 @@
 
                 IF  ( a33(k) .NE. 0._wp ) THEN
-                   fw_yz(k,j) = ( w(k,j,i) / mc_factor - a31(k) * fu_yz(k,j) - &
-                               a32(k) * fv_yz(k,j) ) / a33(k)
+                   fw_yz(k,j) = ( w(k,j,i) -                                   &
+                                  a31(k) * fu_yz(k,j) - a32(k) *               &
+                                  fv_yz(k,j) ) / a33(k)
                 ELSE
                    fw_yz = 0._wp
@@ -929,5 +1007,5 @@
 
                 IF  ( a11(k) .NE. 0._wp ) THEN
-                   fu_xz(k,i) = ( u(k,j,i) / mc_factor - u_init(k) ) / a11(k)
+                   fu_xz(k,i) = ( u(k,j,i) - u_init(k) ) / a11(k)
                 ELSE
                    fu_xz(k,i) = 0._wp
@@ -935,6 +1013,6 @@
 
                 IF  ( a22(k) .NE. 0._wp ) THEN
-                   fv_xz(k,i) = ( v(k,j,i) / mc_factor - a21(k) * fu_xz(k,i) - &
-                               v_init(k) ) / a22(k)
+                   fv_xz(k,i) = ( v(k,j,i) -                                   &
+                                  a21(k) * fu_xz(k,i) - v_init(k) ) / a22(k)
                 ELSE
                    fv_xz(k,i) = 0._wp
@@ -942,6 +1020,7 @@
 
                 IF  ( a33(k) .NE. 0._wp ) THEN
-                   fw_xz(k,i) = ( w(k,j,i) / mc_factor - a31(k) * fu_xz(k,i) - &
-                               a32(k) * fv_xz(k,i) ) / a33(k)
+                   fw_xz(k,i) = ( w(k,j,i) -                                   &
+                                  a31(k) * fu_xz(k,i) -                        &
+                                  a32(k) * fv_xz(k,i) ) / a33(k)
                 ELSE
                    fw_xz = 0._wp
@@ -952,4 +1031,68 @@
        ENDIF
     ENDIF
+!
+!-- Count the number of non-topography grid points at the boundaries where 
+!-- perturbations are imposed. This number is later used for bias corrections
+!-- of the perturbations, i.e. to force that their mean is zero. Please note,
+!-- due to the asymetry of u and v along x and y direction, respectively, 
+!-- different cases must be distinguished. 
+    IF ( myidx == id_stg_left  .OR.  myidx == id_stg_right )  THEN
+!
+!--    Number of grid points where perturbations are imposed on u
+       IF ( myidx == id_stg_left  )  i = nxl
+       IF ( myidx == id_stg_right )  i = nxr+1
+       
+       nr_non_topo_yz_l(1) = SUM(                                              &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzt,nys:nyn,i), 1 ) ) )
+!
+!--    Number of grid points where perturbations are imposed on v and w                                   
+       IF ( myidx == id_stg_left  )  i = nxl-1
+       IF ( myidx == id_stg_right )  i = nxr+1
+       
+       nr_non_topo_yz_l(2) = SUM(                                              &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzt,nysv:nyn,i), 2 ) ) )
+       nr_non_topo_yz_l(3) = SUM(                                              &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzt,nys:nyn,i), 3 ) ) )
+                                    
+#if defined( __parallel )
+       CALL MPI_ALLREDUCE( nr_non_topo_yz_l, nr_non_topo_yz, 3, MPI_INTEGER,   &
+                           MPI_SUM, comm1dy, ierr )
+#else
+       nr_non_topo_yz = nr_non_topo_yz_l
+#endif  
+    ENDIF
+    
+    IF ( myidy == id_stg_south  .OR.  myidy == id_stg_north )  THEN
+!
+!--    Number of grid points where perturbations are imposed on v
+       IF ( myidy == id_stg_south )  j = nys
+       IF ( myidy == id_stg_north )  j = nyn+1
+       
+       nr_non_topo_xz_l(2) = SUM(                                              &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzt,j,nxl:nxr), 2 ) ) )
+!
+!--    Number of grid points where perturbations are imposed on u and w 
+       IF ( myidy == id_stg_south )  j = nys-1
+       IF ( myidy == id_stg_north )  j = nyn+1
+       
+       nr_non_topo_xz_l(1) = SUM(                                              &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzt,j,nxlu:nxr), 1 ) ) )
+       nr_non_topo_xz_l(3) = SUM(                                              &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzt,j,nxl:nxr), 3 ) ) )
+                                    
+#if defined( __parallel )
+       CALL MPI_ALLREDUCE( nr_non_topo_xz_l, nr_non_topo_xz, 3, MPI_INTEGER,   &
+                           MPI_SUM, comm1dx, ierr )
+#else
+       nr_non_topo_xz = nr_non_topo_xz_l
+#endif  
+    ENDIF
+
 
  END SUBROUTINE stg_init
@@ -1059,8 +1202,4 @@
 
 
-    USE control_parameters,                                                    &
-        ONLY: length, restart_string
-
-
     IMPLICIT NONE
 
@@ -1072,7 +1211,11 @@
 
     SELECT CASE ( restart_string(1:length) )
-
-       CASE ( 'mc_factor' )
-          READ ( 13 )  mc_factor
+          
+       CASE ( 'time_stg_adjust' )
+          READ ( 13 )  time_stg_adjust
+          
+       CASE ( 'time_stg_call' )
+          READ ( 13 )  time_stg_call
+          
        CASE ( 'use_syn_turb_gen' )
           READ ( 13 )  use_syn_turb_gen
@@ -1097,7 +1240,10 @@
 
     IMPLICIT NONE
-
-    CALL wrd_write_string( 'mc_factor' )
-    WRITE ( 14 )  mc_factor
+    
+    CALL wrd_write_string( 'time_stg_adjust' )
+    WRITE ( 14 )  time_stg_adjust
+    
+    CALL wrd_write_string( 'time_stg_call' )
+    WRITE ( 14 )  time_stg_call
 
     CALL wrd_write_string( 'use_syn_turb_gen' )
@@ -1118,23 +1264,4 @@
  SUBROUTINE stg_main
 
-
-    USE arrays_3d,                                                             &
-        ONLY:  dzw
-
-    USE control_parameters,                                                    &
-        ONLY:  child_domain, dt_3d,                                            &
-               nesting_offline, rans_mode, time_since_reference_point,         &
-               volume_flow_initial
-
-    USE grid_variables,                                                        &
-        ONLY:  dx, dy
-
-    USE indices,                                                               &
-        ONLY:  wall_flags_0
-
-    USE pmc_interface,                                                         &
-        ONLY : rans_mode_parent
-
-
     IMPLICIT NONE
 
@@ -1142,10 +1269,10 @@
     INTEGER(iwp) :: j           !< loop index in y-direction
     INTEGER(iwp) :: k           !< loop index in z-direction
+    
+    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) :: mc_factor_l     !< local mass flux correction factor
-    REAL(wp) :: volume_flow     !< mass flux through lateral boundary
-    REAL(wp) :: volume_flow_l   !< local mass flux through lateral boundary
-
+    
+    REAL(wp), DIMENSION(3) :: mc_factor_l   !< local mass flux correction factor
 
     IF ( debug_output_timestep )  CALL debug_message( 'stg_main', 'start' )
@@ -1153,4 +1280,14 @@
 !-- Calculate time step which is needed for filter functions
     dt_stg = MAX( dt_3d, dt_stg_call )
+!
+!-- Check if synthetic turbulence generator needs to be called and new 
+!-- perturbations need to be created or if old disturbances can be imposed
+!-- again. 
+    IF ( time_stg_call >= dt_stg_call  .AND.                                  &
+         intermediate_timestep_count == intermediate_timestep_count_max  )  THEN
+       stg_call = .TRUE.
+    ELSE
+       stg_call = .FALSE.
+    ENDIF
 !
 !-- Initial value of fu, fv, fw
@@ -1183,25 +1320,30 @@
 !
 !-- New set of fu, fv, fw
-    CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_left )
-    CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_left )
-    CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_left )
-    
-    IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent         &
-                          .AND.  .NOT.  rans_mode ) )  THEN
-!
-!--    Generate turbulence at right boundary
-       CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_right )
-       CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_right )
-       CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_right )
-!
-!--    Generate turbulence at north boundary
-       CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_north )
-       CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_north )
-       CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_north )
-!
-!--    Generate turbulence at south boundary
-       CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_south )
-       CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_south )
-       CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_south )
+    IF ( stg_call )  THEN
+!      
+!--    Generate turbulence at left boundary (required in all applications, 
+!--    i.e. offline nesting and turbulent inflow from the left)
+       CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_left )
+       CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_left )
+       CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_left )
+       
+       IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent      &
+                             .AND.  .NOT.  rans_mode ) )  THEN
+!      
+!--       Generate turbulence at right boundary
+          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_right )
+          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_right )
+          CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_right )
+!      
+!--       Generate turbulence at north boundary
+          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_north )
+          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_north )
+          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_north )
+!      
+!--       Generate turbulence at south boundary
+          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_south )
+          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_south )
+          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_south )
+       ENDIF
     ENDIF
     
@@ -1209,139 +1351,180 @@
 !-- Turbulence generation at left and or right boundary
     IF ( myidx == id_stg_left  .OR.  myidx == id_stg_right )  THEN
-    
-       DO  j = nysg, nyng
-          DO  k = nzb, nzt + 1
-!
-!--          Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
-             IF ( tu(k) == 0.0_wp )  THEN
-                fu_yz(k,j) = fuo_yz(k,j)
-             ELSE
-                fu_yz(k,j) = fu_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) + &
-                         fuo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
-             ENDIF
-
-             IF ( tv(k) == 0.0_wp )  THEN
-                fv_yz(k,j) = fvo_yz(k,j)
-             ELSE
-                fv_yz(k,j) = fv_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) + &
-                         fvo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
-             ENDIF
-
-             IF ( tw(k) == 0.0_wp )  THEN
-                fw_yz(k,j) = fwo_yz(k,j)
-             ELSE
-                fw_yz(k,j) = fw_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) + &
-                         fwo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
-             ENDIF
-!
-!--          Lund rotation following Eq. 17 in Xie and Castro (2008).
-!--          Additional factors are added to improve the variance of v and w
-             IF( k == 0 )  THEN
-                dist_yz(k,j,1) = 0.0_wp
-                dist_yz(k,j,2) = 0.0_wp
-                dist_yz(k,j,3) = 0.0_wp
-             ELSE
-                dist_yz(k,j,1) = MIN( a11(k) * fu_yz(k,j), 3.0_wp )
-                !experimental test of 1.2
+!
+!--    Calculate new set of perturbations. Do this only at last RK3-substep and
+!--    when dt_stg_call is exceeded. Else the old set of perturbations is 
+!--    imposed
+       IF ( stg_call )  THEN
+       
+          DO  j = nysg, nyng
+             DO  k = nzb, nzt + 1
+!      
+!--             Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
+                IF ( tu(k) == 0.0_wp )  THEN
+                   fu_yz(k,j) = fuo_yz(k,j)
+                ELSE
+                   fu_yz(k,j) = fu_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) + &
+                            fuo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
+                ENDIF
+       
+                IF ( tv(k) == 0.0_wp )  THEN
+                   fv_yz(k,j) = fvo_yz(k,j)
+                ELSE
+                   fv_yz(k,j) = fv_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) + &
+                            fvo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
+                ENDIF
+       
+                IF ( tw(k) == 0.0_wp )  THEN
+                   fw_yz(k,j) = fwo_yz(k,j)
+                ELSE
+                   fw_yz(k,j) = fw_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) + &
+                            fwo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
+                ENDIF
+             ENDDO
+          ENDDO
+          
+          dist_yz(nzb,:,1) = 0.0_wp
+          dist_yz(nzb,:,2) = 0.0_wp
+          dist_yz(nzb,:,3) = 0.0_wp
+          
+          IF ( myidx == id_stg_left  )  i = nxl
+          IF ( myidx == id_stg_right )  i = nxr+1       
+          DO  j = nysg, nyng
+             DO  k = nzb+1, nzt + 1
+!      
+!--             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 ) *          &
+                                    MERGE( 1.0_wp, 0.0_wp,                     &
+                                           BTEST( wall_flags_0(k,j,i), 1 ) )   
+             ENDDO
+          ENDDO
+          
+          IF ( myidx == id_stg_left  )  i = nxl-1
+          IF ( myidx == id_stg_right )  i = nxr+1
+          DO  j = nysg, nyng
+             DO  k = nzb+1, nzt + 1
+!      
+!--             Lund rotation following Eq. 17 in Xie and Castro (2008).
+!--             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 )
+                                         + a22(k) * fv_yz(k,j) ), 3.0_wp ) *   &
+                                    MERGE( 1.0_wp, 0.0_wp,                     &
+                                           BTEST( wall_flags_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 )
-             ENDIF
-
+                                        + a33(k) * fw_yz(k,j) ), 3.0_wp )  *   &
+                                    MERGE( 1.0_wp, 0.0_wp,                     &
+                                           BTEST( wall_flags_0(k,j,i), 3 ) )
+             ENDDO
           ENDDO
-       ENDDO
+       ENDIF
 
 !
 !--    Mass flux correction following Kim et al. (2013)
 !--    This correction factor insures that the mass flux is preserved at the
-!--    inflow boundary
-       IF ( .NOT. nesting_offline  .AND.  .NOT. child_domain )  THEN 
-          mc_factor_l = 0.0_wp
-          mc_factor   = 0.0_wp
-          DO  j = nys, nyn
-             DO  k = nzb+1, nzt
-                mc_factor_l = mc_factor_l + dzw(k)  *                          &
-                              ( mean_inflow_profiles(k,1) + dist_yz(k,j,1) )
-             ENDDO
-          ENDDO
-
-#if defined( __parallel )
-          CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, 1, MPI_REAL, MPI_SUM,    &
-                              comm1dy, ierr )
-#else
-          mc_factor = mc_factor_l
-#endif
-
-          mc_factor = volume_flow_initial(1) / mc_factor
-
-!
-!--       Add disturbance at the inflow
-          DO  j = nysg, nyng
-             DO  k = nzb, nzt+1
-                 u(k,j,-nbgp+1:0) = ( mean_inflow_profiles(k,1) +              &
-                                      dist_yz(k,j,1)             ) * mc_factor
-                 v(k,j,-nbgp:-1)  = ( mean_inflow_profiles(k,2) +              &
-                                      dist_yz(k,j,2)             ) * mc_factor
-                 w(k,j,-nbgp:-1)  =   dist_yz(k,j,3)               * mc_factor
-             ENDDO
-          ENDDO
-
-       ELSE
-!
-!--       First, calculate volume flow at yz boundary
+!--    inflow boundary. First, calculate mean value of the imposed 
+!--    perturbations at yz boundary. 
+!--    Note, this needs to be done only after the last RK3-substep, else 
+!--    the boundary values won't be accessed. 
+       IF ( intermediate_timestep_count == intermediate_timestep_count_max )  THEN
+          mc_factor_l   = 0.0_wp
+          mc_factor     = 0.0_wp
+!       
+!--       Sum up the original volume flows (with and without perturbations). 
+!--       Note, for non-normal components (here v and w) it is actually no 
+!--       volume flow. 
           IF ( myidx == id_stg_left  )  i = nxl
           IF ( myidx == id_stg_right )  i = nxr+1
-
-          volume_flow_l = 0.0_wp
-          volume_flow   = 0.0_wp
-          mc_factor_l   = 0.0_wp
-          mc_factor     = 0.0_wp
-          DO  j = nys, nyn
-             DO  k = nzb+1, nzt
-                volume_flow_l = volume_flow_l + u(k,j,i) * dzw(k) * dy         &
-                                     * MERGE( 1.0_wp, 0.0_wp,                  &
-                                              BTEST( wall_flags_0(k,j,i), 1 ) )
-
-                mc_factor_l = mc_factor_l     + ( u(k,j,i) + dist_yz(k,j,1) )  &
-                                                         * dzw(k) * dy         &
-                                     * MERGE( 1.0_wp, 0.0_wp,                  &
-                                              BTEST( wall_flags_0(k,j,i), 1 ) )
-             ENDDO
-          ENDDO
+         
+          mc_factor_l(1) = SUM( dist_yz(nzb:nzt,nys:nyn,1)                     &
+                         * MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(nzb:nzt,nys:nyn,i), 1 ) ) )
+        
+          IF ( myidx == id_stg_left  )  i = nxl-1
+          IF ( myidx == id_stg_right )  i = nxr+1
+          
+          mc_factor_l(2) = SUM( dist_yz(nzb:nzt,nysv:nyn,2)                    &
+                          * MERGE( 1.0_wp, 0.0_wp,                             &
+                                   BTEST( wall_flags_0(nzb:nzt,nysv:nyn,i), 2 ) ) )
+          mc_factor_l(3) = SUM( dist_yz(nzb:nzt,nys:nyn,3)                     &
+                          * MERGE( 1.0_wp, 0.0_wp,                             &
+                                   BTEST( wall_flags_0(nzb:nzt,nys:nyn,i), 3 ) ) )
+          
 #if defined( __parallel )
-          CALL MPI_ALLREDUCE( volume_flow_l, volume_flow,                      &
-                              1, MPI_REAL, MPI_SUM, comm1dy, ierr )
           CALL MPI_ALLREDUCE( mc_factor_l, mc_factor,                          &
-                              1, MPI_REAL, MPI_SUM, comm1dy, ierr )
-#else
-          volume_flow = volume_flow_l
-          mc_factor   = mc_factor_l
+                              3, MPI_REAL, MPI_SUM, comm1dy, ierr )            
+#else                                                                          
+          mc_factor   = mc_factor_l                                            
 #endif
-
-          mc_factor = volume_flow / mc_factor
-
-!
+!      
+!--       Calculate correction factor and force zero mean perturbations. 
+          mc_factor = mc_factor / REAL( nr_non_topo_yz, KIND = wp )            
+                                                                               
+          IF ( myidx == id_stg_left  )  i = nxl                                
+          IF ( myidx == id_stg_right )  i = nxr+1                              
+                                                                               
+          dist_yz(:,:,1) = ( dist_yz(:,:,1) - mc_factor(1) )                   &
+                        * MERGE( 1.0_wp, 0.0_wp,                               &
+                                 BTEST( wall_flags_0(:,:,i), 1 ) )             
+                                                                               
+                                                                               
+          IF ( myidx == id_stg_left  )  i = nxl-1                              
+          IF ( myidx == id_stg_right )  i = nxr+1                              
+                                                                               
+          dist_yz(:,:,2) = ( dist_yz(:,:,2) - mc_factor(2) )                   &
+                        * MERGE( 1.0_wp, 0.0_wp,                               &
+                                 BTEST( wall_flags_0(:,:,i), 2 ) )             
+                                                                               
+          dist_yz(:,:,3) = ( dist_yz(:,:,3) - mc_factor(3) )                   &
+                        * MERGE( 1.0_wp, 0.0_wp,                               &
+                                 BTEST( wall_flags_0(:,:,i), 3 ) )
+!      
 !--       Add disturbances
           IF ( myidx == id_stg_left  )  THEN
-             DO  j = nys, nyn
-                DO  k = nzb+1, nzt
-                   u(k,j,0) = ( u(k,j,0) + dist_yz(k,j,1) )                    &
-                       * mc_factor * MERGE( 1.0_wp, 0.0_wp,                    &
-                                            BTEST( wall_flags_0(k,j,0), 1 ) )
-                   u(k,j,-1) = u(k,j,0)
-                   v(k,j,-1)  = ( v(k,j,-1)  + dist_yz(k,j,2)  )               &
-                       * mc_factor * MERGE( 1.0_wp, 0.0_wp,                    &
-                                            BTEST( wall_flags_0(k,j,-1), 2 ) )
-                   w(k,j,-1)  = ( w(k,j,-1)  + dist_yz(k,j,3)  )               &
-                       * mc_factor * MERGE( 1.0_wp, 0.0_wp,                    &
+!      
+!--          For the left boundary distinguish between mesoscale offline / self
+!--          nesting and turbulent inflow at the left boundary only. In the latter
+!--          case turbulence is imposed on the mean inflow profiles. 
+             IF ( .NOT. nesting_offline  .AND.  .NOT. child_domain )  THEN 
+!            
+!--             Add disturbance at the inflow
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      u(k,j,-nbgp+1:0) = ( mean_inflow_profiles(k,1) +         &
+                                           dist_yz(k,j,1)             )        &
+                                   * MERGE( 1.0_wp, 0.0_wp,                    &
+                                            BTEST( wall_flags_0(k,j,0), 1 ) )  
+                      v(k,j,-nbgp:-1)  = ( mean_inflow_profiles(k,2) +         &
+                                           dist_yz(k,j,2)             )        &
+                                   * MERGE( 1.0_wp, 0.0_wp,                    &
+                                            BTEST( wall_flags_0(k,j,-1), 2 ) ) 
+                      w(k,j,-nbgp:-1)  =   dist_yz(k,j,3)                      &
+                                   * MERGE( 1.0_wp, 0.0_wp,                    &
                                             BTEST( wall_flags_0(k,j,-1), 3 ) )
+                   ENDDO
                 ENDDO
-             ENDDO
+             ELSE
+             
+                DO  j = nys, nyn
+                   DO  k = nzb+1, nzt
+                      u(k,j,0)   = ( u(k,j,0) + dist_yz(k,j,1) )               &
+                                 * MERGE( 1.0_wp, 0.0_wp,                      &
+                                          BTEST( wall_flags_0(k,j,0), 1 ) )
+                      u(k,j,-1)  = u(k,j,0)
+                      v(k,j,-1)  = ( v(k,j,-1)  + dist_yz(k,j,2)  )            &
+                                 * MERGE( 1.0_wp, 0.0_wp,                      &
+                                          BTEST( wall_flags_0(k,j,-1), 2 ) )
+                      w(k,j,-1)  = ( w(k,j,-1)  + dist_yz(k,j,3) )             &
+                                 * MERGE( 1.0_wp, 0.0_wp,                      &
+                                          BTEST( wall_flags_0(k,j,-1), 3 ) )
+                   ENDDO
+                ENDDO
+             ENDIF
           ENDIF
           IF ( myidx == id_stg_right  )  THEN
@@ -1349,146 +1532,190 @@
                 DO  k = nzb+1, nzt
                    u(k,j,nxr+1) = ( u(k,j,nxr+1) + dist_yz(k,j,1) )            &
-                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
-                                          BTEST( wall_flags_0(k,j,nxr+1), 1 ) )
+                                  * MERGE( 1.0_wp, 0.0_wp,                     &
+                                           BTEST( wall_flags_0(k,j,nxr+1), 1 ) ) 
                    v(k,j,nxr+1) = ( v(k,j,nxr+1) + dist_yz(k,j,2) )            &
-                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
-                                          BTEST( wall_flags_0(k,j,nxr+1), 2 ) )
+                                  * MERGE( 1.0_wp, 0.0_wp,                     &
+                                           BTEST( wall_flags_0(k,j,nxr+1), 2 ) )
                    w(k,j,nxr+1) = ( w(k,j,nxr+1) + dist_yz(k,j,3) )            &
-                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
-                                          BTEST( wall_flags_0(k,j,nxr+1), 3 ) )
+                                  * MERGE( 1.0_wp, 0.0_wp,                     &
+                                           BTEST( wall_flags_0(k,j,nxr+1), 3 ) )
                 ENDDO
              ENDDO
           ENDIF
        ENDIF
-
     ENDIF
 !
 !-- Turbulence generation at north and south boundary
     IF ( myidy == id_stg_north  .OR.  myidy == id_stg_south )  THEN
-   
-       DO  i = nxlg, nxrg
-          DO  k = nzb, nzt + 1
-!
-!--          Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
-             IF ( tu(k) == 0.0_wp )  THEN
-                fu_xz(k,i) = fuo_xz(k,i)
-             ELSE
-                fu_xz(k,i) = fu_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) +     &
-                         fuo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
-             ENDIF
-
-             IF ( tv(k) == 0.0_wp )  THEN
-                fv_xz(k,i) = fvo_xz(k,i)
-             ELSE
-                fv_xz(k,i) = fv_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) +     &
-                      fvo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
-             ENDIF
-            
-             IF ( tw(k) == 0.0_wp )  THEN
-                fw_xz(k,i) = fwo_xz(k,i)
-             ELSE
-                fw_xz(k,i) = fw_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) +     &
-                      fwo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
-             ENDIF
-!
-!--          Lund rotation following Eq. 17 in Xie and Castro (2008).
-!--          Additional factors are added to improve the variance of v and w
-             IF( k == 0 )  THEN
-                dist_xz(k,i,1) = 0.0_wp
-                dist_xz(k,i,2) = 0.0_wp
-                dist_xz(k,i,3) = 0.0_wp
-
-             ELSE
-                dist_xz(k,i,1) = MIN( a11(k) * fu_xz(k,i), 3.0_wp )
-                !experimental test of 1.2
+!
+!--    Calculate new set of perturbations. Do this only at last RK3-substep and
+!--    when dt_stg_call is exceeded. Else the old set of perturbations is 
+!--    imposed
+       IF ( stg_call )  THEN
+          DO  i = nxlg, nxrg
+             DO  k = nzb, nzt + 1
+!         
+!--             Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
+                IF ( tu(k) == 0.0_wp )  THEN
+                   fu_xz(k,i) = fuo_xz(k,i)
+                ELSE
+                   fu_xz(k,i) = fu_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) +     &
+                            fuo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
+                ENDIF
+          
+                IF ( tv(k) == 0.0_wp )  THEN
+                   fv_xz(k,i) = fvo_xz(k,i)
+                ELSE
+                   fv_xz(k,i) = fv_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) +     &
+                         fvo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
+                ENDIF
+               
+                IF ( tw(k) == 0.0_wp )  THEN
+                   fw_xz(k,i) = fwo_xz(k,i)
+                ELSE
+                   fw_xz(k,i) = fw_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) +     &
+                         fwo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
+                ENDIF
+             ENDDO
+          ENDDO
+          
+          
+          dist_xz(nzb,:,1) = 0.0_wp
+          dist_xz(nzb,:,2) = 0.0_wp
+          dist_xz(nzb,:,3) = 0.0_wp
+          
+          IF ( myidy == id_stg_south  ) j = nys
+          IF ( myidy == id_stg_north )  j = nyn+1
+          DO  i = nxlg, nxrg
+             DO  k = nzb+1, nzt + 1
+!         
+!--             Lund rotation following Eq. 17 in Xie and Castro (2008).
+!--             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 )
+                                         + a22(k) * fv_xz(k,i) ), 3.0_wp ) *   &
+                                    MERGE( 1.0_wp, 0.0_wp,                     &
+                                           BTEST( wall_flags_0(k,j,i), 2 ) )   
+             ENDDO
+          ENDDO
+          
+          IF ( myidy == id_stg_south  ) j = nys-1
+          IF ( myidy == id_stg_north )  j = nyn+1
+          DO  i = nxlg, nxrg
+             DO  k = nzb+1, nzt + 1
+!         
+!--             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 ) *          &
+                                    MERGE( 1.0_wp, 0.0_wp,                     &
+                                           BTEST( wall_flags_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 )
-             ENDIF
-
-          ENDDO
-       ENDDO
-
-!
-!--    Mass flux correction following Kim et al. (2013)
-!--    This correction factor insures that the mass flux is preserved at the
-!--    inflow boundary.
-!--    First, calculate volume flow at xz boundary
-       IF ( myidy == id_stg_south  ) j = nys
-       IF ( myidy == id_stg_north )  j = nyn+1
-
-       volume_flow_l = 0.0_wp
-       volume_flow   = 0.0_wp
-       mc_factor_l   = 0.0_wp
-       mc_factor     = 0.0_wp
-       DO  i = nxl, nxr
-          DO  k = nzb+1, nzt
-             volume_flow_l = volume_flow_l + v(k,j,i) * dzw(k) * dx            &
-                                  * MERGE( 1.0_wp, 0.0_wp,                     &
-                                           BTEST( wall_flags_0(k,j,i), 2 ) )
-
-             mc_factor_l = mc_factor_l     + ( v(k,j,i) + dist_xz(k,i,2) )     &
-                                                      * dzw(k) * dx            &
-                                  * MERGE( 1.0_wp, 0.0_wp,                     &
-                                           BTEST( wall_flags_0(k,j,i), 2 ) )
-          ENDDO
-       ENDDO
-#if defined( __parallel )
-       CALL MPI_ALLREDUCE( volume_flow_l, volume_flow,                         &
-                           1, MPI_REAL, MPI_SUM, comm1dx, ierr )
-       CALL MPI_ALLREDUCE( mc_factor_l, mc_factor,                             &
-                           1, MPI_REAL, MPI_SUM, comm1dx, ierr )
-#else
-       volume_flow = volume_flow_l
-       mc_factor   = mc_factor_l
-#endif
-
-       mc_factor = volume_flow / mc_factor
-
-!
-!--    Add disturbances
-       IF ( myidy == id_stg_south  )  THEN
-
-          DO  i = nxl, nxr
-             DO  k = nzb+1, nzt
-                u(k,-1,i) = ( u(k,-1,i) + dist_xz(k,i,1) )                     &
-                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                     &
-                                           BTEST( wall_flags_0(k,-1,i), 1 ) )
-                v(k,0,i)  = ( v(k,0,i)  + dist_xz(k,i,2)  )                    &
-                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                     &
-                                           BTEST( wall_flags_0(k,0,i), 2 ) )
-                v(k,-1,i) = v(k,0,i)
-                w(k,-1,i) = ( w(k,-1,i) + dist_xz(k,i,3)  )                    &
-                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                     &
-                                           BTEST( wall_flags_0(k,-1,i), 3 ) )
+                                        + a33(k) * fw_xz(k,i) ), 3.0_wp )  *   &
+                                    MERGE( 1.0_wp, 0.0_wp,                     &
+                                           BTEST( wall_flags_0(k,j,i), 3 ) ) 
              ENDDO
           ENDDO
        ENDIF
-       IF ( myidy == id_stg_north  )  THEN
-
-          DO  i = nxl, nxr
-             DO  k = nzb+1, nzt
-                u(k,nyn+1,i) = ( u(k,nyn+1,i) + dist_xz(k,i,1) )               &
-                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
-                                          BTEST( wall_flags_0(k,nyn+1,i), 1 ) )
-                v(k,nyn+1,i) = ( v(k,nyn+1,i) + dist_xz(k,i,2) )               &
-                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
-                                          BTEST( wall_flags_0(k,nyn+1,i), 2 ) )
-                w(k,nyn+1,i) = ( w(k,nyn+1,i) + dist_xz(k,i,3) )               &
-                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
-                                          BTEST( wall_flags_0(k,nyn+1,i), 3 ) )
+
+!
+!--    Mass flux correction following Kim et al. (2013)
+!--    This correction factor insures that the mass flux is preserved at the
+!--    inflow boundary. First, calculate mean value of the imposed 
+!--    perturbations at yz boundary. 
+!--    Note, this needs to be done only after the last RK3-substep, else 
+!--    the boundary values won't be accessed. 
+       IF ( intermediate_timestep_count == intermediate_timestep_count_max )  THEN
+          mc_factor_l   = 0.0_wp
+          mc_factor     = 0.0_wp
+         
+          IF ( myidy == id_stg_south  ) j = nys
+          IF ( myidy == id_stg_north )  j = nyn+1
+          
+          mc_factor_l(2) = SUM( dist_xz(nzb:nzt,nxl:nxr,2)                     &
+                         * MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(nzb:nzt,j,nxl:nxr), 2 ) ) )
+          
+          IF ( myidy == id_stg_south  ) j = nys-1
+          IF ( myidy == id_stg_north )  j = nyn+1
+          
+          mc_factor_l(1) = SUM( dist_xz(nzb:nzt,nxlu:nxr,1)                    &
+                         * MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(nzb:nzt,j,nxlu:nxr), 1 ) ) )
+          mc_factor_l(3) = SUM( dist_xz(nzb:nzt,nxl:nxr,3)                     &
+                         * MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(nzb:nzt,j,nxl:nxr), 3 ) ) )
+          
+#if defined( __parallel )
+          CALL MPI_ALLREDUCE( mc_factor_l, mc_factor,                          &
+                              3, MPI_REAL, MPI_SUM, comm1dx, ierr )
+#else     
+          mc_factor   = mc_factor_l
+#endif
+          
+          mc_factor = mc_factor / REAL( nr_non_topo_xz, KIND = wp )
+          
+          IF ( myidy == id_stg_south  ) j = nys
+          IF ( myidy == id_stg_north )  j = nyn+1
+          
+          dist_xz(:,:,2)   = ( dist_xz(:,:,2) - mc_factor(2) )                 &
+                           * MERGE( 1.0_wp, 0.0_wp,                            &
+                                    BTEST( wall_flags_0(:,j,:), 2 ) )          
+                                                                               
+                                                                               
+          IF ( myidy == id_stg_south  ) j = nys-1                              
+          IF ( myidy == id_stg_north )  j = nyn+1                              
+                                                                               
+          dist_xz(:,:,1)   = ( dist_xz(:,:,1) - mc_factor(1) )                 &
+                           * MERGE( 1.0_wp, 0.0_wp,                            &
+                                    BTEST( wall_flags_0(:,j,:), 1 ) )          
+                                                                               
+          dist_xz(:,:,3)   = ( dist_xz(:,:,3) - mc_factor(3) )                 &
+                           * MERGE( 1.0_wp, 0.0_wp,                            &
+                                    BTEST( wall_flags_0(:,j,:), 3 ) )      
+!         
+!--       Add disturbances
+          IF ( myidy == id_stg_south  )  THEN
+             DO  i = nxl, nxr
+                DO  k = nzb+1, nzt
+                   u(k,-1,i) = ( u(k,-1,i) + dist_xz(k,i,1) )                  &
+                                        * MERGE( 1.0_wp, 0.0_wp,               &
+                                                 BTEST( wall_flags_0(k,-1,i), 1 ) )
+                   v(k,0,i)  = ( v(k,0,i)  + dist_xz(k,i,2)  )                 &
+                                        * MERGE( 1.0_wp, 0.0_wp,               &
+                                                 BTEST( wall_flags_0(k,0,i), 2 ) )
+                   v(k,-1,i) = v(k,0,i)
+                   w(k,-1,i) = ( w(k,-1,i) + dist_xz(k,i,3)  )                 &
+                                        * MERGE( 1.0_wp, 0.0_wp,               &
+                                                 BTEST( wall_flags_0(k,-1,i), 3 ) )
+                ENDDO
              ENDDO
-          ENDDO
+          ENDIF
+          IF ( myidy == id_stg_north  )  THEN
+          
+             DO  i = nxl, nxr
+                DO  k = nzb+1, nzt
+                   u(k,nyn+1,i) = ( u(k,nyn+1,i) + dist_xz(k,i,1) )            &
+                                     * MERGE( 1.0_wp, 0.0_wp,                  &
+                                             BTEST( wall_flags_0(k,nyn+1,i), 1 ) )
+                   v(k,nyn+1,i) = ( v(k,nyn+1,i) + dist_xz(k,i,2) )            &
+                                     * MERGE( 1.0_wp, 0.0_wp,                  &
+                                              BTEST( wall_flags_0(k,nyn+1,i), 2 ) )
+                   w(k,nyn+1,i) = ( w(k,nyn+1,i) + dist_xz(k,i,3) )            &
+                                     * MERGE( 1.0_wp, 0.0_wp,                  &
+                                              BTEST( wall_flags_0(k,nyn+1,i), 3 ) )
+                ENDDO
+             ENDDO
+          ENDIF
        ENDIF
     ENDIF
 !
 !-- Finally, set time counter for calling STG to zero
-    time_stg_call = 0.0_wp
+    IF ( stg_call )  time_stg_call = 0.0_wp
 
     IF ( debug_output_timestep )  CALL debug_message( 'stg_main', 'end' )
@@ -1506,8 +1733,4 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE stg_generate_seed_yz( n_y, n_z, b_y, b_z, f_n, id )
-
-
-    USE indices,                                                               &
-        ONLY: ny
 
     IMPLICIT NONE
@@ -1633,9 +1856,4 @@
  SUBROUTINE stg_generate_seed_xz( n_x, n_z, b_x, b_z, f_n, id )
 
-
-    USE indices,                                                               &
-        ONLY: nx
-
-
     IMPLICIT NONE
 
@@ -1758,7 +1976,4 @@
  SUBROUTINE parametrize_reynolds_stress
 
-    USE arrays_3d,                                                             &
-        ONLY:  zu
-
     IMPLICIT NONE
 
@@ -1767,53 +1982,41 @@
     REAL(wp)     ::  zzi         !< ratio of z/zi
     
-!
-!-- 
     DO  k = nzb+1, nzt+1
-          
-       IF ( zu(k) <= zi_ribulk )  THEN
-!
-!--       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 merge-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                                    &
-                                 )
-                                 
-          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            &                     
-                                 )
-!
-!--       u'w' and v'w'. Assume isotropy.
-          r31(k) = - scale_us**2 * (                                           &
-                         1.0_wp - EXP( 3.0_wp * ( zzi - 1.0_wp ) )             &
-                                   )
-                              
-          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)
-!
-!--    Above the boundary layer, assmume laminar flow conditions.
-       ELSE
-          r11(k) = 10E-8_wp
-          r22(k) = 10E-8_wp
-          r33(k) = 10E-8_wp
-          r21(k) = 10E-8_wp
-          r31(k) = 10E-8_wp
-          r32(k) = 10E-8_wp
-       ENDIF
+!
+!--    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 merge-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)
     ENDDO
 
@@ -1845,5 +2048,5 @@
 !-- Calculate coefficient matrix. Split loops to allow for loop vectorization.
     DO  k = nzb+1, nzt+1
-       IF ( r11(k) > 0.0_wp )  THEN
+       IF ( r11(k) > 10E-6_wp )  THEN
           a11(k) = SQRT( r11(k) ) 
           a21(k) = r21(k) / a11(k)
@@ -1857,5 +2060,5 @@
     DO  k = nzb+1, nzt+1
        a22(k) = r22(k) - a21(k)**2 
-       IF ( a22(k) > 0.0_wp )  THEN
+       IF ( a22(k) > 10E-6_wp )  THEN
           a22(k) = SQRT( a22(k) )
           a32(k) = r32(k) - a21(k) * a31(k) / a22(k)
@@ -1867,5 +2070,5 @@
     DO  k = nzb+1, nzt+1
        a33(k) = r33(k) - a31(k)**2 - a32(k)**2
-       IF ( a33(k) > 0.0_wp )  THEN
+       IF ( a33(k) > 10E-6_wp )  THEN
           a33(k) =  SQRT( a33(k) )
        ELSE
@@ -1897,7 +2100,13 @@
     IF ( debug_output_timestep )  CALL debug_message( 'stg_adjust', 'start' )
 !
-!-- Compute mean boundary layer height according to Richardson-Bulk
-!-- criterion using the inflow profiles. Further velocity scale as well as
-!-- mean friction velocity are calculated.
+!-- In case of dirichlet inflow boundary conditions only at one lateral 
+!-- boundary, i.e. in the case no offline or self nesting is applied but 
+!-- synthetic turbulence shall be parametrized nevertheless, the 
+!-- boundary-layer depth need to determined first.
+    IF ( .NOT. nesting_offline  .AND.  .NOT.  child_domain )                   &
+       CALL nesting_offl_calc_zi    
+!
+!-- Compute scaling parameters (domain-averaged), such as friction velocity
+!-- are calculated.
     CALL calc_scaling_variables
 !
@@ -1942,16 +2151,10 @@
  SUBROUTINE calc_length_and_time_scale
 
-    USE arrays_3d,                                                             &
-        ONLY:  dzw, ddzw, u_init, v_init
-
-    USE grid_variables,                                                        &
-        ONLY:  ddx, ddy, dx, dy
-
     IMPLICIT NONE
 
 
-    INTEGER(iwp) :: k            !< loop index in z-direction
-
-    REAL(wp) ::  length_scale    !< typical length scale
+    INTEGER(iwp) ::  k !< loop index in z-direction
+    
+    REAL(wp) ::  length_scale_dum     !< effectively used length scale
     
 !
@@ -1970,51 +2173,44 @@
 !-- using the length scales and the mean profiles of u- and v-component.
     DO  k = nzb+1, nzt+1
-    
-       length_scale = 8.0_wp * MIN( dx, dy, dzw(k) )
-          
-       IF ( zu(k) <= zi_ribulk )  THEN  
-!
-!--       Assume isotropic turbulence length scales
-          nux(k) = MAX( INT( length_scale * ddx     ), 1 )
-          nuy(k) = MAX( INT( length_scale * ddy     ), 1 )
-          nuz(k) = MAX( INT( length_scale * ddzw(k) ), 1 )
-          nvx(k) = MAX( INT( length_scale * ddx     ), 1 )
-          nvy(k) = MAX( INT( length_scale * ddy     ), 1 )
-          nvz(k) = MAX( INT( length_scale * ddzw(k) ), 1 )
-          nwx(k) = MAX( INT( length_scale * ddx     ), 1 )
-          nwy(k) = MAX( INT( length_scale * ddy     ), 1 )
-          nwz(k) = MAX( INT( length_scale * ddzw(k) ), 1 )
-!
-!--       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 ) )
-          tv(k)  = MIN( dt_stg_adjust, length_scale /                          &
-                                  ( ABS( v_init(k) ) + 0.1_wp ) )
-!
-!--       Time scale of w-component is a mixture from u- and v-component. 
-          tw(k)  = SQRT( tu(k)**2 + tv(k)**2 )
-!
-!--    Above the boundary layer length scales are zero, i.e. imposed turbulence
-!--    is not correlated in space and time, just white noise. This saves
-!--    computations power. 
-       ELSE
-          nux(k) = 0.0_wp
-          nuy(k) = 0.0_wp
-          nuz(k) = 0.0_wp
-          nvx(k) = 0.0_wp
-          nvy(k) = 0.0_wp
-          nvz(k) = 0.0_wp
-          nwx(k) = 0.0_wp
-          nwy(k) = 0.0_wp
-          nwz(k) = 0.0_wp
-          
-          tu(k)  = 0.0_wp
-          tv(k)  = 0.0_wp
-          tw(k)  = 0.0_wp
-       ENDIF
+!
+!--    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
 !
@@ -2045,10 +2241,4 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE calc_scaling_variables
-             
-    USE arrays_3d,                                                          &
-        ONLY:  drho_air         
-             
-    USE indices,                                                               &
-        ONLY:  nx, ny  
         
     USE surface_mod,                                                           &
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 4020)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 4022)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! Call synthetic turbulence generator at last RK3 substep right after boundary
+! conditions are updated in offline nesting in order to assure that 
+! perturbations are always imposed 
+! 
+! 4017 2019-06-06 12:16:46Z schwenkel
 ! Mass (volume) flux correction included to ensure global mass conservation for child domains.
 ! 
@@ -1240,8 +1245,7 @@
              CALL nesting_offl_bc
 !
-!--       Impose a turbulent inflow using synthetic generated turbulence, 
-!--       only once per time step.
-          IF ( use_syn_turb_gen  .AND.  time_stg_call >= dt_stg_call  .AND.                        &
-             intermediate_timestep_count == intermediate_timestep_count_max )  THEN
+!--       Impose a turbulent inflow using synthetic generated turbulence.
+          IF ( use_syn_turb_gen  .AND.                                                             &
+               intermediate_timestep_count == intermediate_timestep_count_max )  THEN
              CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
              CALL stg_main
