Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 2934)
+++ palm/trunk/SOURCE/Makefile	(revision 2938)
@@ -25,4 +25,9 @@
 # -----------------
 # $Id$
+# No initialization of child domains via dynamic input file, except for soil
+# moisture and temperature
+# Apply turbulence generator at non-cyclic lateral boundary in nesting case
+# 
+# 2936 2018-03-27 14:49:27Z suehring
 # Added dependencies for parent and child synchronization
 # 
@@ -1008,4 +1013,5 @@
 	modules.o \
 	netcdf_data_input_mod.o \
+	pmc_interface_mod.o \
 	radiation_model_mod.o \
 	surface_mod.o
@@ -1434,5 +1440,6 @@
 	cpulog_mod.o \
 	mod_kinds.o \
-	modules.o
+	modules.o \
+	pmc_interface_mod.o
 temperton_fft_mod.o: \
 	mod_kinds.o \
@@ -1510,4 +1517,5 @@
 	modules.o \
 	plant_canopy_model_mod.o \
+	pmc_interface_mod.o \
 	surface_mod.o \
 	user_actions.o
Index: palm/trunk/SOURCE/boundary_conds.f90
===================================================================
--- palm/trunk/SOURCE/boundary_conds.f90	(revision 2934)
+++ palm/trunk/SOURCE/boundary_conds.f90	(revision 2938)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! Set boundary condition for TKE and TKE dissipation rate in case of nesting
+! and if parent model operates in RANS mode but child model in LES mode.
+! mode
+! 
+! 2793 2018-02-07 10:54:33Z suehring
 ! Removed preprocessor directive __chem
 ! 
@@ -193,7 +198,7 @@
                inflow_s, intermediate_timestep_count,                          &
                microphysics_morrison, microphysics_seifert, nest_domain,       &
-               nest_bound_l, nest_bound_s, nudging, ocean, outflow_l,          &
-               outflow_n, outflow_r, outflow_s, passive_scalar, rans_tke_e,    &
-               tsc, use_cmax
+               nest_bound_l, nest_bound_n, nest_bound_r, nest_bound_s, nudging,&
+               ocean, outflow_l, outflow_n, outflow_r, outflow_s,              &
+               passive_scalar, rans_mode, rans_tke_e, tsc, use_cmax
 
     USE grid_variables,                                                        &
@@ -209,5 +214,5 @@
 
     USE pmc_interface,                                                         &
-        ONLY : nesting_mode
+        ONLY : nesting_mode, rans_mode_parent
 
     USE surface_mod,                                                           &
@@ -321,6 +326,6 @@
 
 !
-!-- Boundary conditions for TKE
-!-- Generally Neumann conditions with de/dz=0 are assumed
+!-- Boundary conditions for TKE.
+!-- Generally Neumann conditions with de/dz=0 are assumed.
     IF ( .NOT. constant_diffusion )  THEN
 
@@ -343,9 +348,28 @@
        IF ( .NOT. nest_domain )  THEN
           e_p(nzt+1,:,:) = e_p(nzt,:,:)
-       ENDIF
-    ENDIF
-
-!
-!-- Boundary conditions for TKE dissipation rate
+!
+!--    Nesting case: if parent operates in RANS mode and child in LES mode,
+!--    no TKE is transfered. This case, set Neumann conditions at lateral and 
+!--    top child boundaries. 
+!--    If not ( both either in RANS or in LES mode ), TKE boundary condition
+!--    is treated in the nesting. 
+       ELSE
+
+          IF ( rans_mode_parent  .AND.  .NOT. rans_mode )  THEN
+
+
+
+             e_p(nzt+1,:,:) = e_p(nzt,:,:)
+             IF ( nest_bound_l )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
+             IF ( nest_bound_r )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
+             IF ( nest_bound_s )  e_p(:,nys-1,:) = e_p(:,nys,:)
+             IF ( nest_bound_n )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
+
+          ENDIF
+       ENDIF
+    ENDIF
+
+!
+!-- Boundary conditions for TKE dissipation rate. 
     IF ( rans_tke_e .AND. .NOT. nest_domain )  THEN
        diss_p(nzt+1,:,:) = diss_p(nzt,:,:)
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 2934)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 2938)
@@ -25,4 +25,10 @@
 ! -----------------
 ! $Id$
+! - Revise start and end indices for imposing random disturbances in case of
+!   nesting or large-scale forcing.
+! - Remove check for inifor initialization in case of nested runs. 
+! - Adapt call to check for synthetic turbulence geneartor settings.
+! 
+! 2936 2018-03-27 14:49:27Z suehring
 ! Check spinup in case of nested runs, spinup time and timestep must be 
 ! identical to assure synchronuous simulations. 
@@ -1382,26 +1388,4 @@
     ENDIF
 !
-!-- In case of nested run assure that all domains are initialized the same
-!-- way, i.e. if at least at one domain is initialized with soil and 
-!-- atmospheric data provided by COSMO, all domains must be initialized the 
-!-- same way, to assure that soil and atmospheric quantities are 
-!-- consistent. 
-    IF ( nested_run )  THEN
-       check_nest = .TRUE.
-#if defined( __parallel )
-       CALL MPI_ALLREDUCE( TRIM( initializing_actions ) == 'inifor',           &
-                           check_nest, 1, MPI_LOGICAL,                         &
-                           MPI_LAND, MPI_COMM_WORLD, ierr )
-
-       IF ( TRIM( initializing_actions ) == 'inifor'  .AND.                    &
-            .NOT.  check_nest )  THEN
-          message_string = 'In case of nesting, if at least in one ' //        &
-                           'domain initializing_actions = inifor, '  //        &
-                           'all domains need to be initialized that way.'
-          CALL message( 'netcdf_data_input_mod', 'PA0430', 3, 2, 0, 6, 0 )
-       ENDIF
-#endif
-    ENDIF
-!
 !-- In case of spinup and nested run, spinup end time must be identical 
 !-- in order to have synchronously running simulations. 
@@ -1436,6 +1420,6 @@
 
 !
-!-- When synthetic turbulence generator is used, peform addtional checks
-    IF ( syn_turb_gen )  CALL stg_check_parameters
+!-- Check for synthetic turbulence generator settings
+    CALL stg_check_parameters
 !
 !-- When plant canopy model is used, peform addtional checks
@@ -3898,5 +3882,6 @@
           dist_nxr(1) = MIN( inflow_disturbance_end, nxr )
        ELSEIF ( bc_lr == 'nested'  .OR.  bc_lr == 'forcing' )  THEN
-          dist_nxl = MAX( inflow_disturbance_begin, nxl )
+          dist_nxl    = MAX( inflow_disturbance_begin, nxl )
+          dist_nxr    = MIN( nx - inflow_disturbance_begin, nxr )
        ENDIF
        IF ( bc_ns == 'dirichlet/radiation' )  THEN
@@ -3907,5 +3892,6 @@
           dist_nyn(1) = MIN( inflow_disturbance_end, nyn )
        ELSEIF ( bc_ns == 'nested'  .OR.  bc_ns == 'forcing' )  THEN
-          dist_nys = MAX( inflow_disturbance_begin, nys )
+          dist_nys    = MAX( inflow_disturbance_begin, nys )
+          dist_nyn    = MIN( ny - inflow_disturbance_begin, nyn )
        ENDIF
     ELSE
@@ -3918,4 +3904,7 @@
           dist_nxl    = inflow_disturbance_begin
           dist_nxr(1) = inflow_disturbance_end
+       ELSEIF ( bc_lr == 'nested'  .OR.  bc_lr == 'forcing' )  THEN
+          dist_nxr    = nx - inflow_disturbance_begin
+          dist_nxl    = inflow_disturbance_begin
        ENDIF
        IF ( bc_ns == 'dirichlet/radiation' )  THEN
@@ -3925,4 +3914,7 @@
           dist_nys    = inflow_disturbance_begin
           dist_nyn(1) = inflow_disturbance_end
+       ELSEIF ( bc_ns == 'nested'  .OR.  bc_ns == 'forcing' )  THEN
+          dist_nyn    = ny - inflow_disturbance_begin
+          dist_nys    = inflow_disturbance_begin
        ENDIF
     ENDIF
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 2934)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 2938)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! - Revise Inifor initialization for geostrophic wind components
+! - Initialize synthetic turbulence generator in case of Inifor initialization  
+!
+! 2936 2018-03-27 14:49:27Z suehring
 ! Synchronize parent and child models after initialization.
 ! Remove obsolete masking of topography grid points for Runge-Kutta weighted
@@ -1189,6 +1193,14 @@
 !--           are not identical to the initial wind profiles in this case. 
 !--           This need to be further revised.
-!           ug(:) = u_init(:)
-!           vg(:) = v_init(:)
+          IF ( init_3d%from_file_ug )  THEN
+             ug(:) = init_3d%ug_init(:)
+          ENDIF
+          IF ( init_3d%from_file_vg )  THEN
+             vg(:) = init_3d%vg_init(:)
+          ENDIF
+
+          ug(nzt+1) = ug(nzt)
+          vg(nzt+1) = vg(nzt)
+
 !
 !--       Set inital w to 0
@@ -1243,4 +1255,7 @@
 !--       fluxes, etc. 
           CALL init_surfaces
+!
+!--       Initialize turbulence generator
+          IF( use_syn_turb_gen )  CALL stg_init
 
           CALL location_message( 'finished', .TRUE. )
@@ -1335,7 +1350,5 @@
 !
 !--       Overwrite initial profiles in case of synthetic turbulence generator
-          IF( use_syn_turb_gen ) THEN
-             CALL stg_init
-          ENDIF
+          IF( use_syn_turb_gen )  CALL stg_init
 
 !
Index: palm/trunk/SOURCE/init_pegrid.f90
===================================================================
--- palm/trunk/SOURCE/init_pegrid.f90	(revision 2934)
+++ palm/trunk/SOURCE/init_pegrid.f90	(revision 2938)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! - No checks for domain decomposition in case of turbulence generator 
+!  (is done in stg module)
+! - Introduce ids to indicate lateral processors for turbulence generator
+! 
+! 2936 2018-03-27 14:49:27Z suehring
 ! Variable use_synthetic_turbulence_generator has been abbreviated
 ! 
@@ -251,5 +256,6 @@
 
     USE synthetic_turbulence_generator_mod,                                    &
-        ONLY:  use_syn_turb_gen
+        ONLY:  id_stg_left, id_stg_north, id_stg_right, id_stg_south,          &
+               use_syn_turb_gen
 
     USE transpose_indices,                                                     &
@@ -267,4 +273,8 @@
     INTEGER(iwp) ::  id_outflow_source_l      !< local value of id_outflow_source
     INTEGER(iwp) ::  id_recycling_l           !<
+    INTEGER(iwp) ::  id_stg_left_l            !< left lateral boundary local core id in case of turbulence generator  
+    INTEGER(iwp) ::  id_stg_north_l           !< north lateral boundary local core id in case of turbulence generator  
+    INTEGER(iwp) ::  id_stg_right_l           !< right lateral boundary local core id in case of turbulence generator  
+    INTEGER(iwp) ::  id_stg_south_l           !< south lateral boundary local core id in case of turbulence generator  
     INTEGER(iwp) ::  ind(5)                   !<
     INTEGER(iwp) ::  j                        !<
@@ -518,6 +528,5 @@
 !-- 1. transposition  z --> x
 !-- This transposition is not neccessary in case of a 1d-decomposition along x
-    IF ( psolver == 'poisfft'  .OR.  calculate_spectra  .OR.                   &
-         use_syn_turb_gen )  THEN
+    IF ( psolver == 'poisfft'  .OR.  calculate_spectra )  THEN
 
        IF ( pdims(2) /= 1 )  THEN
@@ -1263,5 +1272,48 @@
        ENDIF
     ENDIF
-
+!
+!-- In case of synthetic turbulence geneartor determine ids. 
+!-- Please note, if no forcing or nesting is applied, the generator is applied
+!-- only at the left lateral boundary.
+    IF ( use_syn_turb_gen )  THEN
+       IF ( force_bound_l  .OR.  nest_bound_l  .OR.  inflow_l )  THEN
+          id_stg_left_l = myidx
+       ELSE
+          id_stg_left_l = 0
+       ENDIF
+       IF ( force_bound_r  .OR.  nest_bound_r )  THEN
+          id_stg_right_l = myidx
+       ELSE
+          id_stg_right_l = 0
+       ENDIF
+       IF ( force_bound_s  .OR.  nest_bound_s )  THEN
+          id_stg_south_l = myidy
+       ELSE
+          id_stg_south_l = 0
+       ENDIF
+       IF ( force_bound_n  .OR.  nest_bound_n )  THEN
+          id_stg_north_l = myidy
+       ELSE
+          id_stg_north_l = 0
+       ENDIF
+
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( id_stg_left_l, id_stg_left,   1, MPI_INTEGER,       &
+                           MPI_SUM, comm1dx, ierr )
+
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( id_stg_right_l, id_stg_right, 1, MPI_INTEGER,       &
+                           MPI_SUM, comm1dx, ierr )
+
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( id_stg_south_l, id_stg_south, 1, MPI_INTEGER,       &
+                           MPI_SUM, comm1dy, ierr )
+
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( id_stg_north_l, id_stg_north, 1, MPI_INTEGER,       &
+                           MPI_SUM, comm1dy, ierr )
+
+    ENDIF 
+ 
 !
 !-- Broadcast the id of the inflow PE
Index: palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2934)
+++ palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2938)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! Initialization of soil moisture and temperature via Inifor-provided data also
+! in nested child domains, even if no dynamic input file is available for 
+! child domain. 1D soil profiles are received from parent model.  
+! 
+! 2936 2018-03-27 14:49:27Z suehring
 ! renamed lsm_par to land_surface_parameters. Bugfix in message calls
 ! 
@@ -2365,6 +2370,14 @@
        USE control_parameters,                                                 &
            ONLY:  message_string
+
+       USE indices,                                                            &
+           ONLY:  nx, ny
+
+       USE pmc_interface,                                                      &
+           ONLY:  nested_run
     
        IMPLICIT NONE
+
+       LOGICAL      ::  init_soil_dynamically_in_child !< flag controlling initialization of soil in child domains
 
        INTEGER(iwp) ::  i                       !< running index
@@ -2384,4 +2397,5 @@
 
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  bound, bound_root_fr  !< temporary arrays for storing index bounds
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pr_soil_init !< temporary array used for averaging soil profiles
 
 !
@@ -3894,4 +3908,19 @@
        IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
 !
+!--       In case of nested runs, check if soil is initialized via dynamic
+!--       input file in root domain and distribute this information
+!--       to all embedded child domains. This case, soil moisture and
+!--       temperature will be initialized via root domain.  
+          init_soil_dynamically_in_child = .FALSE. 
+          IF ( nested_run )  THEN
+#if defined( __parallel )
+             CALL MPI_ALLREDUCE( init_3d%from_file_tsoil  .OR.                 &
+                                 init_3d%from_file_msoil,                      &
+                                 init_soil_dynamically_in_child,               &
+                                 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr )
+#endif
+          ENDIF
+
+!
 !--       First, initialize soil temperature and moisture. 
 !--       According to the initialization for surface and soil parameters, 
@@ -3923,10 +3952,97 @@
           ENDDO
 !
-!--       Level 2, if soil moisture and/or temperature  are 
-!--       provided from file, interpolate / extrapolate the provided data
-!--       onto the respective soil layers. Please note, both, zs as well as 
-!--       init_3d%z_soil indicate a depth with positive values, so that no 
-!--       distinction between atmosphere is required concerning interpolation.
-!--       Start with soil moisture
+!--       Initialization of soil moisture and temperature from file. 
+!--       In case of nested runs, only root parent reads dynamic input file. 
+!--       This case, transfer respective soil information provide
+!--       by dynamic input file from root parent domain onto all other domains.
+          IF ( init_soil_dynamically_in_child )  THEN
+!
+!--          Child domains will be only initialized with horizontally 
+!--          averaged soil profiles in parent domain (for sake of simplicity). 
+!--          If required, average soil data on root parent domain before 
+!--          distribute onto child domains.
+             IF ( init_3d%from_file_msoil  .AND.  init_3d%lod_msoil == 2 )     &
+             THEN
+                ALLOCATE( pr_soil_init(0:init_3d%nzs-1) )
+
+                DO  k = 0, init_3d%nzs-1
+                   pr_soil_init(k) = SUM( init_3d%msoil(k,nys:nyn,nxl:nxr)  )
+                ENDDO
+!
+!--             Allocate 1D array for soil-moisture profile (will not be 
+!--             allocated in lod==2 case). 
+                ALLOCATE( init_3d%msoil_init(0:init_3d%nzs-1) )
+                init_3d%msoil_init = 0.0_wp
+#if defined( __parallel )
+                CALL MPI_ALLREDUCE( pr_soil_init(0), init_3d%msoil_init(0),    &
+                                    SIZE(pr_soil_init),                        &
+                                    MPI_REAL, MPI_SUM, comm2d, ierr )
+#endif
+                init_3d%msoil_init = init_3d%msoil_init /                      &
+                                        REAL( ( nx + 1 ) * ( ny + 1), KIND=wp )
+                DEALLOCATE( pr_soil_init )
+             ENDIF
+             IF ( init_3d%from_file_tsoil  .AND.  init_3d%lod_tsoil == 2 )  THEN
+                ALLOCATE( pr_soil_init(0:init_3d%nzs-1) )
+
+                DO  k = 0, init_3d%nzs-1
+                   pr_soil_init(k) = SUM( init_3d%tsoil(k,nys:nyn,nxl:nxr)  )
+                ENDDO
+!
+!--             Allocate 1D array for soil-temperature profile (will not be 
+!--             allocated in lod==2 case). 
+                ALLOCATE( init_3d%tsoil_init(0:init_3d%nzs-1) )
+                init_3d%tsoil_init = 0.0_wp
+#if defined( __parallel )
+                CALL MPI_ALLREDUCE( pr_soil_init(0), init_3d%tsoil_init(0),    &
+                                    SIZE(pr_soil_init),                        &
+                                    MPI_REAL, MPI_SUM, comm2d, ierr )
+#endif
+                init_3d%tsoil_init = init_3d%tsoil_init /                      &
+                                        REAL( ( nx + 1 ) * ( ny + 1), KIND=wp )
+                DEALLOCATE( pr_soil_init )
+
+             ENDIF
+
+#if defined( __parallel )
+!
+!--          Distribute soil grid information on file from root to all childs.
+!--          Only process with rank 0 sends the information. 
+             CALL MPI_BCAST( init_3d%nzs,    1,                                &
+                             MPI_INTEGER, 0, MPI_COMM_WORLD, ierr )
+
+             IF ( .NOT.  ALLOCATED( init_3d%z_soil ) )                         &
+                ALLOCATE( init_3d%z_soil(1:init_3d%nzs) )
+
+             CALL MPI_BCAST( init_3d%z_soil, SIZE(init_3d%z_soil),             &
+                             MPI_REAL, 0, MPI_COMM_WORLD, ierr )
+#endif
+!
+!--          ALLOCATE arrays on child domains and set control attributes.
+             IF ( .NOT. init_3d%from_file_msoil )  THEN
+                ALLOCATE( init_3d%msoil_init(0:init_3d%nzs-1) )
+                init_3d%lod_msoil = 1
+                init_3d%from_file_msoil = .TRUE.
+             ENDIF
+             IF ( .NOT. init_3d%from_file_tsoil )  THEN
+                ALLOCATE( init_3d%tsoil_init(0:init_3d%nzs-1) )
+                init_3d%lod_tsoil = 1
+                init_3d%from_file_tsoil = .TRUE.
+             ENDIF
+
+
+#if defined( __parallel )
+!
+!--          Distribute soil profiles from root to all childs
+             CALL MPI_BCAST( init_3d%msoil_init, SIZE(init_3d%msoil_init),     &
+                             MPI_REAL, 0, MPI_COMM_WORLD, ierr )
+             CALL MPI_BCAST( init_3d%tsoil_init, SIZE(init_3d%tsoil_init),     &
+                             MPI_REAL, 0, MPI_COMM_WORLD, ierr )
+#endif
+
+          ENDIF
+!
+!--       Proceed with Level 2 initialization. Information from dynamic input
+!--       is now available on all processes. 
           IF ( init_3d%from_file_msoil )  THEN
 
@@ -3983,5 +4099,4 @@
                 ENDDO
              ENDIF
-
           ENDIF
 !
Index: palm/trunk/SOURCE/large_scale_forcing_nudging_mod.f90
===================================================================
--- palm/trunk/SOURCE/large_scale_forcing_nudging_mod.f90	(revision 2934)
+++ palm/trunk/SOURCE/large_scale_forcing_nudging_mod.f90	(revision 2938)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Further improvements for nesting in larger-scale model
+! 
+! 2863 2018-03-08 11:36:25Z suehring
 ! Corrected "Former revisions" section
 ! 
@@ -210,5 +213,5 @@
        ENDIF
 !
-!--       Top boundary
+!--    Top boundary
        k = nzt
        DO  i = nxl, nxr
@@ -235,4 +238,7 @@
           ENDDO
        ENDDO
+
+write(9,*) "w correction", w_correct
+flush(9)
 
     END SUBROUTINE forcing_bc_mass_conservation
@@ -581,50 +587,91 @@
 !--    identical to the Inifor grid. At the top boundary an extrapolation is 
 !--    not possible. 
-       IF ( ANY( zu(1:nzt+1) /= force%zu_atmos(1:force%nzu) ) )  THEN
-          DO  i = nxlu, nxr
-             DO  j = nys, nyn
-                u(nzt+1,j,i) = force%u_top(0,j,i) + ddt_lsf * t_ref *          &
-                           ( force%u_top(1,j,i) - force%u_top(0,j,i) ) *       &
-                              MERGE( 1.0_wp, 0.0_wp,                           &
-                                     BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
-             ENDDO
-          ENDDO
-
-          DO  i = nxl, nxr
-             DO  j = nysv, nyn
-                v(nzt+1,j,i) = force%v_top(0,j,i) + ddt_lsf * t_ref *          &
-                           ( force%v_top(1,j,i) - force%v_top(0,j,i) ) *       &
-                              MERGE( 1.0_wp, 0.0_wp,                           &
-                                     BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
-             ENDDO
-          ENDDO
-
+       DO  i = nxlu, nxr
+          DO  j = nys, nyn
+             u(nzt+1,j,i) = force%u_top(0,j,i) + ddt_lsf * t_ref *             &
+                        ( force%u_top(1,j,i) - force%u_top(0,j,i) ) *          &
+                           MERGE( 1.0_wp, 0.0_wp,                             &
+                                  BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
+          ENDDO
+       ENDDO
+
+       DO  i = nxl, nxr
+          DO  j = nysv, nyn
+             v(nzt+1,j,i) = force%v_top(0,j,i) + ddt_lsf * t_ref *             &
+                        ( force%v_top(1,j,i) - force%v_top(0,j,i) ) *          &
+                           MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
+          ENDDO
+       ENDDO
+
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             w(nzt:nzt+1,j,i) = force%w_top(0,j,i) + ddt_lsf * t_ref *         &
+                        ( force%w_top(1,j,i) - force%w_top(0,j,i) ) *          &
+                           MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(nzt:nzt+1,j,i), 3 ) )
+          ENDDO
+       ENDDO
+
+
+       IF ( .NOT. neutral )  THEN
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                w(nzt:nzt+1,j,i) = force%w_top(0,j,i) + ddt_lsf * t_ref *      &
-                           ( force%w_top(1,j,i) - force%w_top(0,j,i) ) *       &
-                              MERGE( 1.0_wp, 0.0_wp,                           &
-                                     BTEST( wall_flags_0(nzt:nzt+1,j,i), 3 ) )
-             ENDDO
-          ENDDO
-
-
-          IF ( .NOT. neutral )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   pt(nzt+1,j,i) = force%pt_top(0,j,i) + ddt_lsf * t_ref *     &
-                           ( force%pt_top(1,j,i) - force%pt_top(0,j,i) )
-                ENDDO
-             ENDDO
-          ENDIF
-
-          IF ( humidity )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   q(nzt+1,j,i) = force%q_top(0,j,i) + ddt_lsf * t_ref *       &
-                           ( force%q_top(1,j,i) - force%q_top(0,j,i) )
-                ENDDO
-             ENDDO
-          ENDIF
+                pt(nzt+1,j,i) = force%pt_top(0,j,i) + ddt_lsf * t_ref *        &
+                        ( force%pt_top(1,j,i) - force%pt_top(0,j,i) )
+             ENDDO
+          ENDDO
+       ENDIF
+
+       IF ( humidity )  THEN
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                q(nzt+1,j,i) = force%q_top(0,j,i) + ddt_lsf * t_ref *          &
+                        ( force%q_top(1,j,i) - force%q_top(0,j,i) )
+             ENDDO
+          ENDDO
+       ENDIF
+!
+!--    At the edges( left-south, left-north, right-south and right-north) set 
+!--    data on ghost points. 
+       IF ( force_bound_l  .AND.  force_bound_s )  THEN
+          DO  i = 1, nbgp
+             u(:,nys-i,nxlg:nxl)   = u(:,nys,nxlg:nxl)
+             w(:,nys-i,nxlg:nxl-1) = w(:,nys,nxlg:nxl-1)
+             IF ( .NOT. neutral )  pt(:,nys-i,nxlg:nxl-1) = pt(:,nys,nxlg:nxl-1)
+             IF ( humidity )       q(:,nys-i,nxlg:nxl-1)  = q(:,nys,nxlg:nxl-1)
+          ENDDO
+          DO  i = 1, nbgp+1
+             v(:,nysv-i,nxlg:nxl-1) = v(:,nysv,nxlg:nxl-1)
+          ENDDO
+       ENDIF
+       IF ( force_bound_l  .AND.  force_bound_n )  THEN
+          DO  i = 1, nbgp
+             u(:,nyn+i,nxlg:nxl)   = u(:,nyn,nxlg:nxl)
+             v(:,nyn+i,nxlg:nxl-1) = v(:,nyn,nxlg:nxl-1)
+             w(:,nyn+i,nxlg:nxl-1) = w(:,nyn,nxlg:nxl-1)
+             IF ( .NOT. neutral )  pt(:,nyn+i,nxlg:nxl-1) = pt(:,nyn,nxlg:nxl-1)
+             IF ( humidity )       q(:,nyn+i,nxlg:nxl-1)  = q(:,nyn,nxlg:nxl-1)
+          ENDDO
+       ENDIF
+       IF ( force_bound_r  .AND.  force_bound_s )  THEN
+          DO  i = 1, nbgp
+             u(:,nys-i,nxr+1:nxrg) = u(:,nys,nxr+1:nxrg)
+             w(:,nys-i,nxr+1:nxrg) = w(:,nys,nxr+1:nxrg)
+             IF ( .NOT. neutral )  pt(:,nys-i,nxr+1:nxrg) = pt(:,nys,nxr+1:nxrg)
+             IF ( humidity )       q(:,nys-i,nxr+1:nxrg)  = q(:,nys,nxr+1:nxrg)
+          ENDDO
+          DO  i = 1, nbgp+1
+             v(:,nysv-i,nxr+1:nxrg) = v(:,nysv,nxr+1:nxrg)
+          ENDDO
+       ENDIF
+       IF ( force_bound_r  .AND.  force_bound_n )  THEN
+          DO  i = 1, nbgp
+             u(:,nyn+i,nxr+1:nxrg) = u(:,nyn,nxr+1:nxrg)
+             v(:,nyn+i,nxr+1:nxrg) = v(:,nyn,nxr+1:nxrg)
+             w(:,nyn+i,nxr+1:nxrg) = w(:,nyn,nxr+1:nxrg)
+             IF ( .NOT. neutral )  pt(:,nyn+i,nxr+1:nxrg) = pt(:,nyn,nxr+1:nxrg)
+             IF ( humidity )       q(:,nyn+i,nxr+1:nxrg)  = q(:,nyn,nxr+1:nxrg)
+          ENDDO
        ENDIF
 !
@@ -941,4 +988,7 @@
     SUBROUTINE lsf_init
 
+       USE control_parameters,                                                 &
+           ONLY:  bc_lr_cyc, bc_ns_cyc
+
        USE netcdf_data_input_mod,                                              &
            ONLY:  netcdf_data_input_lsf  
@@ -975,4 +1025,13 @@
 
        IF ( forcing )  THEN
+!
+!--       Allocate arrays for geostrophic wind components. Arrays will 
+!--       incorporate 2 time levels in order to interpolate in between. Please
+!--       note, forcing using geostrophic wind components is only required in 
+!--       case of cyclic boundary conditions.
+          IF ( bc_lr_cyc  .AND.  bc_ns_cyc )  THEN
+             ALLOCATE( force%ug(0:1,nzb:nzt+1) )
+             ALLOCATE( force%vg(0:1,nzb:nzt+1) )
+          ENDIF
 !
 !--       Allocate arrays for reading boundary values. Arrays will incorporate 2 
Index: palm/trunk/SOURCE/netcdf_data_input_mod.f90
===================================================================
--- palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 2934)
+++ palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 2938)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Initial read of geostrophic wind components from dynamic driver. 
+! 
+! 2773 2018-01-30 14:12:54Z suehring
 ! Revise checks for surface_fraction.
 ! 
@@ -126,4 +129,7 @@
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  zw_atmos         !< vertical levels at w grid in dynamic input file
 
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ug         !< domain-averaged geostrophic component
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vg         !< domain-averaged geostrophic component
+
        REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u_left   !< u-component at left boundary
        REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  v_left   !< v-component at left boundary
@@ -180,6 +186,9 @@
        LOGICAL ::  from_file_tsoil  = .FALSE. !< flag indicating whether soil temperature is already initialized from file 
        LOGICAL ::  from_file_u      = .FALSE. !< flag indicating whether u is already initialized from file 
+       LOGICAL ::  from_file_ug     = .FALSE. !< flag indicating whether ug is already initialized from file 
        LOGICAL ::  from_file_v      = .FALSE. !< flag indicating whether v is already initialized from file 
-       LOGICAL ::  from_file_w      = .FALSE. !< flag indicating whether w is already initialized from file 
+       LOGICAL ::  from_file_vg     = .FALSE. !< flag indicating whether ug is already initialized from file 
+       LOGICAL ::  from_file_w      = .FALSE. !< flag indicating whether w is already initialized from file
+
 
        REAL(wp) ::  fill_msoil       !< fill value for soil moisture
@@ -422,4 +431,8 @@
        MODULE PROCEDURE get_variable_4d_real
     END INTERFACE get_variable
+
+    INTERFACE get_variable_pr
+       MODULE PROCEDURE get_variable_pr
+    END INTERFACE get_variable_pr
 
     INTERFACE get_attribute
@@ -2025,12 +2038,22 @@
              ENDIF
 !
-!--          Read geostrophic wind components
-!              IF ( check_existence( var_names, 'ls_forcing_ug' ) )  THEN
-! 
-!              ENDIF
-!              IF ( check_existence( var_names, 'ls_forcing_vg' ) )  THEN
-! 
-!              ENDIF
-
+!--          Read initial geostrophic wind components at 
+!--          t = 0 (index 1 in file). 
+             IF ( check_existence( var_names, 'ls_forcing_ug' ) )  THEN
+                ALLOCATE( init_3d%ug_init(nzb:nzt+1) )
+                CALL get_variable_pr( id_dynamic, 'ls_forcing_ug', 1,          &
+                                      init_3d%ug_init )
+                init_3d%from_file_ug = .TRUE.
+             ELSE
+                init_3d%from_file_ug = .FALSE.
+             ENDIF
+             IF ( check_existence( var_names, 'ls_forcing_vg' ) )  THEN
+                ALLOCATE( init_3d%vg_init(nzb:nzt+1) )
+                CALL get_variable_pr( id_dynamic, 'ls_forcing_vg', 1,          &
+                                      init_3d%vg_init )
+                init_3d%from_file_vg = .TRUE.
+             ELSE
+                init_3d%from_file_vg = .FALSE.
+             ENDIF
 !
 !--          Read inital 3D data of u, v, w, pt and q, 
@@ -2222,4 +2245,5 @@
 !--          Read soil moisture
              IF ( land_surface )  THEN
+
                 IF ( check_existence( var_names, 'init_soil_m' ) )  THEN
 !
@@ -2240,5 +2264,5 @@
 !
 !--                level-of-detail 2 - read 3D initialization data
-                   ELSEIF ( init_3d%lod_msoil == 2 )  THEN  ! need to be corrected
+                   ELSEIF ( init_3d%lod_msoil == 2 )  THEN
                       ALLOCATE ( init_3d%msoil(0:init_3d%nzs-1,nys:nyn,nxl:nxr) )
                       DO  i = nxl, nxr
@@ -2272,5 +2296,5 @@
 !
 !--                level-of-detail 2 - read 3D initialization data
-                   ELSEIF ( init_3d%lod_tsoil == 2 )  THEN  ! need to be corrected
+                   ELSEIF ( init_3d%lod_tsoil == 2 )  THEN
                       ALLOCATE ( init_3d%tsoil(0:init_3d%nzs-1,nys:nyn,nxl:nxr) )
                       DO  i = nxl, nxr
@@ -2346,6 +2370,8 @@
 
        USE control_parameters,                                                 &
-           ONLY:  force_bound_l, force_bound_n, force_bound_r, force_bound_s,  &
+           ONLY:  bc_lr_cyc, bc_ns_cyc, force_bound_l, force_bound_n,          &
+                  force_bound_r, force_bound_s,                                &
                   forcing, humidity, message_string, neutral, simulated_time
+
 
        USE indices,                                                            &
@@ -2430,5 +2456,16 @@
              force%tind = MINLOC( ABS( force%time - simulated_time ), DIM = 1 )&
                           - 1
-             force%tind_p = force%tind + 1  
+             force%tind_p = force%tind + 1
+!
+!--          Read geostrophic wind components. In case of forcing, this is only
+!--          required if cyclic boundary conditions are applied.
+             IF ( bc_lr_cyc  .AND.  bc_ns_cyc )  THEN
+                DO  t = force%tind, force%tind_p
+                   CALL get_variable_pr( id_dynamic, 'ls_forcing_ug', t+1,     &
+                                         force%ug(t-force%tind,:) )
+                   CALL get_variable_pr( id_dynamic, 'ls_forcing_vg', t+1,     &
+                                         force%ug(t-force%tind,:) )
+                ENDDO
+             ENDIF  
 !
 !--          Read data at lateral and top boundaries. Please note, at left and 
@@ -3832,4 +3869,45 @@
     END SUBROUTINE get_variable_1d_real
 
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Reads a time-dependent 1D float variable from file. 
+!------------------------------------------------------------------------------!
+    SUBROUTINE get_variable_pr( id, variable_name, t, var )
+#if defined( __netcdf )
+
+       USE pegrid
+
+       IMPLICIT NONE
+
+       CHARACTER(LEN=*)                      ::  variable_name    !< variable name 
+
+       INTEGER(iwp), INTENT(IN)              ::  id               !< file id
+       INTEGER(iwp), DIMENSION(1:2)          ::  id_dim           !< dimension ids
+       INTEGER(iwp)                          ::  id_var           !< dimension id
+       INTEGER(iwp)                          ::  n_file           !< number of data-points in file along z dimension
+       INTEGER(iwp), INTENT(IN)              ::  t                !< timestep number
+
+       REAL(wp), DIMENSION(:), INTENT(INOUT) ::  var  !< variable to be read
+
+!
+!--    First, inquire variable ID
+       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
+!
+!--    Inquire dimension size of vertical dimension
+       nc_stat = NF90_INQUIRE_VARIABLE( id, id_var, DIMIDS = id_dim )
+       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim(1), LEN = n_file )
+!
+!--    Read variable. 
+       nc_stat = NF90_GET_VAR( id, id_var, var,                                &
+                               start = (/ 1,      t     /),                    &
+                               count = (/ n_file, 1     /) )
+       CALL handle_error( 'get_variable_pr', 527 ) 
+
+#endif
+    END SUBROUTINE get_variable_pr
+
+
 !------------------------------------------------------------------------------!
 ! Description:
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 2934)
+++ palm/trunk/SOURCE/parin.f90	(revision 2938)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Change initialization in case child domain should be initialized with Inifor.
+! 
+! 2936 2018-03-27 14:49:27Z suehring
 ! inipar renamed to initialization_parameters.
 ! d3par renamed to runtime_parameters.
@@ -907,9 +910,7 @@
 !--       are set properly. An exception is made in case of restart runs and
 !--       if user decides to do everything by its own.
-!--       MS: is this really necessary?
           IF ( nest_domain  .AND.  .NOT. (                                     &
                TRIM( initializing_actions ) == 'read_restart_data'  .OR.       &
-               TRIM( initializing_actions ) == 'by_user'            .OR.       &
-               TRIM( initializing_actions ) == 'inifor' ) )  THEN
+               TRIM( initializing_actions ) == 'by_user' ) )  THEN
              initializing_actions = 'set_constant_profiles'
           ENDIF
Index: palm/trunk/SOURCE/pmc_interface_mod.f90
===================================================================
--- palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 2934)
+++ palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 2938)
@@ -25,4 +25,16 @@
 ! -----------------
 ! $Id$
+! - Nesting for RANS mode implemented
+!    - Interpolation of TKE onto child domain only if both parent and child are 
+!      either in LES mode or in RANS mode
+!    - Interpolation of dissipation if both parent and child are in RANS mode
+!      and TKE-epsilon closure is applied
+!    - Enable anterpolation of TKE and dissipation rate in case parent and 
+!      child operate in RANS mode
+!
+! - Some unused variables removed from ONLY list
+! - Some formatting adjustments for particle nesting
+! 
+! 2936 2018-03-27 14:49:27Z suehring
 ! Control logics improved to allow nesting also in cases with
 ! constant_flux_layer = .F. or constant_diffusion = .T.
@@ -254,11 +266,11 @@
 #if defined( __nopointer )
     USE arrays_3d,                                                             &
-        ONLY:  dzu, dzw, e, e_p, nc, nr, pt, pt_p, q, q_p, qc, qr, s, u, u_p,  &
+        ONLY:  diss, dzu, dzw, e, e_p, nc, nr, pt, q, qc, qr, s, u, u_p,       &
                v, v_p, w, w_p, zu, zw
 #else
    USE arrays_3d,                                                              &
-        ONLY:  dzu, dzw, e, e_p, e_1, e_2, nc, nc_2, nc_p, nr, nr_2, nr_p, pt, &
-               pt_p, pt_1, pt_2, q, q_p, q_1, q_2, qc, qc_2, qr, qr_2, s, s_2, &
-               u, u_p, u_1, u_2, v, v_p, v_1, v_2, w, w_p, w_1, w_2, zu, zw
+        ONLY:  diss, diss_2, dzu, dzw, e, e_p, e_2, nc, nc_2, nc_p, nr, nr_2,  &
+               pt, pt_2, q, q_2, qc, qc_2, qr, qr_2, s, s_2,                       &
+               u, u_p, u_2, v, v_p, v_2, w, w_p, w_2, zu, zw
 #endif
 
@@ -269,6 +281,6 @@
                microphysics_morrison, microphysics_seifert,                    &
                nest_bound_l, nest_bound_r, nest_bound_s, nest_bound_n,         &
-               nest_domain, neutral, passive_scalar, roughness_length,         & 
-               simulated_time, topography, volume_flow
+               nest_domain, neutral, passive_scalar, rans_mode, rans_tke_e,    &
+               roughness_length, simulated_time, topography, volume_flow
 
     USE chem_modules,                                                          &
@@ -359,4 +371,5 @@
 
     LOGICAL, SAVE ::  nested_run = .FALSE.  !< general switch
+    LOGICAL       ::  rans_mode_parent = .FALSE. !< mode of parent model (.F. - LES mode, .T. - RANS mode)
 
     REAL(wp), SAVE ::  anterp_relax_length_l = -1.0_wp   !<
@@ -387,4 +400,5 @@
     REAL(wp), SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_t    !<
 
+    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  dissc !< coarse grid array on child domain - dissipation rate
     REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ec   !<
     REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ptc  !<
@@ -401,6 +415,5 @@
     INTEGER(idp), SAVE, DIMENSION(:,:), ALLOCATABLE, TARGET, PUBLIC ::  part_adrc   !<
 
-
-    REAL(wp), SAVE, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  chem_spec_c   !< child coarse data array for chemical species 
+    REAL(wp), SAVE, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  chem_spec_c !< coarse grid array on child domain - chemical species
 
 !
@@ -616,5 +629,5 @@
            anterp_relax_length_t, child_to_parent, comm_world_nesting,         &
            cpl_id, nested_run, nesting_datatransfer_mode, nesting_mode,        &
-           parent_to_child
+           parent_to_child, rans_mode_parent
 
     PUBLIC pmci_boundary_conds
@@ -804,4 +817,5 @@
 !   Initialize the pmc parent
     CALL pmc_parentinit
+
 !
 !-- Corners of all children of the present parent
@@ -821,5 +835,6 @@
 
        child_id = pmc_parent_for_child(m)
-       IF ( myid == 0 )  THEN        
+
+       IF ( myid == 0 )  THEN
 
           CALL pmc_recv_from_child( child_id, val,  size(val),  0, 123, ierr )
@@ -932,4 +947,9 @@
           DEALLOCATE( cl_coord_x )
           DEALLOCATE( cl_coord_y )
+
+!
+!--       Send information about operating mode (LES or RANS) to child. This will be 
+!--       used to control TKE nesting and setting boundary conditions properly.
+          CALL pmc_send_to_child( child_id, rans_mode, 1, 0, 19, ierr )  
 !
 !--       Send coarse grid information to child
@@ -990,4 +1010,5 @@
           ENDIF
        ENDDO
+
        CALL pmc_s_setind_and_allocmem( child_id )
     ENDDO
@@ -1195,5 +1216,15 @@
        CALL pmc_set_dataarray_name( 'coarse', 'v'  ,'fine', 'v',  ierr )
        CALL pmc_set_dataarray_name( 'coarse', 'w'  ,'fine', 'w',  ierr )
+!
+!--    Set data array name for TKE. Please note, nesting of TKE is actually
+!--    only done if both parent and child are in LES or in RANS mode. Due to 
+!--    design of model coupler, however, data array names must be already 
+!--    available at this point.
        CALL pmc_set_dataarray_name( 'coarse', 'e'  ,'fine', 'e',  ierr )
+!
+!--    Nesting of dissipation rate only if both parent and child are in RANS
+!--    mode and TKE-epsilo closure is applied. Please so also comment for TKE
+!--    above.
+       CALL pmc_set_dataarray_name( 'coarse', 'diss'  ,'fine', 'diss',  ierr )
 
        IF ( .NOT. neutral )  THEN
@@ -1260,4 +1291,7 @@
           CALL pmc_send_to_parent( coord_x, nx + 1 + 2 * nbgp, 0, 11, ierr )
           CALL pmc_send_to_parent( coord_y, ny + 1 + 2 * nbgp, 0, 12, ierr )
+
+          CALL pmc_recv_from_parent( rans_mode_parent, 1, 0, 19, ierr )
+!
 !
 !--       Receive Coarse grid information.
@@ -1319,5 +1353,6 @@
        CALL MPI_BCAST( cg%zu, cg%nz+2,  MPI_REAL, 0, comm2d, ierr )
        CALL MPI_BCAST( cg%zw, cg%nz+2,  MPI_REAL, 0, comm2d, ierr )
-       
+       CALL MPI_BCAST( rans_mode_parent, 1, MPI_LOGICAL, 0, comm2d, ierr )
+  
 !
 !--    Find the index bounds for the nest domain in the coarse-grid index space
@@ -1357,8 +1392,11 @@
        ENDIF
 !
-!--    Define the SGS-TKE scaling factor based on the grid-spacing ratio
-       IF ( .NOT. constant_diffusion )  THEN
-          CALL pmci_init_tkefactor
-       ENDIF
+!--    Define the SGS-TKE scaling factor based on the grid-spacing ratio. Only
+!--    if both parent and child are in LES mode or in RANS mode. 
+!--    Please note, in case parent and child are in RANS mode, TKE weighting 
+!--    factor is simply one. 
+       IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.         &
+            (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.          &
+               .NOT. constant_diffusion ) )  CALL pmci_init_tkefactor
 !
 !--    Two-way coupling for general and vertical nesting.
@@ -2995,4 +3033,6 @@
 !--    energy spectrum. Near the surface, the reduction of TKE is made
 !--    smaller than further away from the surface.
+!--    Please note, in case parent and child model operate in RANS mode, 
+!--    TKE is not grid depenedent and weighting factor is one. 
 
        IMPLICIT NONE
@@ -3012,33 +3052,96 @@
        REAL(wp), PARAMETER ::  p23 = 2.0_wp/3.0_wp   !<        
 
-       IF ( nest_bound_l )  THEN
-          ALLOCATE( tkefactor_l(nzb:nzt+1,nysg:nyng) )
-          tkefactor_l = 0.0_wp
-          i = nxl - 1
-          DO  j = nysg, nyng
-             k_wall = get_topography_top_index_ji( j, i, 's' )
-
-             DO  k = k_wall + 1, nzt
-
-                kc     = kco(k) + 1
-                glsf   = ( dx * dy * dzu(k) )**p13
-                glsc   = ( cg%dx * cg%dy *cg%dzu(kc) )**p13
-                height = zu(k) - zu(k_wall)
-                fw     = EXP( -cfw * height / glsf )
-                tkefactor_l(k,j) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *     &
-                                              ( glsf / glsc )**p23 )
+!
+       IF ( .NOT. rans_mode  .AND.  .NOT. rans_mode_parent )  THEN 
+          IF ( nest_bound_l )  THEN
+             ALLOCATE( tkefactor_l(nzb:nzt+1,nysg:nyng) )
+             tkefactor_l = 0.0_wp
+             i = nxl - 1
+             DO  j = nysg, nyng
+                k_wall = get_topography_top_index_ji( j, i, 's' )
+
+                DO  k = k_wall + 1, nzt
+                   kc     = kco(k) + 1
+                   glsf   = ( dx * dy * dzu(k) )**p13
+                   glsc   = ( cg%dx * cg%dy *cg%dzu(kc) )**p13
+                   height = zu(k) - zu(k_wall)
+                   fw     = EXP( -cfw * height / glsf )
+                   tkefactor_l(k,j) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *      &
+                                                 ( glsf / glsc )**p23 )
+                ENDDO
+                tkefactor_l(k_wall,j) = c_tkef * fw0
              ENDDO
-             tkefactor_l(k_wall,j) = c_tkef * fw0
-          ENDDO
-       ENDIF
-
-       IF ( nest_bound_r )  THEN
-          ALLOCATE( tkefactor_r(nzb:nzt+1,nysg:nyng) )
-          tkefactor_r = 0.0_wp
-          i = nxr + 1
-          DO  j = nysg, nyng
-             k_wall = get_topography_top_index_ji( j, i, 's' )
-
-             DO  k = k_wall + 1, nzt
+          ENDIF
+
+          IF ( nest_bound_r )  THEN
+             ALLOCATE( tkefactor_r(nzb:nzt+1,nysg:nyng) )
+             tkefactor_r = 0.0_wp
+             i = nxr + 1
+             DO  j = nysg, nyng
+                k_wall = get_topography_top_index_ji( j, i, 's' )
+
+                DO  k = k_wall + 1, nzt
+                   kc     = kco(k) + 1
+                   glsf   = ( dx * dy * dzu(k) )**p13
+                   glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
+                   height = zu(k) - zu(k_wall)
+                   fw     = EXP( -cfw * height / glsf )
+                   tkefactor_r(k,j) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *      &
+                                                 ( glsf / glsc )**p23 )
+                ENDDO
+                tkefactor_r(k_wall,j) = c_tkef * fw0
+             ENDDO
+          ENDIF
+
+          IF ( nest_bound_s )  THEN
+             ALLOCATE( tkefactor_s(nzb:nzt+1,nxlg:nxrg) )
+             tkefactor_s = 0.0_wp
+             j = nys - 1
+             DO  i = nxlg, nxrg
+                k_wall = get_topography_top_index_ji( j, i, 's' )
+                
+                DO  k = k_wall + 1, nzt
+    
+                   kc     = kco(k) + 1
+                   glsf   = ( dx * dy * dzu(k) )**p13
+                   glsc   = ( cg%dx * cg%dy * cg%dzu(kc) ) ** p13
+                   height = zu(k) - zu(k_wall)
+                   fw     = EXP( -cfw*height / glsf )
+                   tkefactor_s(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *      &
+                        ( glsf / glsc )**p23 )
+                ENDDO
+                tkefactor_s(k_wall,i) = c_tkef * fw0
+             ENDDO
+          ENDIF
+
+          IF ( nest_bound_n )  THEN
+             ALLOCATE( tkefactor_n(nzb:nzt+1,nxlg:nxrg) )
+             tkefactor_n = 0.0_wp
+             j = nyn + 1
+             DO  i = nxlg, nxrg
+                k_wall = get_topography_top_index_ji( j, i, 's' )
+
+                DO  k = k_wall + 1, nzt
+
+                   kc     = kco(k) + 1
+                   glsf   = ( dx * dy * dzu(k) )**p13
+                   glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
+                   height = zu(k) - zu(k_wall)
+                   fw     = EXP( -cfw * height / glsf )
+                   tkefactor_n(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *     &
+                                                 ( glsf / glsc )**p23 )
+                ENDDO
+                tkefactor_n(k_wall,i) = c_tkef * fw0
+             ENDDO
+          ENDIF
+
+          ALLOCATE( tkefactor_t(nysg:nyng,nxlg:nxrg) )
+          k = nzt
+
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
+!
+!--             Determine vertical index for local topography top
+                k_wall = get_topography_top_index_ji( j, i, 's' )
 
                 kc     = kco(k) + 1
@@ -3047,71 +3150,32 @@
                 height = zu(k) - zu(k_wall)
                 fw     = EXP( -cfw * height / glsf )
-                tkefactor_r(k,j) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *     &
+                tkefactor_t(j,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *        &
                                               ( glsf / glsc )**p23 )
              ENDDO
-             tkefactor_r(k_wall,j) = c_tkef * fw0
           ENDDO
-       ENDIF
-
-       IF ( nest_bound_s )  THEN
-          ALLOCATE( tkefactor_s(nzb:nzt+1,nxlg:nxrg) )
-          tkefactor_s = 0.0_wp
-          j = nys - 1
-          DO  i = nxlg, nxrg
-             k_wall = get_topography_top_index_ji( j, i, 's' )
-             
-             DO  k = k_wall + 1, nzt
- 
-                kc     = kco(k) + 1
-                glsf   = ( dx * dy * dzu(k) )**p13
-                glsc   = ( cg%dx * cg%dy * cg%dzu(kc) ) ** p13
-                height = zu(k) - zu(k_wall)
-                fw     = EXP( -cfw*height / glsf )
-                tkefactor_s(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *     &
-                     ( glsf / glsc )**p23 )
-             ENDDO
-             tkefactor_s(k_wall,i) = c_tkef * fw0
-          ENDDO
-       ENDIF
-
-       IF ( nest_bound_n )  THEN
-          ALLOCATE( tkefactor_n(nzb:nzt+1,nxlg:nxrg) )
-          tkefactor_n = 0.0_wp
-          j = nyn + 1
-          DO  i = nxlg, nxrg
-             k_wall = get_topography_top_index_ji( j, i, 's' )
-
-             DO  k = k_wall + 1, nzt
-
-                kc     = kco(k) + 1
-                glsf   = ( dx * dy * dzu(k) )**p13
-                glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
-                height = zu(k) - zu(k_wall)
-                fw     = EXP( -cfw * height / glsf )
-                tkefactor_n(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *     &
-                                              ( glsf / glsc )**p23 )
-             ENDDO
-             tkefactor_n(k_wall,i) = c_tkef * fw0
-          ENDDO
-       ENDIF
-
-       ALLOCATE( tkefactor_t(nysg:nyng,nxlg:nxrg) )
-       k = nzt
-
-       DO  i = nxlg, nxrg
-          DO  j = nysg, nyng
-!
-!--          Determine vertical index for local topography top
-             k_wall = get_topography_top_index_ji( j, i, 's' )
-
-             kc     = kco(k) + 1
-             glsf   = ( dx * dy * dzu(k) )**p13
-             glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
-             height = zu(k) - zu(k_wall)
-             fw     = EXP( -cfw * height / glsf )
-             tkefactor_t(j,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *        &
-                                           ( glsf / glsc )**p23 )
-          ENDDO
-       ENDDO
+!
+!--    RANS mode
+       ELSE
+          IF ( nest_bound_l )  THEN
+             ALLOCATE( tkefactor_l(nzb:nzt+1,nysg:nyng) )
+             tkefactor_l = 1.0_wp
+          ENDIF
+          IF ( nest_bound_r )  THEN
+             ALLOCATE( tkefactor_r(nzb:nzt+1,nysg:nyng) )
+             tkefactor_r = 1.0_wp
+          ENDIF
+          IF ( nest_bound_s )  THEN
+             ALLOCATE( tkefactor_s(nzb:nzt+1,nxlg:nxrg) )
+             tkefactor_s = 1.0_wp
+          ENDIF
+          IF ( nest_bound_n )  THEN
+             ALLOCATE( tkefactor_n(nzb:nzt+1,nxlg:nxrg) )
+             tkefactor_n = 1.0_wp
+          ENDIF
+
+          ALLOCATE( tkefactor_t(nysg:nyng,nxlg:nxrg) )
+          tkefactor_t = 1.0_wp
+
+       ENDIF
       
     END SUBROUTINE pmci_init_tkefactor
@@ -3177,17 +3241,18 @@
 !-- List of array names, which can be coupled.
 !-- In case of 3D please change also the second array for the pointer version
-    IF ( TRIM(name) == "u"  )  p_3d => u
-    IF ( TRIM(name) == "v"  )  p_3d => v
-    IF ( TRIM(name) == "w"  )  p_3d => w
-    IF ( TRIM(name) == "e"  )  p_3d => e
-    IF ( TRIM(name) == "pt" )  p_3d => pt
-    IF ( TRIM(name) == "q"  )  p_3d => q
-    IF ( TRIM(name) == "qc" )  p_3d => qc
-    IF ( TRIM(name) == "qr" )  p_3d => qr
-    IF ( TRIM(name) == "nr" )  p_3d => nr
-    IF ( TRIM(name) == "nc" )  p_3d => nc
-    IF ( TRIM(name) == "s"  )  p_3d => s
-    IF ( TRIM(name) == "nr_part"  )   i_2d => nr_part
-    IF ( TRIM(name) == "part_adr"  )  i_2d => part_adr
+    IF ( TRIM(name) == "u"          )  p_3d => u
+    IF ( TRIM(name) == "v"          )  p_3d => v
+    IF ( TRIM(name) == "w"          )  p_3d => w
+    IF ( TRIM(name) == "e"          )  p_3d => e
+    IF ( TRIM(name) == "pt"         )  p_3d => pt
+    IF ( TRIM(name) == "q"          )  p_3d => q
+    IF ( TRIM(name) == "qc"         )  p_3d => qc
+    IF ( TRIM(name) == "qr"         )  p_3d => qr
+    IF ( TRIM(name) == "nr"         )  p_3d => nr
+    IF ( TRIM(name) == "nc"         )  p_3d => nc
+    IF ( TRIM(name) == "s"          )  p_3d => s
+    IF ( TRIM(name) == "diss"       )  p_3d => diss   
+    IF ( TRIM(name) == "nr_part"    )   i_2d => nr_part
+    IF ( TRIM(name) == "part_adr"   )  i_2d => part_adr
     IF ( INDEX( TRIM(name), "chem_" ) /= 0 )  p_3d => chem_species(n)%conc
 
@@ -3219,15 +3284,16 @@
     ENDIF
 #else
-    IF ( TRIM(name) == "u"  )  p_3d_sec => u_2
-    IF ( TRIM(name) == "v"  )  p_3d_sec => v_2
-    IF ( TRIM(name) == "w"  )  p_3d_sec => w_2
-    IF ( TRIM(name) == "e"  )  p_3d_sec => e_2
-    IF ( TRIM(name) == "pt" )  p_3d_sec => pt_2
-    IF ( TRIM(name) == "q"  )  p_3d_sec => q_2
-    IF ( TRIM(name) == "qc" )  p_3d_sec => qc_2
-    IF ( TRIM(name) == "qr" )  p_3d_sec => qr_2
-    IF ( TRIM(name) == "nr" )  p_3d_sec => nr_2
-    IF ( TRIM(name) == "nc" )  p_3d_sec => nc_2
-    IF ( TRIM(name) == "s"  )  p_3d_sec => s_2
+    IF ( TRIM(name) == "u"    )  p_3d_sec => u_2
+    IF ( TRIM(name) == "v"    )  p_3d_sec => v_2
+    IF ( TRIM(name) == "w"    )  p_3d_sec => w_2
+    IF ( TRIM(name) == "e"    )  p_3d_sec => e_2
+    IF ( TRIM(name) == "pt"   )  p_3d_sec => pt_2
+    IF ( TRIM(name) == "q"    )  p_3d_sec => q_2
+    IF ( TRIM(name) == "qc"   )  p_3d_sec => qc_2
+    IF ( TRIM(name) == "qr"   )  p_3d_sec => qr_2
+    IF ( TRIM(name) == "nr"   )  p_3d_sec => nr_2
+    IF ( TRIM(name) == "nc"   )  p_3d_sec => nc_2
+    IF ( TRIM(name) == "s"    )  p_3d_sec => s_2
+    IF ( TRIM(name) == "diss" )  p_3d_sec => diss_2
     IF ( INDEX( TRIM(name), "chem_" ) /= 0 )  p_3d_sec => spec_conc_2(:,:,:,n)
 
@@ -3358,4 +3424,7 @@
        IF ( .NOT. ALLOCATED( ec ) )  ALLOCATE( ec(0:nzc+1,js:je,is:ie) )
        p_3d => ec
+    ELSEIF ( TRIM( name ) == "diss" )  THEN
+       IF ( .NOT. ALLOCATED( dissc ) )  ALLOCATE( dissc(0:nzc+1,js:je,is:ie) )
+       p_3d => dissc
     ELSEIF ( TRIM( name ) == "pt")  THEN
        IF ( .NOT. ALLOCATED( ptc ) )  ALLOCATE( ptc(0:nzc+1,js:je,is:ie) )
@@ -3379,9 +3448,9 @@
        IF ( .NOT. ALLOCATED( sc ) )  ALLOCATE( sc(0:nzc+1,js:je,is:ie) )
        p_3d => sc
-    ELSEIF (trim(name) == "nr_part") then
-       IF (.not.allocated(nr_partc))  allocate(nr_partc(js:je, is:ie))
+    ELSEIF ( TRIM( name ) == "nr_part") THEN
+       IF ( .NOT. ALLOCATED( nr_partc ) )  ALLOCATE( nr_partc(js:je,is:ie) )
        i_2d => nr_partc
-    ELSEIF (trim(name) == "part_adr") then
-       IF (.not.allocated(part_adrc))  allocate(part_adrc(js:je, is:ie))
+    ELSEIF ( TRIM( name ) == "part_adr") THEN
+       IF ( .NOT. ALLOCATED( part_adrc ) )  ALLOCATE( part_adrc(js:je,is:ie) )
        i_2d => part_adrc
     ELSEIF ( TRIM( name(1:5) ) == "chem_" )  THEN
@@ -3484,9 +3553,16 @@
                                    r2yo, r1zw, r2zw, 'w' )
 
-       IF ( .NOT. constant_diffusion )  THEN
-          CALL pmci_interp_tril_all ( e,  ec,  ico, jco, kco, r1xo, r2xo,      &
-                                      r1yo, r2yo, r1zo, r2zo, 'e' )
-       ENDIF
-       
+       IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.          &
+            (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.           &
+               .NOT. constant_diffusion ) )  THEN
+          CALL pmci_interp_tril_all ( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo, &
+                                      r2yo, r1zo, r2zo, 'e' )
+       ENDIF
+
+       IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
+          CALL pmci_interp_tril_all ( diss,  dissc,  ico, jco, kco, r1xo, r2xo,&
+                                      r1yo, r2yo, r1zo, r2zo, 's' )
+       ENDIF
+
        IF ( .NOT. neutral )  THEN
           CALL pmci_interp_tril_all ( pt, ptc, ico, jco, kco, r1xo, r2xo,      &
@@ -4207,5 +4283,7 @@
                                        nzt_topo_nestbc_l, 'l', 'w' )
 
-             IF ( .NOT. constant_diffusion )  THEN
+             IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.   &
+                  (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.    &
+                     .NOT. constant_diffusion ) )  THEN
                 CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
                                           r1yo, r2yo, r1zo, r2zo,              &
@@ -4214,5 +4292,13 @@
                                           nzt_topo_nestbc_l, 'l', 'e' )
              ENDIF
-             
+
+             IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
+                CALL pmci_interp_tril_lr( diss,  dissc,  ico, jco, kco, r1xo,  &
+                                          r2xo, r1yo, r2yo, r1zo, r2zo,        &
+                                          logc_w_l, logc_ratio_w_l,            &
+                                          logc_kbounds_w_l,                    &
+                                          nzt_topo_nestbc_l, 'l', 's' )
+             ENDIF
+
              IF ( .NOT. neutral )  THEN
                 CALL pmci_interp_tril_lr( pt, ptc, ico, jco, kco, r1xo, r2xo,  &
@@ -4304,5 +4390,7 @@
                                        nzt_topo_nestbc_r, 'r', 'w' )
 
-             IF ( .NOT. constant_diffusion )  THEN
+             IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.   &
+                  (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.    &
+                     .NOT. constant_diffusion ) )  THEN
                 CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
                                           r1yo,r2yo, r1zo, r2zo,               &
@@ -4310,4 +4398,14 @@
                                           logc_kbounds_w_r,                    &
                                           nzt_topo_nestbc_r, 'r', 'e' )
+
+             ENDIF
+
+             IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
+                CALL pmci_interp_tril_lr( diss,  dissc,  ico, jco, kco, r1xo,  &
+                                          r2xo, r1yo,r2yo, r1zo, r2zo,         &
+                                          logc_w_r, logc_ratio_w_r,            &
+                                          logc_kbounds_w_r,                    &
+                                          nzt_topo_nestbc_r, 'r', 's' )
+
              ENDIF
 
@@ -4382,5 +4480,4 @@
                 ENDDO
              ENDIF
-
           ENDIF
 !
@@ -4406,5 +4503,7 @@
                                        nzt_topo_nestbc_s, 's','w' )
 
-             IF ( .NOT. constant_diffusion )  THEN
+             IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.   &
+                  (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.    &
+                     .NOT. constant_diffusion ) )  THEN
                 CALL pmci_interp_tril_sn( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
                                           r1yo, r2yo, r1zo, r2zo,              &
@@ -4412,6 +4511,16 @@
                                           logc_kbounds_w_s,                    &
                                           nzt_topo_nestbc_s, 's', 'e' )
+
              ENDIF
-          
+
+             IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
+                CALL pmci_interp_tril_sn( diss, dissc,  ico, jco, kco, r1xo,   &
+                                          r2xo, r1yo, r2yo, r1zo, r2zo,        &
+                                          logc_w_s, logc_ratio_w_s,            &
+                                          logc_kbounds_w_s,                    &
+                                          nzt_topo_nestbc_s, 's', 's' )
+
+             ENDIF
+
              IF ( .NOT. neutral )  THEN
                 CALL pmci_interp_tril_sn( pt, ptc, ico, jco, kco, r1xo, r2xo,  &
@@ -4482,5 +4591,4 @@
                 ENDDO
              ENDIF
-
           ENDIF
 !
@@ -4505,5 +4613,8 @@
                                        logc_kbounds_w_n,                       &
                                        nzt_topo_nestbc_n, 'n', 'w' )
-             IF ( .NOT. constant_diffusion )  THEN
+
+             IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.   & 
+                  (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.    &
+                     .NOT. constant_diffusion ) )  THEN
                 CALL pmci_interp_tril_sn( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
                                           r1yo, r2yo, r1zo, r2zo,              &
@@ -4511,6 +4622,16 @@
                                           logc_kbounds_w_n,                    &
                                           nzt_topo_nestbc_n, 'n', 'e' )
+
              ENDIF
-          
+
+             IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
+                CALL pmci_interp_tril_sn( diss, dissc,  ico, jco, kco, r1xo,   &
+                                          r2xo, r1yo, r2yo, r1zo, r2zo,        &
+                                          logc_w_n, logc_ratio_w_n,            &
+                                          logc_kbounds_w_n,                    &
+                                          nzt_topo_nestbc_n, 'n', 's' )
+
+             ENDIF
+
              IF ( .NOT. neutral )  THEN
                 CALL pmci_interp_tril_sn( pt, ptc, ico, jco, kco, r1xo, r2xo,  &
@@ -4581,5 +4702,4 @@
                 ENDDO
              ENDIF
-             
           ENDIF
 
@@ -4593,8 +4713,17 @@
        CALL pmci_interp_tril_t( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo,      &
                                 r2yo, r1zw, r2zw, 'w' )
-       IF ( .NOT. constant_diffusion )  THEN
+
+       IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.         &
+            (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.          &
+               .NOT. constant_diffusion ) )  THEN
           CALL pmci_interp_tril_t( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo,   &
                                    r2yo, r1zo, r2zo, 'e' )
        ENDIF
+
+       IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
+          CALL pmci_interp_tril_t( diss, dissc, ico, jco, kco, r1xo, r2xo,     &
+                                   r1yo, r2yo, r1zo, r2zo, 's' )
+       ENDIF
+
        IF ( .NOT. neutral )  THEN
           CALL pmci_interp_tril_t( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo,   &
@@ -4644,5 +4773,5 @@
           ENDDO
        ENDIF
-    
+
    END SUBROUTINE pmci_interpolation
 
@@ -4666,4 +4795,18 @@
       CALL pmci_anterp_tophat( w,  wc,  kctw, iflo, ifuo, jflo, jfuo, kflw,    &
                                kfuw, ijfc_s, kfc_w, 'w' )
+!
+!--   Anterpolation of TKE and dissipation rate if parent and child are in 
+!--   RANS mode.
+      IF ( rans_mode_parent  .AND.  rans_mode )  THEN
+         CALL pmci_anterp_tophat( e, ec, kctu, iflo, ifuo, jflo, jfuo, kflo,   &
+                                  kfuo, ijfc_s, kfc_s, 'e' )
+!
+!--      Anterpolation of dissipation rate only if TKE-e closure is applied.
+         IF ( rans_tke_e )  THEN
+            CALL pmci_anterp_tophat( diss, dissc, kctu, iflo, ifuo, jflo, jfuo,&
+                                     kflo, kfuo, ijfc_s, kfc_s, 'diss' )
+         ENDIF
+
+      ENDIF
 
       IF ( .NOT. neutral )  THEN
@@ -5569,3 +5712,4 @@
  END SUBROUTINE pmci_boundary_conds
 
+
 END MODULE pmc_interface
Index: palm/trunk/SOURCE/pmc_mpi_wrapper_mod.f90
===================================================================
--- palm/trunk/SOURCE/pmc_mpi_wrapper_mod.f90	(revision 2934)
+++ palm/trunk/SOURCE/pmc_mpi_wrapper_mod.f90	(revision 2938)
@@ -26,4 +26,7 @@
 ! -----------------
 ! $Id$
+! Extent interface by logical buffer 
+! 
+! 2936 2018-03-27 14:49:27Z suehring
 ! Bugfix: wrong placement of include 'mpif.h' corrected
 ! 
@@ -103,4 +106,5 @@
 
     INTERFACE pmc_recv_from_parent
+       MODULE PROCEDURE pmc_recv_from_parent_logical
        MODULE PROCEDURE pmc_recv_from_parent_integer
        MODULE PROCEDURE pmc_recv_from_parent_real_r1
@@ -110,4 +114,5 @@
 
     INTERFACE pmc_send_to_child
+       MODULE PROCEDURE pmc_send_to_child_logical
        MODULE PROCEDURE pmc_send_to_child_integer
        MODULE PROCEDURE pmc_send_to_child_real_r1
@@ -149,4 +154,41 @@
 
 
+ SUBROUTINE pmc_recv_from_parent_logical( buf, n, parent_rank, tag, ierr )
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(IN)                ::  n            !<
+    INTEGER, INTENT(IN)                ::  parent_rank  !<
+    INTEGER, INTENT(IN)                ::  tag          !<
+    INTEGER, INTENT(OUT)               ::  ierr         !<
+
+    LOGICAL, INTENT(OUT)               ::  buf          !<
+    
+    ierr = 0
+    CALL MPI_RECV( buf, n, MPI_LOGICAL, parent_rank, tag, m_to_parent_comm,    &
+                   MPI_STATUS_IGNORE, ierr )
+
+ END SUBROUTINE pmc_recv_from_parent_logical
+
+ SUBROUTINE pmc_send_to_child_logical( child_id, buf, n, child_rank, tag,      &
+                                       ierr )
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(IN)               ::  child_id     !<
+    INTEGER, INTENT(IN)               ::  n            !<
+    INTEGER, INTENT(IN)               ::  child_rank   !<
+    INTEGER, INTENT(IN)               ::  tag          !<
+    INTEGER, INTENT(OUT)              ::  ierr         !<
+
+    LOGICAL, INTENT(IN)               ::  buf          !<
+   
+    ierr = 0
+    CALL MPI_SEND( buf, n, MPI_LOGICAL, child_rank, tag,                       &
+                   m_to_child_comm(child_id), ierr )
+
+ END SUBROUTINE pmc_send_to_child_logical
+
+
  SUBROUTINE pmc_send_to_parent_integer( buf, n, parent_rank, tag, ierr )
 
@@ -167,5 +209,4 @@
 
 
-
  SUBROUTINE pmc_recv_from_parent_integer( buf, n, parent_rank, tag, ierr )
 
@@ -312,5 +353,4 @@
 
  END SUBROUTINE pmc_recv_from_parent_real_r3
-
 
 
Index: palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90
===================================================================
--- palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90	(revision 2934)
+++ palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90	(revision 2938)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Apply turbulence generator at all non-cyclic lateral boundaries in case of 
+! realistic Inifor large-scale forcing or RANS-LES nesting 
+! 
+! 2936 2018-03-27 14:49:27Z suehring
 ! variable named found has been introduced for checking if restart data was found, 
 ! reading of restart strings has been moved completely to read_restart_data_mod, 
@@ -101,6 +105,5 @@
 
     USE control_parameters,                                                    &
-        ONLY:  initializing_actions, message_string,                           &
-               syn_turb_gen
+        ONLY:  initializing_actions, message_string, syn_turb_gen
 
     USE cpulog,                                                                &
@@ -108,5 +111,5 @@
 
     USE indices,                                                               &
-        ONLY:  nbgp, nzb, nzt, nyng, nysg
+        ONLY:  nbgp, nzb, nzt, nxl, nxlg, nxr, nxrg, nys, nyn, nyng, nysg
 
     USE kinds
@@ -117,5 +120,5 @@
 
     USE pegrid,                                                                &
-        ONLY:  comm1dx, comm1dy, comm2d, id_inflow, ierr, myidx, pdims
+        ONLY:  comm1dx, comm1dy, comm2d, ierr, myidx, myidy, pdims
 
     USE transpose_indices,                                                     &
@@ -132,19 +135,34 @@
     LOGICAL :: use_syn_turb_gen = .FALSE.           !< switch to use synthetic turbulence generator
 
-    INTEGER(iwp) :: stg_type_yz        !< MPI type for full z range
-    INTEGER(iwp) :: stg_type_yz_small  !< MPI type for small z range
-    INTEGER(iwp) :: merg               !< maximum length scale (in gp)
-    INTEGER(iwp) :: mergp              !< merg + nbgp
+    INTEGER(iwp) ::  id_stg_left        !< left lateral boundary core id in case of turbulence generator  
+    INTEGER(iwp) ::  id_stg_north       !< north lateral boundary core id in case of turbulence generator  
+    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) ::  stg_type_xz        !< MPI type for full z range
+    INTEGER(iwp) ::  stg_type_xz_small  !< MPI type for small z range
+    INTEGER(iwp) ::  stg_type_yz        !< MPI type for full z range
+    INTEGER(iwp) ::  stg_type_yz_small  !< MPI type for small z range
+    INTEGER(iwp) ::  merg               !< maximum length scale (in gp)
+    INTEGER(iwp) ::  mergp              !< merg + nbgp
+    INTEGER(iwp) ::  nzb_x_stg          !< 
+    INTEGER(iwp) ::  nzt_x_stg          !< 
+    INTEGER(iwp) ::  nzb_y_stg          !< 
+    INTEGER(iwp) ::  nzt_y_stg          !< 
 
     REAL(wp) :: mc_factor    !< mass flux correction factor
 
-    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displs      !< displacement for MPI_GATHERV
-    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: recv_count  !< receive count for MPI_GATHERV
-    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuy         !< length scale of u in y direction (in gp)
-    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuz         !< length scale of u in z direction (in gp)
-    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvy         !< length scale of v in y direction (in gp)
-    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvz         !< length scale of v in z direction (in gp)
-    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwy         !< length scale of w in y direction (in gp)
-    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwz         !< length scale of w in z direction (in gp)
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displs_xz      !< displacement for MPI_GATHERV
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: recv_count_xz  !< receive count for MPI_GATHERV
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displs_yz      !< displacement for MPI_GATHERV
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: recv_count_yz  !< receive count for MPI_GATHERV
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nux            !< length scale of u in x direction (in gp)
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuy            !< length scale of u in y direction (in gp)
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuz            !< length scale of u in z direction (in gp)
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvx            !< length scale of v in x direction (in gp)
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvy            !< length scale of v in y direction (in gp)
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvz            !< length scale of v in z direction (in gp)
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwx            !< length scale of w in x direction (in gp)
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwy            !< length scale of w in y direction (in gp)
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwz            !< length scale of w in z direction (in gp)
 
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: seed        !< seed of random number for rn-generator
@@ -160,16 +178,25 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE :: tw              !< Lagrangian time scale of w
 
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bux           !< filter function for u in x direction
     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buy           !< filter function for u in y direction
     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buz           !< filter function for u in z direction
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvx           !< filter function for v in x direction
     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvy           !< filter function for v in y direction
     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvz           !< filter function for v in z direction
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwx           !< filter function for w in y direction
     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwy           !< filter function for w in y direction
     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwz           !< filter function for w in z direction
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu            !< velocity seed for u
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo           !< velocity seed for u with new random number
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv            !< velocity seed for v
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo           !< velocity seed for v with new random number
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw            !< velocity seed for w
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo           !< velocity seed for w with new random number
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu_xz         !< velocity seed for u at xz plane
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo_xz        !< velocity seed for u at xz plane with new random number
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu_yz         !< velocity seed for u at yz plane
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo_yz        !< velocity seed for u at yz plane with new random number
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv_xz         !< velocity seed for v at xz plane
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo_xz        !< velocity seed for v at xz plane with new random number
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv_yz         !< velocity seed for v at yz plane
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo_yz        !< velocity seed for v at yz plane with new random number
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw_xz         !< velocity seed for w at xz plane
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo_xz        !< velocity seed for w at xz plane with new random number
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw_yz         !< velocity seed for w at yz plane
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo_yz        !< velocity seed for w at yz plane with new random number
 
 
@@ -188,5 +215,10 @@
 
 !
-!-- Generate velocity seeds
+!-- Generate velocity seeds at south and north domain boundary
+    INTERFACE stg_generate_seed_xz
+       MODULE PROCEDURE stg_generate_seed_xz
+    END INTERFACE stg_generate_seed_xz
+!
+!-- Generate velocity seeds at left and/or right domain boundary
     INTERFACE stg_generate_seed_yz
        MODULE PROCEDURE stg_generate_seed_yz
@@ -240,5 +272,6 @@
 !
 !-- Public variables
-    PUBLIC  use_syn_turb_gen
+    PUBLIC  id_stg_left, id_stg_north, id_stg_right, id_stg_south,             &
+            use_syn_turb_gen
 
 
@@ -255,33 +288,55 @@
 
     USE control_parameters,                                                    &
-        ONLY:  bc_lr, bc_ns, turbulent_inflow
+        ONLY:  bc_lr, bc_ns, forcing, nest_domain, rans_mode, turbulent_inflow
+
+    USE pmc_interface,                                                         &
+        ONLY : rans_mode_parent
 
 
     IMPLICIT NONE
 
+    IF ( .NOT. use_syn_turb_gen  .AND.  .NOT. rans_mode  .AND.  forcing )  THEN
+       message_string = 'Synthetic turbulence generator has to be applied ' // &
+                        'when forcing is used and model operates in LES mode.'
+       CALL message( 'stg_check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
+    ENDIF
+
+    IF ( .NOT. use_syn_turb_gen  .AND.  nest_domain                            &
+         .AND. rans_mode_parent  .AND.  .NOT. rans_mode )  THEN
+       message_string = 'Synthetic turbulence generator has to be applied ' // &
+                        'when nesting is applied and parent operates in '  //  &
+                        'RANS-mode but current child in LES mode.'
+       CALL message( 'stg_check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
+    ENDIF
+
     IF ( use_syn_turb_gen )  THEN
 
-       IF ( INDEX( initializing_actions, 'set_constant_profiles' ) == 0  .AND. &
-            INDEX( initializing_actions, 'read_restart_data' ) == 0 )  THEN
-          message_string = 'Using synthetic turbulence generator ' //          &
-                           'requires &initializing_actions = '            //   &
-                           '"set_constant_profiles" or "read_restart_data"'
-          CALL message( 'stg_check_parameters', 'PA0015', 1, 2, 0, 6, 0 )
+       IF ( .NOT. forcing  .AND.  .NOT. nest_domain )  THEN 
+
+          IF ( INDEX( initializing_actions, 'set_constant_profiles' ) == 0     &
+        .AND.  INDEX( initializing_actions, 'read_restart_data' ) == 0 )  THEN
+             message_string = 'Using synthetic turbulence generator ' //       &
+                              'requires &initializing_actions = '         //   &
+                              '"set_constant_profiles" or "read_restart_data"'
+             CALL message( 'stg_check_parameters', 'PA0015', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( bc_lr /= 'dirichlet/radiation' )  THEN
+             message_string = 'Using synthetic turbulence generator ' //       &
+                              'requires &bc_lr = "dirichlet/radiation"'
+             CALL message( 'stg_check_parameters', 'PA0035', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF ( bc_ns /= 'cyclic' )  THEN
+             message_string = 'Using synthetic turbulence generator ' //       &
+                              'requires &bc_ns = "cyclic"'
+             CALL message( 'stg_check_parameters', 'PA0037', 1, 2, 0, 6, 0 )
+          ENDIF
+
        ENDIF
 
-       IF ( bc_lr /= 'dirichlet/radiation' )  THEN
-          message_string = 'Using synthetic turbulence generator ' //          &
-                           'requires &bc_lr = "dirichlet/radiation"'
-          CALL message( 'stg_check_parameters', 'PA0035', 1, 2, 0, 6, 0 )
-       ENDIF
-       IF ( bc_ns /= 'cyclic' )  THEN
-          message_string = 'Using synthetic turbulence generator ' //          &
-                           'requires &bc_ns = "cyclic"'
-          CALL message( 'stg_check_parameters', 'PA0037', 1, 2, 0, 6, 0 )
-       ENDIF
        IF ( turbulent_inflow )  THEN
           message_string = 'Using synthetic turbulence generator ' //          &
-                           'in combination &with turbulent_inflow = .T. ' //   &
-                           'is not allowed'
+                           'in combination &with turbulent_inflow = .T. '//    &
+                              'is not allowed'
           CALL message( 'stg_check_parameters', 'PA0039', 1, 2, 0, 6, 0 )
        ENDIF
@@ -330,14 +385,19 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  u_init, v_init
+        ONLY:  ddzw, u_init, v_init, zu
 
     USE control_parameters,                                                    &
-        ONLY:  coupling_char, dz, e_init
+        ONLY:  coupling_char, dz, e_init, forcing, nest_domain
 
     USE grid_variables,                                                        &
-        ONLY:  ddy
+        ONLY:  ddx, ddy
+
+    USE indices,                                                               &
+        ONLY:  nz
 
 
     IMPLICIT NONE
+
+    LOGICAL ::  file_stg_exist = .FALSE. !< flag indication whether parameter file for Reynolds stress and length scales exist
 
 #if defined( __parallel )
@@ -346,8 +406,10 @@
 #endif
 
+    INTEGER(iwp) :: i                        !> grid index in x-direction
     INTEGER(iwp) :: j                        !> loop index
     INTEGER(iwp) :: k                        !< index
     INTEGER(iwp) :: newtype                  !< dummy MPI type
     INTEGER(iwp) :: realsize                 !< size of REAL variables
+    INTEGER(iwp) :: nnz                      !< increment used to determine processor decomposition of z-axis along x and y direction  
     INTEGER(iwp) :: nseed                    !< dimension of random number seed
     INTEGER(iwp) :: startseed = 1234567890   !< start seed for random number generator
@@ -362,4 +424,5 @@
     REAL(wp) :: d21     !< vertical interpolation for a21
     REAL(wp) :: d22     !< vertical interpolation for a22
+    REAL(wp) :: dum_exp !< dummy variable used for exponential vertical decrease of turbulent length scales
     REAL(wp) :: luy     !< length scale for u in y direction
     REAL(wp) :: luz     !< length scale for u in z direction
@@ -370,4 +433,6 @@
     REAL(wp) :: zz      !< height
 
+    REAL(wp) :: length_scale_surface, r_ii_0, time_scale, length_scale_z
+
     REAL(wp),DIMENSION(nzb:nzt+1) :: r11  !< Reynolds parameter
     REAL(wp),DIMENSION(nzb:nzt+1) :: r21  !< Reynolds parameter
@@ -389,9 +454,28 @@
     ALLOCATE ( a11(nzb:nzt+1), a21(nzb:nzt+1), a22(nzb:nzt+1),                 &
                a31(nzb:nzt+1), a32(nzb:nzt+1), a33(nzb:nzt+1),                 &
-               nuy(nzb:nzt+1), nuz(nzb:nzt+1), nvy(nzb:nzt+1),                 &
-               nvz(nzb:nzt+1), nwy(nzb:nzt+1), nwz(nzb:nzt+1),                 &
+               nux(nzb:nzt+1), nuy(nzb:nzt+1), nuz(nzb:nzt+1),                 &
+               nvx(nzb:nzt+1), nvy(nzb:nzt+1), nvz(nzb:nzt+1),                 &
+               nwx(nzb:nzt+1), nwy(nzb:nzt+1), nwz(nzb:nzt+1),                 &
                tu(nzb:nzt+1),  tv(nzb:nzt+1),  tw(nzb:nzt+1)   )
 
 #if defined( __parallel )
+!
+!-- Determine processor decomposition of z-axis along x- and y-direction
+    nnz = nz / pdims(1)
+    nzb_x_stg = 1 + myidx * nnz
+    nzt_x_stg = ( myidx + 1 ) * nnz
+
+    IF ( MOD( nz , pdims(1) ) /= 0  .AND.  myidx == id_stg_right )             &  
+       nzt_x_stg = myidx * nnz + MOD( nz , pdims(1) )
+
+    IF ( forcing  .OR.  nest_domain )  THEN
+       nnz = nz / pdims(2)
+       nzb_y_stg = 1 + myidy * nnz
+       nzt_y_stg = ( myidy + 1 ) * nnz
+
+       IF ( MOD( nz , pdims(2) ) /= 0  .AND.  myidy == id_stg_north )          &  
+          nzt_y_stg = myidy * nnz + MOD( nz , pdims(2) )
+    ENDIF
+
 !
 !-- Define MPI type used in stg_generate_seed_yz to gather vertical splitted
@@ -399,6 +483,8 @@
     CALL MPI_TYPE_SIZE( MPI_REAL, realsize, ierr )
     extent = 1 * realsize
-
-    ! stg_type_yz: yz-slice with vertical bounds nzb:nzt+1
+!
+!-- Set-up MPI datatyp to involve all cores for turbulence generation at yz 
+!-- layer 
+!-- stg_type_yz: yz-slice with vertical bounds nzb:nzt+1
     CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt-nzb+2,nyng-nysg+1],                 &
             [1,nyng-nysg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
@@ -407,6 +493,6 @@
     CALL MPI_TYPE_FREE( newtype, ierr )
 
-    ! stg_type_yz_small: yz-slice with vertical bounds nzb_x:nzt_x+1
-    CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_x-nzb_x+2,nyng-nysg+1],             &
+    ! stg_type_yz_small: yz-slice with vertical bounds nzb_x_stg:nzt_x_stg+1
+    CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_x_stg-nzb_x_stg+2,nyng-nysg+1],             &
             [1,nyng-nysg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
     CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_yz_small, ierr )
@@ -415,14 +501,43 @@
 
     ! receive count and displacement for MPI_GATHERV in stg_generate_seed_yz
-    ALLOCATE( recv_count(pdims(1)), displs(pdims(1)) )
-
-    recv_count           = nzt_x-nzb_x + 1
-    recv_count(pdims(1)) = recv_count(pdims(1)) + 1
+    ALLOCATE( recv_count_yz(pdims(1)), displs_yz(pdims(1)) )
+
+    recv_count_yz           = nzt_x_stg-nzb_x_stg + 1
+    recv_count_yz(pdims(1)) = recv_count_yz(pdims(1)) + 1
 
     DO  j = 1, pdims(1)
-       displs(j) = 0 + (nzt_x-nzb_x+1) * (j-1)
-    ENDDO
+       displs_yz(j) = 0 + (nzt_x_stg-nzb_x_stg+1) * (j-1)
+    ENDDO
+!
+!-- Set-up MPI datatyp to involve all cores for turbulence generation at xz 
+!-- layer 
+!-- stg_type_xz: xz-slice with vertical bounds nzb:nzt+1
+    IF ( forcing  .OR.  nest_domain)  THEN
+       CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt-nzb+2,nxrg-nxlg+1],                 &
+               [1,nxrg-nxlg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
+       CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_xz, ierr )
+       CALL MPI_TYPE_COMMIT( stg_type_xz, ierr )
+       CALL MPI_TYPE_FREE( newtype, ierr )
+
+       ! stg_type_yz_small: xz-slice with vertical bounds nzb_x_stg:nzt_x_stg+1
+       CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_y_stg-nzb_y_stg+2,nxrg-nxlg+1],             &
+               [1,nxrg-nxlg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
+       CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_xz_small, ierr )
+       CALL MPI_TYPE_COMMIT( stg_type_xz_small, ierr )
+       CALL MPI_TYPE_FREE( newtype, ierr )
+
+       ! receive count and displacement for MPI_GATHERV in stg_generate_seed_yz
+       ALLOCATE( recv_count_xz(pdims(2)), displs_xz(pdims(2)) )
+
+       recv_count_xz           = nzt_y_stg-nzb_y_stg + 1
+       recv_count_xz(pdims(2)) = recv_count_xz(pdims(2)) + 1
+
+       DO  j = 1, pdims(2)
+          displs_xz(j) = 0 + (nzt_y_stg-nzb_y_stg+1) * (j-1)
+       ENDDO
+
+    ENDIF
+
 #endif
-
 !
 !-- Define seed of random number
@@ -435,5 +550,4 @@
     CALL RANDOM_SEED(put=seed)
 
-!
 !-- Read inflow profile
 !-- Height levels of profiles in input profile is as follows:
@@ -441,43 +555,118 @@
 !-- zw: lwy, lwz, tw, r31, r32, r33, d3
 !-- WARNING: zz is not used at the moment
-    OPEN( 90, FILE='STG_PROFILES'//TRIM( coupling_char ), STATUS='OLD',        &
-                   FORM='FORMATTED')
-
-    ! Skip header
-    READ( 90, * )
-
-    DO  k = nzb, nzt+1
-       READ( 90, * ) zz, luy, luz, tu(k), lvy, lvz, tv(k), lwy, lwz, tw(k),    &
-                     r11(k), r21(k), r22(k), r31(k), r32(k), r33(k),           &
-                     d1, d2, d3, d5
-
-!
-!--    Convert length scales from meter to number of grid points
-       nuy(k) = INT( luy * ddy )
-       nuz(k) = INT( luz / dz  )
-       nvy(k) = INT( lvy * ddy )
-       nvz(k) = INT( lvz / dz  )
-       nwy(k) = INT( lwy * ddy )
-       nwz(k) = INT( lwz / dz  )
-
-!
-!--    Save Mean inflow profiles
-       IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
-          mean_inflow_profiles(k,1) = d1
-          mean_inflow_profiles(k,2) = d2
-         !  mean_inflow_profiles(k,4) = d4
-          mean_inflow_profiles(k,5) = d5
-       ENDIF
-    ENDDO
-
-    CLOSE( 90 )
-
-!
-!--    Assign initial profiles
-    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) )
-
+    INQUIRE( FILE = 'STG_PROFILES' // TRIM( coupling_char ),                   &
+             EXIST = file_stg_exist  )
+
+    IF ( file_stg_exist )  THEN
+
+       OPEN( 90, FILE='STG_PROFILES'//TRIM( coupling_char ), STATUS='OLD',     &
+                      FORM='FORMATTED')
+!
+!--    Skip header
+       READ( 90, * )
+
+       DO  k = nzb, nzt+1
+          READ( 90, * ) zz, luy, luz, tu(k), lvy, lvz, tv(k), lwy, lwz, tw(k), &
+                        r11(k), r21(k), r22(k), r31(k), r32(k), r33(k),        &
+                        d1, d2, d3, d5
+
+!
+!--       Convert length scales from meter to number of grid points
+          nuy(k) = INT( luy * ddy )
+          nuz(k) = INT( luz / dz  )
+          nvy(k) = INT( lvy * ddy )
+          nvz(k) = INT( lvz / dz  )
+          nwy(k) = INT( lwy * ddy )
+          nwz(k) = INT( lwz / dz  )
+!
+!--       Workaround, assume isotropic turbulence
+          nwx(k) = nwy(k)
+          nvx(k) = nvy(k)
+          nux(k) = nuy(k)
+!
+!--       Save Mean inflow profiles
+          IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
+             mean_inflow_profiles(k,1) = d1
+             mean_inflow_profiles(k,2) = d2
+            !  mean_inflow_profiles(k,4) = d4
+             mean_inflow_profiles(k,5) = d5
+          ENDIF
+       ENDDO
+
+       CLOSE( 90 )
+   
+    ELSE
+!
+!--    Set-up defaul length scales. Assume exponentially decreasing length 
+!--    scales and isotropic turbulence. 
+!--    Typical length (time) scales of 100 m (s) should be a goog compromise 
+!--    between all stratrifications. Near-surface variances are fixed to 
+!--    0.1 m2/s2, vertical fluxes are one order of magnitude smaller.
+!--    Vertical fluxes  
+       length_scale_surface = 100.0_wp
+       r_ii_0               = 0.1_wp
+       time_scale           = 100.0_wp
+       DO  k = nzb+1, nzt+1
+          dum_exp        = MERGE( -zu(k) / ( 0.3* zu(nzt) ),                   &
+                                  0.0_wp,                                      &
+                                  zu(k) > 0.3 * zu(nzt)                        &
+                                )
+          length_scale_z = length_scale_surface * EXP( dum_exp )
+
+          nux(k) = MAX( INT( length_scale_z * ddx     ), 1 )
+          nuy(k) = MAX( INT( length_scale_z * ddy     ), 1 )
+          nuz(k) = MAX( INT( length_scale_z * ddzw(k) ), 1 )
+          nvx(k) = MAX( INT( length_scale_z * ddx     ), 1 )
+          nvy(k) = MAX( INT( length_scale_z * ddy     ), 1 )
+          nvz(k) = MAX( INT( length_scale_z * ddzw(k) ), 1 )
+          nwx(k) = MAX( INT( length_scale_z * ddx     ), 1 )
+          nwy(k) = MAX( INT( length_scale_z * ddy     ), 1 )
+          nwz(k) = MAX( INT( length_scale_z * ddzw(k) ), 1 )
+
+          r11(k) = r_ii_0 * EXP( dum_exp )
+          r22(k) = r_ii_0 * EXP( dum_exp )
+          r33(k) = r_ii_0 * EXP( dum_exp )
+
+          r21(k) = 0.1_wp * r_ii_0 * EXP( dum_exp )
+          r31(k) = 0.1_wp * r_ii_0 * EXP( dum_exp )
+          r32(k) = 0.1_wp * r_ii_0 * EXP( dum_exp )
+
+          tu(k)  = time_scale
+          tv(k)  = time_scale
+          tw(k)  = time_scale
+
+       ENDDO
+       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)
+
+       r11(nzb) = r11(nzb+1)
+       r22(nzb) = r22(nzb+1)
+       r33(nzb) = r33(nzb+1)
+
+       r21(nzb) = r11(nzb+1)
+       r31(nzb) = r31(nzb+1)
+       r32(nzb) = r32(nzb+1)
+
+       tu(nzb)  = time_scale
+       tv(nzb)  = time_scale
+       tw(nzb)  = time_scale
+
+    ENDIF
+
+!
+!-- Assign initial profiles
+    IF ( .NOT. forcing  .AND.  .NOT.  nest_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
 !
 !-- Calculate coefficient matrix from Reynolds stress tensor (Lund rotation)
@@ -534,5 +723,4 @@
 
     ENDDO
-
 !
 !-- Define the size of the filter functions and allocate them.
@@ -541,7 +729,7 @@
     ! arrays must be large enough to cover the largest length scale
     DO  k = nzb, nzt+1
-       j = MAX( ABS(nuy(k)), ABS(nuz(k)), &
-                ABS(nvy(k)), ABS(nvz(k)), &
-                ABS(nwy(k)), ABS(nwz(k))  )
+       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 > merg )  merg = j
     ENDDO
@@ -550,27 +738,49 @@
     mergp = merg + nbgp
 
-    ALLOCATE ( buy(-merg:merg,nzb:nzt+1), buz(-merg:merg,nzb:nzt+1), &
-               bvy(-merg:merg,nzb:nzt+1), bvz(-merg:merg,nzb:nzt+1), &
-               bwy(-merg:merg,nzb:nzt+1), bwz(-merg:merg,nzb:nzt+1)  )
-
-!
-!-- Allocate velocity seeds
-    ALLOCATE ( fu( nzb:nzt+1,nysg:nyng), fuo(nzb:nzt+1,nysg:nyng), &
-               fv( nzb:nzt+1,nysg:nyng), fvo(nzb:nzt+1,nysg:nyng), &
-               fw( nzb:nzt+1,nysg:nyng), fwo(nzb:nzt+1,nysg:nyng)  )
-
-    fu  = 0._wp
-    fuo = 0._wp
-    fv  = 0._wp
-    fvo = 0._wp
-    fw  = 0._wp
-    fwo = 0._wp
+    ALLOCATE ( bux(-merg:merg,nzb:nzt+1),                                      &
+               buy(-merg:merg,nzb:nzt+1),                                      &
+               buz(-merg:merg,nzb:nzt+1),                                      &
+               bvx(-merg:merg,nzb:nzt+1),                                      &
+               bvy(-merg:merg,nzb:nzt+1),                                      &
+               bvz(-merg:merg,nzb:nzt+1),                                      &
+               bwx(-merg:merg,nzb:nzt+1),                                      &
+               bwy(-merg:merg,nzb:nzt+1),                                      &
+               bwz(-merg:merg,nzb:nzt+1)  )
+
+!
+!-- Allocate velocity seeds for turbulence at xz-layer
+    ALLOCATE ( fu_xz( nzb:nzt+1,nxlg:nxrg), fuo_xz(nzb:nzt+1,nxlg:nxrg),       &
+               fv_xz( nzb:nzt+1,nxlg:nxrg), fvo_xz(nzb:nzt+1,nxlg:nxrg),       &
+               fw_xz( nzb:nzt+1,nxlg:nxrg), fwo_xz(nzb:nzt+1,nxlg:nxrg)  )
+
+!
+!-- Allocate velocity seeds for turbulence at yz-layer
+    ALLOCATE ( fu_yz( nzb:nzt+1,nysg:nyng), fuo_yz(nzb:nzt+1,nysg:nyng),       &
+               fv_yz( nzb:nzt+1,nysg:nyng), fvo_yz(nzb:nzt+1,nysg:nyng),       &
+               fw_yz( nzb:nzt+1,nysg:nyng), fwo_yz(nzb:nzt+1,nysg:nyng)  )
+
+    fu_xz  = 0.0_wp
+    fuo_xz = 0.0_wp
+    fv_xz  = 0.0_wp
+    fvo_xz = 0.0_wp
+    fw_xz  = 0.0_wp
+    fwo_xz = 0.0_wp
+
+    fu_yz  = 0.0_wp
+    fuo_yz = 0.0_wp
+    fv_yz  = 0.0_wp
+    fvo_yz = 0.0_wp
+    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
@@ -581,32 +791,35 @@
 
 !
-!--    In case of restart, calculate velocity seeds fu, fv, fw from former
-!      time step.
-!      Bug: fu, fv, fw are different in those heights where a11, a22, a33
-!           are 0 compared to the prerun. This is mostly for k=nzt+1.
+!-- In case of restart, calculate velocity seeds fu, fv, fw from former
+!   time step.
+!   Bug: fu, fv, fw are different in those heights where a11, a22, a33
+!        are 0 compared to the prerun. This is mostly for k=nzt+1.
     IF ( TRIM( initializing_actions ) == 'read_restart_data' )  THEN
-       IF ( myidx == id_inflow )  THEN
+       IF ( myidx == id_stg_left  .OR.  myidx == id_stg_right )  THEN
+
+          IF ( myidx == id_stg_left  )  i = -1
+          IF ( myidx == id_stg_right )  i = nxr+1
+
           DO  j = nysg, nyng
              DO  k = nzb, nzt+1
 
                 IF  ( a11(k) .NE. 0._wp ) THEN
-                fu(k,j) = ( u(k,j,-1) / mc_factor -                    &
-                            mean_inflow_profiles(k,1) ) / a11(k)
+                   fu_yz(k,j) = ( u(k,j,i) / mc_factor - u_init(k) ) / a11(k)
                 ELSE
-                   fu(k,j) = 0._wp
+                   fu_yz(k,j) = 0._wp
                 ENDIF
 
                 IF  ( a22(k) .NE. 0._wp ) THEN
-                fv(k,j) = ( v(k,j,-1) / mc_factor - a21(k) * fu(k,j) - &
-                            mean_inflow_profiles(k,2) ) / a22(k)
+                   fv_yz(k,j) = ( v(k,j,i) / mc_factor - a21(k) * fu_yz(k,j) - &
+                               v_init(k) ) / a22(k)
                 ELSE
-                   fv(k,j) = 0._wp
+                   fv_yz(k,j) = 0._wp
                 ENDIF
 
                 IF  ( a33(k) .NE. 0._wp ) THEN
-                fw(k,j) = ( w(k,j,-1) / mc_factor - a31(k) * fu(k,j) - &
-                            a32(k) * fv(k,j) ) / a33(k)
+                   fw_yz(k,j) = ( w(k,j,i) / mc_factor - a31(k) * fu_yz(k,j) - &
+                               a32(k) * fv_yz(k,j) ) / a33(k)
                 ELSE
-                   fw = 0._wp
+                   fw_yz = 0._wp
                 ENDIF
 
@@ -759,5 +972,4 @@
     IMPLICIT NONE
 
-
     CALL wrd_write_string( 'mc_factor' )
     WRITE ( 14 )  mc_factor
@@ -785,9 +997,12 @@
 
     USE control_parameters,                                                    &
-        ONLY:  dt_3d, intermediate_timestep_count, simulated_time,             &
-               volume_flow_initial
+        ONLY:  dt_3d, forcing, intermediate_timestep_count,  nest_domain,      &
+               simulated_time, volume_flow_initial
+
+    USE grid_variables,                                                        &
+        ONLY:  dx, dy
 
     USE indices,                                                               &
-        ONLY:  nyn, nys
+        ONLY:  wall_flags_0
 
     USE statistics,                                                            &
@@ -797,4 +1012,5 @@
     IMPLICIT NONE
 
+    INTEGER(iwp) :: i           !< grid index in x-direction
     INTEGER(iwp) :: j           !< loop index in y-direction
     INTEGER(iwp) :: k           !< loop index in z-direction
@@ -802,6 +1018,9 @@
     REAL(wp) :: dt_stg          !< wheighted subtimestep
     REAL(wp) :: mc_factor_l     !< local mass flux correction factor
-
-    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,5,nbgp) :: inflow_dist !< disturbance for inflow profiles
+    REAL(wp) :: volume_flow     !< mass flux through lateral boundary
+    REAL(wp) :: volume_flow_l   !< local mass flux through lateral boundary
+
+    REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg,5) :: dist_xz !< imposed disturbances at north/south boundary
+    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,5) :: dist_yz !< imposed disturbances at left/right boundary
 
 
@@ -815,16 +1034,53 @@
 !-- Initial value of fu, fv, fw
     IF ( simulated_time == 0.0_wp .AND. .NOT. velocity_seed_initialized )  THEN
-       CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu )
-       CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv )
-       CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw )
+       CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_left )
+       CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_left )
+       CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, id_stg_left )
+
+       IF ( forcing  .OR.  nest_domain )  THEN
+!
+!--       Generate turbulence at right boundary
+          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_right )
+          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_right )
+          CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, id_stg_right )
+!
+!--       Generate turbulence at north boundary
+          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fu_xz, id_stg_north )
+          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fv_xz, id_stg_north )
+          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fw_xz, id_stg_north )
+!
+!--       Generate turbulence at south boundary
+          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fu_xz, id_stg_south )
+          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fv_xz, id_stg_south )
+          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fw_xz, id_stg_south )
+       ENDIF
        velocity_seed_initialized = .TRUE.
     ENDIF
 !
 !-- New set of fu, fv, fw
-    CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo )
-    CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo )
-    CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo )
-
-    IF ( myidx == id_inflow )  THEN
+    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 ( forcing  .OR.  nest_domain )  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
+!
+!-- Turbulence generation at left and or right boundary
+    IF ( myidx == id_stg_left  .OR.  myidx == id_stg_right )  THEN
 
        DO  j = nysg, nyng
@@ -833,22 +1089,22 @@
 !--          Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
              IF ( tu(k) == 0.0_wp )  THEN
-                fu(k,j) = fuo(k,j)
+                fu_yz(k,j) = fuo_yz(k,j)
              ELSE
-                fu(k,j) = fu(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) +     &
-                         fuo(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
+                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(k,j) = fvo(k,j)
+                fv_yz(k,j) = fvo_yz(k,j)
              ELSE
-                fv(k,j) = fv(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) +     &
-                         fvo(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
+                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(k,j) = fwo(k,j)
+                fw_yz(k,j) = fwo_yz(k,j)
              ELSE
-                fw(k,j) = fw(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) +     &
-                         fwo(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
+                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
 !
@@ -856,24 +1112,24 @@
 !--          Additional factors are added to improve the variance of v and w
              IF( k == 0 )  THEN
-                inflow_dist(k,j,1,1) = 0.0_wp
-                inflow_dist(k,j,2,1) = 0.0_wp
-                inflow_dist(k,j,3,1) = 0.0_wp
-                inflow_dist(k,j,4,1) = 0.0_wp
-                inflow_dist(k,j,5,1) = 0.0_wp
+                dist_yz(k,j,1) = 0.0_wp
+                dist_yz(k,j,2) = 0.0_wp
+                dist_yz(k,j,3) = 0.0_wp
+!                 dist_yz(k,j,4) = 0.0_wp
+!                 dist_yz(k,j,5) = 0.0_wp
              ELSE
-                inflow_dist(k,j,1,1) = a11(k) * fu(k,j)
+                dist_yz(k,j,1) = a11(k) * fu_yz(k,j)
                 !experimental test of 1.2
-                inflow_dist(k,j,2,1) = ( SQRT( a22(k) / MAXVAL(a22) )          &
+                dist_yz(k,j,2) = ( SQRT( a22(k) / MAXVAL(a22) )                   &
                                          * 1.2_wp )                            &
-                                       * (   a21(k) * fu(k,j)                  &
-                                           + a22(k) * fv(k,j) )
-                inflow_dist(k,j,3,1) = ( SQRT(a33(k) / MAXVAL(a33) )           &
+                                       * (   a21(k) * fu_yz(k,j)                  &
+                                           + a22(k) * fv_yz(k,j) )
+                dist_yz(k,j,3) = ( SQRT(a33(k) / MAXVAL(a33) )                    &
                                          * 1.3_wp )                            &
-                                       * (   a31(k) * fu(k,j)                  &
-                                           + a32(k) * fv(k,j)                  &
-                                           + a33(k) * fw(k,j) )
+                                       * (   a31(k) * fu_yz(k,j)                  &
+                                           + a32(k) * fv_yz(k,j)                  &
+                                           + a33(k) * fw_yz(k,j) )
                 ! Calculation for pt and e not yet implemented
-                inflow_dist(k,j,4,1) = 0.0_wp
-                inflow_dist(k,j,5,1) = 0.0_wp
+!                 dist_yz(k,j,4) = 0.0_wp
+!                 dist_yz(k,j,5) = 0.0_wp
              ENDIF
 
@@ -885,34 +1141,216 @@
 !--    This correction factor insures that the mass flux is preserved at the
 !--    inflow boundary
-       mc_factor_l = 0.0_wp
-       mc_factor   = 0.0_wp
-       DO  j = nys, nyn
+       IF ( .NOT. forcing  .AND.  .NOT. nest_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
+          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
+#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
+#endif
+
+          mc_factor = volume_flow / mc_factor
+
+!
+!--       Add disturbances
+          IF ( myidx == id_stg_left  )  THEN
+
+             DO  j = nysg, nyng
+                DO  k = nzb, nzt+1
+                   u(k,j,-nbgp+1:0) = ( u(k,j,-nbgp+1:0) + dist_yz(k,j,1) )    &
+                                        * mc_factor
+                   v(k,j,-nbgp:-1)  = ( v(k,j,-nbgp:-1)  + dist_yz(k,j,2)  )   &
+                                        * mc_factor
+                   w(k,j,-nbgp:-1)  = ( w(k,j,-nbgp:-1)  + dist_yz(k,j,3)  )   &
+                                        * mc_factor
+                ENDDO
+             ENDDO
+          ENDIF
+          IF ( myidx == id_stg_right  )  THEN
+
+             DO  j = nysg, nyng
+                DO  k = nzb, nzt+1
+                   u(k,j,nxr+1:nxr+nbgp) = ( u(k,j,nxr+1:nxr+nbgp) +           &
+                                             dist_yz(k,j,1) ) * mc_factor
+                   v(k,j,nxr+1:nxr+nbgp) = ( v(k,j,nxr+1:nxr+nbgp) +           &
+                                             dist_yz(k,j,2) ) * mc_factor
+                   w(k,j,nxr+1:nxr+nbgp) = ( w(k,j,nxr+1:nxr+nbgp) +           &
+                                             dist_yz(k,j,3) ) * mc_factor
+                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) = a11(k) * fu_xz(k,i)
+                !experimental test of 1.2
+                dist_xz(k,i,2) = ( SQRT( a22(k) / MAXVAL(a22) )                &
+                                         * 1.2_wp )                            &
+                                       * (   a21(k) * fu_xz(k,i)               &
+                                           + a22(k) * fv_xz(k,i) )
+                dist_xz(k,i,3) = ( 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) )
+             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
-             mc_factor_l = mc_factor_l + dzw(k)  *                             &
-                           ( mean_inflow_profiles(k,1) + inflow_dist(k,j,1,1) )
+             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( mc_factor_l, mc_factor,  &
-                           1, MPI_REAL, MPI_SUM, comm1dy, ierr )
+       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
-       mc_factor = mc_factor_l
+       volume_flow = volume_flow_l
+       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) +                 &
-                                   inflow_dist(k,j,1,1)      ) * mc_factor
-              v(k,j,-nbgp:-1)  = ( mean_inflow_profiles(k,2) +                 &
-                                   inflow_dist(k,j,2,1)      ) * mc_factor
-              w(k,j,-nbgp:-1)  =   inflow_dist(k,j,3,1)        * mc_factor
+       mc_factor = volume_flow / mc_factor
+
+!
+!--    Add disturbances
+       IF ( myidy == id_stg_south  )  THEN
+
+          DO  i = nxlg, nxrg
+             DO  k = nzb, nzt+1
+                u(k,-nbgp:-1,i) = ( u(k,-nbgp:-1,i) + dist_xz(k,i,1) )         &
+                                     * mc_factor
+                v(k,-nbgp:0,i)  = ( v(k,-nbgp:0,i)  + dist_xz(k,i,2)  )        &
+                                     * mc_factor
+                w(k,-nbgp:-1,i) = ( w(k,-nbgp:-1,i) + dist_xz(k,i,3)  )        &
+                                     * mc_factor
+             ENDDO
           ENDDO
-       ENDDO
-
+       ENDIF
+       IF ( myidy == id_stg_north  )  THEN
+
+          DO  i = nxlg, nxrg
+             DO  k = nzb, nzt+1
+                u(k,nyn+1:nyn+nbgp,i) = ( u(k,nyn+1:nyn+nbgp,i) +              &
+                                          dist_xz(k,i,1) ) * mc_factor
+                v(k,nyn+1:nyn+nbgp,i) = ( v(k,nyn+1:nyn+nbgp,i) +              &
+                                          dist_xz(k,i,2) ) * mc_factor
+                w(k,nyn+1:nyn+nbgp,i) = ( w(k,nyn+1:nyn+nbgp,i) +              &
+                                          dist_xz(k,i,3) ) * mc_factor
+             ENDDO
+          ENDDO
+       ENDIF
     ENDIF
 
@@ -920,5 +1358,4 @@
 
  END SUBROUTINE stg_main
-
 
 !------------------------------------------------------------------------------!
@@ -931,5 +1368,5 @@
 !> parts are collected to form the full array.
 !------------------------------------------------------------------------------!
- SUBROUTINE stg_generate_seed_yz( n_y, n_z, b_y, b_z, f_n )
+ SUBROUTINE stg_generate_seed_yz( n_y, n_z, b_y, b_z, f_n, id )
 
 
@@ -940,4 +1377,5 @@
     IMPLICIT NONE
 
+    INTEGER(iwp) :: id          !< core ids at respective boundaries
     INTEGER(iwp) :: j           !< loop index in y-direction
     INTEGER(iwp) :: jj          !< loop index in y-direction
@@ -957,5 +1395,5 @@
     REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_y     !< filter func in y-dir
     REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter func in z-dir
-    REAL(wp), DIMENSION(nzb_x:nzt_x+1,nysg:nyng) :: f_n_l   !<  local velocity seed
+    REAL(wp), DIMENSION(nzb_x_stg:nzt_x_stg+1,nysg:nyng) :: f_n_l   !<  local velocity seed
     REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng)     :: f_n     !<  velocity seed
     REAL(wp), DIMENSION(:,:), ALLOCATABLE        :: rand_it !<  random number
@@ -1018,5 +1456,5 @@
 
     DO  j = nysg, nyng
-       DO  k = nzb_x, nzt_x+1
+       DO  k = nzb_x_stg, nzt_x_stg+1
           DO  jj = -n_y2(k), n_y2(k)
              DO  kk = -n_z2(k), n_z2(k)
@@ -1032,13 +1470,13 @@
 !
 !-- Gather velocity seeds of full subdomain
-    send_count = nzt_x - nzb_x + 1
-    IF ( nzt_x == nzt )  send_count = send_count + 1
+    send_count = nzt_x_stg - nzb_x_stg + 1
+    IF ( nzt_x_stg == nzt )  send_count = send_count + 1
 
 #if defined( __parallel )
-    CALL MPI_GATHERV( f_n_l(nzb_x,nysg), send_count, stg_type_yz_small,        &
-                      f_n(nzb+1,nysg), recv_count, displs, stg_type_yz,        &
-                      id_inflow, comm1dx, ierr )
+    CALL MPI_GATHERV( f_n_l(nzb_x_stg,nysg), send_count, stg_type_yz_small,        &
+                      f_n(nzb+1,nysg), recv_count_yz, displs_yz, stg_type_yz,        &
+                      id, comm1dx, ierr )
 #else
-    f_n(nzb+1:nzt+1,nysg:nyng) = f_n_l(nzb_x:nzt_x+1,nysg:nyng)
+    f_n(nzb+1:nzt+1,nysg:nyng) = f_n_l(nzb_x_stg:nzt_x_stg+1,nysg:nyng)
 #endif
 
@@ -1046,3 +1484,132 @@
  END SUBROUTINE stg_generate_seed_yz
 
+
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Generate a set of random number rand_it wich is equal on each PE
+!> and calculate the velocity seed f_n.
+!> f_n is splitted in vertical direction by the number of PEs in y-direction and
+!> and each PE calculates a vertical subsection of f_n. At the the end, all
+!> parts are collected to form the full array.
+!------------------------------------------------------------------------------!
+ SUBROUTINE stg_generate_seed_xz( n_x, n_z, b_x, b_z, f_n, id )
+
+
+    USE indices,                                                               &
+        ONLY: nx
+
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) :: id          !< core ids at respective boundaries
+    INTEGER(iwp) :: i           !< loop index in x-direction
+    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
+    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x2   !< n_y*2
+    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z2   !< n_z*2
+
+    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(-merg:merg,nzb:nzt+1)    :: b_x     !< filter func in y-dir
+    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter func in z-dir
+    REAL(wp), DIMENSION(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg) :: f_n_l   !<  local velocity seed
+    REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg)     :: f_n     !<  velocity seed
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE        :: rand_it !<  random number
+
+
+!
+!-- Generate random numbers using a seed generated in stg_init.
+!-- The set of random numbers are modified to have an average of 0 and
+!-- unit variance.
+    ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp:nx+mergp) )
+
+    rand_av        = 0.0_wp
+    rand_sigma_inv = 0.0_wp
+    nxz_inv        = 1.0_wp / REAL( ( nzt+1 - nzb+1 ) * ( nx+1 ), KIND=wp )
+
+    DO  i = 0, nx
+       DO  k = nzb, nzt+1
+          CALL RANDOM_NUMBER( rand_it(k,i) )
+          rand_av = rand_av + rand_it(k,i)
+       ENDDO
+    ENDDO
+
+    rand_av = rand_av * nxz_inv
+
+    DO  i = 0, nx
+       DO  k = nzb, nzt+1
+          rand_it(k,i)   = rand_it(k,i) - rand_av
+          rand_sigma_inv = rand_sigma_inv + rand_it(k,i) ** 2
+       ENDDO
+    ENDDO
+
+    rand_sigma_inv = 1.0_wp / SQRT(rand_sigma_inv * nxz_inv)
+
+    DO  i = 0, nx
+       DO  k = nzb, nzt+1
+          rand_it(k,i) = rand_it(k,i) * rand_sigma_inv
+       ENDDO
+    ENDDO
+
+!
+!-- Periodic fill of random number in space
+    DO  i = 0, nx
+       DO  k = 1, mergp
+          rand_it(nzb-k,i)   = rand_it(nzt+2-k,i)    ! bottom margin
+          rand_it(nzt+1+k,i) = rand_it(nzb+k-1,i)    ! top margin
+       ENDDO
+    ENDDO
+    DO  i = 1, mergp
+       DO  k = nzb-mergp, nzt+1+mergp
+          rand_it(k,-i)   = rand_it(k,nx-i+1)        ! left margin
+          rand_it(k,nx+i) = rand_it(k,i-1)           ! right margin
+       ENDDO
+    ENDDO
+
+!
+!-- Generate velocity seed following Eq.6 of Xie and Castro (2008)
+    n_x2 = n_x * 2
+    n_z2 = n_z * 2
+    f_n_l  = 0.0_wp
+
+    DO  i = nxlg, nxrg
+       DO  k = nzb_y_stg, nzt_y_stg+1
+          DO  ii = -n_x2(k), n_x2(k)
+             DO  kk = -n_z2(k), n_z2(k)
+                f_n_l(k,i) = f_n_l(k,i)                                        &
+                           + b_x(ii,k) * b_z(kk,k) * rand_it(k+kk,i+ii)
+             ENDDO
+          ENDDO
+       ENDDO
+    ENDDO
+
+    DEALLOCATE( rand_it )
+
+!
+!-- Gather velocity seeds of full subdomain
+    send_count = nzt_y_stg - nzb_y_stg + 1
+    IF ( nzt_y_stg == nzt )  send_count = send_count + 1
+
+
+#if defined( __parallel )
+    CALL MPI_GATHERV( f_n_l(nzb_y_stg,nxlg), send_count, stg_type_xz_small,        &
+                      f_n(nzb+1,nxlg), recv_count_xz, displs_xz, stg_type_xz,  &
+                      id, comm1dy, ierr )
+#else
+    f_n(nzb+1:nzt+1,nxlg:nxrg) = f_n_l(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg)
+#endif
+
+
+ END SUBROUTINE stg_generate_seed_xz
+
  END MODULE synthetic_turbulence_generator_mod
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 2934)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 2938)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Nesting of dissipation rate in case of RANS mode and TKE-e closure is applied
+! 
+! 2936 2018-03-27 14:49:27Z suehring
 ! Little formatting adjustment.
 ! 
@@ -358,6 +361,6 @@
                microphysics_morrison, microphysics_seifert, mid, nest_domain,  &
                neutral, nr_timesteps_this_run, nudging,                        &
-               ocean, passive_scalar,                                          &
-               prho_reference, pt_reference, pt_slope_offset, random_heatflux, &
+               ocean, passive_scalar, prho_reference, pt_reference,            &
+               pt_slope_offset, random_heatflux, rans_mode,                    &
                rans_tke_e, run_coupled, simulated_time, simulated_time_chr,    &
                skip_time_do2d_xy, skip_time_do2d_xz, skip_time_do2d_yz,        &
@@ -782,4 +785,8 @@
                 IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e, nbgp )
 
+                IF ( .NOT. constant_diffusion  .AND.  rans_mode  .AND.         &
+                                                      rans_tke_e )             &
+                   CALL exchange_horiz( diss, nbgp )
+
                 IF ( air_chemistry )  THEN
                    DO  n = 1, nspec     
@@ -820,7 +827,5 @@
 !
 !--       Impose a turbulent inflow using synthetic generated turbulence
-          IF ( use_syn_turb_gen ) THEN 
-             CALL  stg_main
-          ENDIF
+          IF ( use_syn_turb_gen )  CALL  stg_main
 
 !
@@ -862,8 +867,8 @@
                                intermediate_timestep_count_max  )  THEN
              CALL forcing_bc
-          ENDIF
-!
-!--       Moreover, ensure mass conservation at each RK substep.
-          IF ( forcing )  CALL forcing_bc_mass_conservation
+!
+!--          Moreover, ensure mass conservation
+             CALL forcing_bc_mass_conservation
+          ENDIF
 
 !
Index: palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 2934)
+++ palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 2938)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Further todo's
+! 
+! 2936 2018-03-27 14:49:27Z suehring
 ! - defined l_grid only within this module
 ! - Moved l_wall definition from modules.f90
@@ -69,4 +72,6 @@
 !>       add OpenMP directives whereever possible
 !>       remove debug output variables (dummy1, dummy2, dummy3)
+!> @todo Move initialization of wall-mixing length from init_grid
+!> @todo Check for random disturbances
 !> @note <Enter notes on the module>
 !> @bug  TKE-e closure still crashes due to too small dt
@@ -284,5 +289,6 @@
 
     USE control_parameters,                                                    &
-        ONLY:  message_string, neutral, turbulent_inflow, turbulent_outflow
+        ONLY:  message_string, nest_domain, neutral, turbulent_inflow,         &
+               turbulent_outflow
 
     IMPLICIT NONE
@@ -302,6 +308,6 @@
              rans_tke_e = .TRUE.
 
-             IF ( INDEX( initializing_actions, 'set_1d-model_profiles' )       &
-                  == 0 )  THEN
+             IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) == 0  &
+                  .AND.  .NOT.  nest_domain )  THEN
                 message_string = 'Initializing without 1D model while ' //     &
                                  'using TKE-e closure&' //                     &
@@ -826,4 +832,7 @@
         ONLY:  use_sgs_for_particles, wang_kernel
 
+    USE pmc_interface,                                                         &
+        ONLY:  nested_run
+
     IMPLICIT NONE
 
@@ -847,7 +856,11 @@
     ALLOCATE( e_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 #endif
-
+!
+!-- Allocate arrays required for dissipation. 
+!-- Please note, if it is a nested run, arrays need to be allocated even if
+!-- they do not necessarily need to be transferred, which is attributed to 
+!-- the design of the model coupler which allocates memory for each variable. 
     IF ( rans_tke_e  .OR.  use_sgs_for_particles  .OR.  wang_kernel  .OR.      &
-         collision_turbulence )  THEN
+         collision_turbulence  .OR.  nested_run )  THEN
 #if defined( __nopointer )
        ALLOCATE( diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
@@ -858,5 +871,5 @@
 #else
        ALLOCATE( diss_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-       IF ( rans_tke_e )  THEN
+       IF ( rans_tke_e  .OR.  nested_run )  THEN
           ALLOCATE( diss_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
           ALLOCATE( diss_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
@@ -871,7 +884,7 @@
 
     IF ( rans_tke_e  .OR.  use_sgs_for_particles  .OR.     &
-         wang_kernel  .OR.  collision_turbulence )  THEN
+         wang_kernel  .OR.  collision_turbulence  .OR.  nested_run )  THEN
        diss => diss_1
-       IF ( rans_tke_e )  THEN
+       IF ( rans_tke_e  .OR.  nested_run )  THEN
        diss_p => diss_2; tdiss_m => diss_3
        ENDIF
