Index: /palm/trunk/SOURCE/Makefile
===================================================================
--- /palm/trunk/SOURCE/Makefile	(revision 3346)
+++ /palm/trunk/SOURCE/Makefile	(revision 3347)
@@ -25,4 +25,10 @@
 # -----------------
 # $Id$
+# Add module for offline nesting;
+# Add surface_mod to synthetic-turbulence generator;
+# Bugfix, missing dependency for turbulence generator in init_3d_model;
+# Some formatting ajdustments
+# 
+# 3343 2018-10-15 10:38:52Z suehring
 # (from branch resler)
 # Add biometeorology
@@ -566,4 +572,5 @@
 	netcdf_data_input_mod.f90 \
 	netcdf_interface_mod.f90 \
+	nesting_offl_mod.f90 \
 	ocean_mod.f90 \
 	outflow_turbulence.f90 \
@@ -1005,4 +1012,5 @@
 	modules.o \
 	netcdf_interface_mod.o \
+	nesting_offl_mod.o \
 	ocean_mod.o \
 	plant_canopy_model_mod.o \
@@ -1036,4 +1044,5 @@
 	netcdf_data_input_mod.o \
 	netcdf_interface_mod.o \
+	nesting_offl_mod.o \
 	ocean_mod.o \
 	plant_canopy_model_mod.o \
@@ -1045,4 +1054,5 @@
 	surface_layer_fluxes_mod.o \
 	surface_mod.o \
+	synthetic_turbulence_generator_mod.o \
 	turbulence_closure_mod.o \
 	urban_surface_mod.o \
@@ -1112,5 +1122,4 @@
 	mod_kinds.o \
 	modules.o \
-	netcdf_data_input_mod.o \
 	surface_mod.o
 local_stop.o: \
@@ -1288,4 +1297,10 @@
 	urban_surface_mod.o \
 	uv_exposure_model_mod.o
+nesting_offl_mod.o: \
+        cpulog_mod.o \
+        mod_kinds.o \
+	modules.o \
+	netcdf_data_input_mod.o \
+	surface_mod.o
 ocean_mod.o: \
 	advec_s_pw.o \
@@ -1335,4 +1350,5 @@
 	modules.o \
 	multi_agent_system_mod.o \
+	nesting_offl_mod.o \
 	netcdf_interface_mod.o \
 	ocean_mod.o \
@@ -1577,5 +1593,7 @@
 	mod_kinds.o \
 	modules.o \
-	pmc_interface_mod.o
+	nesting_offl_mod.o \
+	pmc_interface_mod.o \
+	surface_mod.o
 temperton_fft_mod.o: \
 	mod_kinds.o \
@@ -1600,4 +1618,5 @@
 	modules.o \
 	multi_agent_system_mod.o \
+	nesting_offl_mod.o \
 	ocean_mod.o \
 	pmc_interface_mod.o \
Index: /palm/trunk/SOURCE/boundary_conds.f90
===================================================================
--- /palm/trunk/SOURCE/boundary_conds.f90	(revision 3346)
+++ /palm/trunk/SOURCE/boundary_conds.f90	(revision 3347)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Bugfix in setting boundary conditions in offline nesting
+! 
+! 3341 2018-10-15 10:31:27Z suehring
 ! changes concerning modularization of ocean option
 ! 
@@ -293,5 +296,5 @@
 !
 !-- Vertical nesting: Vertical velocity not zero at the top of the fine grid
-    IF (  .NOT.  child_domain  .AND.                                           &
+    IF (  .NOT.  child_domain  .AND.  .NOT.  nesting_offline  .AND.            &
                  TRIM(coupling_mode) /= 'vnested_fine' )  THEN
        w_p(nzt:nzt+1,:,:) = 0.0_wp  !< nzt is not a prognostic level (but cf. pres)
@@ -729,6 +732,4 @@
 !-- in case of nest boundaries. This must not be done in case of vertical nesting 
 !-- mode as in that case the lateral boundaries are actually cyclic. 
-!-- @todo: Is this really needed? Boundary values will be overwritten in 
-!--        coupler or by Inifor data.
     IF ( nesting_mode /= 'vertical'  .OR.  nesting_offline )  THEN
        IF ( bc_dirichlet_s )  THEN
Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 3346)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 3347)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Remove offline nesting in if clause for pressure top boundary condition
+! 
+! 3343 2018-10-15 10:38:52Z suehring
 ! (from branch resler)
 ! Add biometeorology,
@@ -1954,6 +1957,5 @@
        ibc_p_t = 0
 !-- TO_DO: later set bc_p_t to neumann before, in case of nested domain
-    ELSEIF ( bc_p_t == 'neumann' .OR. bc_p_t == 'nested'  .OR.                 &
-             bc_p_t == 'nesting_offline'  )  THEN
+    ELSEIF ( bc_p_t == 'neumann' .OR. bc_p_t == 'nested' )  THEN
        ibc_p_t = 1
     ELSE
Index: /palm/trunk/SOURCE/header.f90
===================================================================
--- /palm/trunk/SOURCE/header.f90	(revision 3346)
+++ /palm/trunk/SOURCE/header.f90	(revision 3347)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Header output concerning offline nesting
+! 
+! 3343 2018-10-15 10:38:52Z suehring
 ! call of ocean_header
 ! 
@@ -438,4 +441,7 @@
     USE netcdf_interface,                                                      &
         ONLY:  netcdf_data_format, netcdf_data_format_string, netcdf_deflate
+        
+    USE nesting_offl_mod,                                                      &
+        ONLY:  nesting_offl_header
 
     USE ocean_mod,                                                             &
@@ -903,4 +909,8 @@
 !-- Large scale forcing and nudging
     IF ( large_scale_forcing )  CALL lsf_nudging_header( io )
+    
+!
+!-- Offline nesting 
+    IF ( nesting_offline )  CALL nesting_offl_header( io )
 
 !
Index: /palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_3d_model.f90	(revision 3346)
+++ /palm/trunk/SOURCE/init_3d_model.f90	(revision 3347)
@@ -25,4 +25,12 @@
 ! -----------------
 ! $Id$
+! - Separate offline nesting from large_scale_nudging_mod
+! - Improve the synthetic turbulence generator
+! 
+! 3298 2018-10-02 12:21:11Z kanani
+! Minor formatting (kanani)
+! Added CALL to the chem_emissions module (Russo)
+! 
+! 3294 2018-10-01 02:37:10Z raasch
 ! allocate and set stokes drift velocity profiles
 ! 
@@ -563,4 +571,7 @@
         ONLY:  chem_emis, chem_emis_att, init_3d,                              &
                netcdf_data_input_init_3d, netcdf_data_input_interpolate
+        
+    USE nesting_offl_mod,                                                      &
+        ONLY:  nesting_offl_init
 
     USE ocean_mod,                                                             &
@@ -601,6 +612,7 @@
 
     USE synthetic_turbulence_generator_mod,                                    &
-        ONLY:  stg_init, use_syn_turb_gen
-
+        ONLY:  parametrize_inflow_turbulence, stg_adjust, stg_init,            &
+               use_syn_turb_gen
+               
     USE surface_layer_fluxes_mod,                                              &
         ONLY:  init_surface_layer_fluxes
@@ -2152,7 +2164,7 @@
     IF ( nudging )  CALL nudge_init
 !
-!-- Initialize 1D/3D offline-nesting with COSMO model and read data from 
-!-- external file.
-    IF ( large_scale_forcing  .OR.  nesting_offline )  CALL lsf_init
+!-- Initialize 1D large-scale forcing and nudging and read data from external
+!-- ASCII file
+    IF ( large_scale_forcing )  CALL lsf_init
 !
 !-- Initialize surface forcing corresponding to large-scale forcing. Therein, 
@@ -2163,4 +2175,8 @@
        ENDIF
     ENDIF
+!
+!-- Initializae 3D offline nesting in COSMO model and read data from 
+!-- external NetCDF file.
+    IF ( nesting_offline )  CALL nesting_offl_init
 !
 !-- Initialize quantities for special advections schemes
@@ -2279,5 +2295,12 @@
        CALL location_message( 'finished', .TRUE. )
     ENDIF
-
+!
+!-- In case the synthetic turbulence generator does not have any information
+!-- about the inflow turbulence, these information will be parametrized 
+!-- depending on the initial atmospheric conditions and surface properties.
+!-- Please note, within pre-determined time intervals these turbulence 
+!-- information can be updated if desired. 
+    IF ( use_syn_turb_gen  .AND.  parametrize_inflow_turbulence )              &
+       CALL stg_adjust
 !
 !-- If required, set chemical emissions
Index: /palm/trunk/SOURCE/init_pegrid.f90
===================================================================
--- /palm/trunk/SOURCE/init_pegrid.f90	(revision 3346)
+++ /palm/trunk/SOURCE/init_pegrid.f90	(revision 3347)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Bugfix in setting number of ghost layers in neutral case 
+! 
+! 3341 2018-10-15 10:31:27Z suehring
 ! unused variables removed
 ! 
@@ -723,5 +726,5 @@
 ! 
 !-- Determine the number of ghost point layers
-    IF ( ( scalar_advec == 'ws-scheme' .AND. .NOT. neutral ) .OR.             &
+    IF ( scalar_advec   == 'ws-scheme'  .OR.                                   &
          momentum_advec == 'ws-scheme' )  THEN
        nbgp = 3
Index: /palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 3346)
+++ /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 3347)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Assign real value instead of integer
+! 
+! 3341 2018-10-15 10:31:27Z suehring
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -2781,7 +2784,7 @@
                          vegetation_type_f%var(j,i)                 = 1 ! bare soil
                          soil_type_f%var_2d(j,i)                    = 1 
-                         surface_fraction_f%frac(ind_veg_wall,j,i)  = 1
-                         surface_fraction_f%frac(ind_pav_green,j,i) = 0
-                         surface_fraction_f%frac(ind_wat_win,j,i)   = 0
+                         surface_fraction_f%frac(ind_veg_wall,j,i)  = 1.0_wp
+                         surface_fraction_f%frac(ind_pav_green,j,i) = 0.0_wp
+                         surface_fraction_f%frac(ind_wat_win,j,i)   = 0.0_wp
                       ENDIF
                       
Index: /palm/trunk/SOURCE/large_scale_forcing_nudging_mod.f90
===================================================================
--- /palm/trunk/SOURCE/large_scale_forcing_nudging_mod.f90	(revision 3346)
+++ /palm/trunk/SOURCE/large_scale_forcing_nudging_mod.f90	(revision 3347)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Offline nesting parts are moved to an independent module nesting_offl_mod
+! 
+! 3341 2018-10-15 10:31:27Z suehring
 ! ocean renamed ocean_mode
 ! 
@@ -87,6 +90,5 @@
 
     USE control_parameters,                                                    &
-        ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, &
-               bc_lr, bc_ns, bc_pt_b, bc_q_b, constant_diffusion,              &
+        ONLY:  bc_lr, bc_ns, bc_pt_b, bc_q_b, constant_diffusion,              &
                constant_heatflux, constant_waterflux,                          &
                data_output_pr, dt_3d, end_time,                                &
@@ -94,8 +96,11 @@
                ibc_pt_b, ibc_q_b,                                              &
                large_scale_forcing, large_scale_subsidence, lsf_surf, lsf_vert,&
-               lsf_exception, message_string, nesting_offline, neutral,        &
+               lsf_exception, message_string, neutral,                         &
                nudging, passive_scalar, pt_surface, ocean_mode, q_surface,     &
                surface_heatflux, surface_pressure, surface_waterflux,          &
                topography, use_subsidence_tendencies
+               
+    USE cpulog,                                                                &
+        ONLY:  cpu_log, log_point
 
     USE grid_variables
@@ -107,7 +112,4 @@
     USE kinds
 
-    USE netcdf_data_input_mod,                                                 &
-        ONLY:  nest_offl
-
     USE pegrid
 
@@ -120,6 +122,4 @@
     INTEGER(iwp) ::  nlsf = 1000                       !< maximum number of profiles in LSF_DATA (large scale forcing)
     INTEGER(iwp) ::  ntnudge = 1000                    !< maximum number of profiles in NUDGING_DATA (nudging)
-
-    REAL(wp) ::  d_area_t
 
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ptnudge     !< vertical profile of pot. temperature interpolated to vertical grid (nudging)
@@ -156,5 +156,4 @@
            lsf_nudging_check_parameters, nudge_init,                           &
            lsf_nudging_check_data_output_pr, lsf_nudging_header,               &
-           lsf_nesting_offline, lsf_nesting_offline_mass_conservation,         & 
            nudge, nudge_ref
            
@@ -177,482 +176,4 @@
  CONTAINS
 
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> In this subroutine a constant mass within the model domain is guaranteed. 
-!> Larger-scale models may be based on a compressible equation system, which is
-!> not consistent with PALMs incompressible equation system. In order to avoid
-!> a decrease or increase of mass during the simulation, non-divergent flow
-!> through the lateral and top boundaries is compensated by the vertical wind 
-!> component at the top boundary. 
-!------------------------------------------------------------------------------!
-    SUBROUTINE lsf_nesting_offline_mass_conservation
-
-       USE control_parameters,                                                 &
-           ONLY:  volume_flow
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i !< grid index in x-direction
-       INTEGER(iwp) ::  j !< grid index in y-direction
-       INTEGER(iwp) ::  k !< grid index in z-direction
-
-       REAL(wp) ::  w_correct                       !< vertical velocity increment required to compensate non-divergent flow through the boundaries
-       REAL(wp), DIMENSION(1:3) ::  volume_flow_l   !< local volume flow
-
-       volume_flow   = 0.0_wp
-       volume_flow_l = 0.0_wp
-
-       d_area_t = 1.0_wp / ( ( nx + 1 ) * dx * ( ny + 1 ) * dy )
-
-       IF ( bc_dirichlet_l )  THEN
-          i = nxl
-          DO  j = nys, nyn
-             DO  k = nzb+1, nzt
-                volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) * dy   &
-                                     * MERGE( 1.0_wp, 0.0_wp,                  &
-                                              BTEST( wall_flags_0(k,j,i), 1 ) )
-             ENDDO
-          ENDDO
-       ENDIF
-       IF ( bc_dirichlet_r )  THEN
-          i = nxr+1
-          DO  j = nys, nyn
-             DO  k = nzb+1, nzt
-                volume_flow_l(1) = volume_flow_l(1) - u(k,j,i) * dzw(k) * dy   &
-                                     * MERGE( 1.0_wp, 0.0_wp,                  &
-                                              BTEST( wall_flags_0(k,j,i), 1 ) )
-             ENDDO
-          ENDDO
-       ENDIF
-       IF ( bc_dirichlet_s )  THEN
-          j = nys
-          DO  i = nxl, nxr
-             DO  k = nzb+1, nzt
-                volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) * dx   &
-                                     * MERGE( 1.0_wp, 0.0_wp,                  &
-                                              BTEST( wall_flags_0(k,j,i), 2 ) )
-             ENDDO
-          ENDDO
-       ENDIF
-       IF ( bc_dirichlet_n )  THEN
-          j = nyn+1
-          DO  i = nxl, nxr
-             DO  k = nzb+1, nzt
-                volume_flow_l(2) = volume_flow_l(2) - v(k,j,i) * dzw(k) * dx   &
-                                     * MERGE( 1.0_wp, 0.0_wp,                  &
-                                              BTEST( wall_flags_0(k,j,i), 2 ) )
-             ENDDO
-          ENDDO
-       ENDIF
-!
-!--    Top boundary
-       k = nzt
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dx * dy
-          ENDDO
-       ENDDO
-
-#if defined( __parallel )
-       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-       CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, 3, MPI_REAL, MPI_SUM,      &
-                           comm2d, ierr )
-#else
-       volume_flow = volume_flow_l
-#endif
-
-       w_correct = SUM( volume_flow ) * d_area_t
-
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             DO  k = nzt, nzt + 1
-                w(k,j,i) = w(k,j,i) + w_correct
-             ENDDO
-          ENDDO
-       ENDDO
-
-    END SUBROUTINE lsf_nesting_offline_mass_conservation
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Set the lateral and top boundary conditions in case the PALM domain is 
-!> nested offline in a mesoscale model. 
-!------------------------------------------------------------------------------!
-    SUBROUTINE lsf_nesting_offline
-
-       USE control_parameters,                                                 &
-           ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,              &
-                  bc_dirichlet_s, humidity, neutral, passive_scalar, rans_mode,&
-                  rans_tke_e, time_since_reference_point                      
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i !< running index x-direction
-       INTEGER(iwp) ::  j !< running index y-direction
-       INTEGER(iwp) ::  k !< running index z-direction
-
-       REAL(wp) ::  ddt_lsf !< inverse value of time resolution of forcing data
-       REAL(wp) ::  t_ref   !< time past since last reference step
-     
-!
-!--    Calculate time interval of forcing data       
-       ddt_lsf = 1.0_wp / ( nest_offl%time(nest_offl%tind_p) -                 &
-                            nest_offl%time(nest_offl%tind) )
-!
-!--    Calculate reziproke time past since last reference step. Please note, 
-!--    the time coordinate is still not updated, so that the actual time need
-!--    to be incremented by dt_3d. Moreover, note that the simulation time
-!--    passed since simulation start is time_since_reference_point, not 
-!--    simulated_time!
-       t_ref = time_since_reference_point + dt_3d -                            &
-                                            nest_offl%time(nest_offl%tind)
-                                            
-       IF ( bc_dirichlet_l )  THEN
-
-          DO  j = nys, nyn
-             DO  k = nzb+1, nzt
-                u(k,j,nxlg:nxl) = nest_offl%u_left(0,k,j) + ddt_lsf * t_ref  * &
-                       ( nest_offl%u_left(1,k,j) - nest_offl%u_left(0,k,j) ) * &
-                         MERGE( 1.0_wp, 0.0_wp,                                &
-                                BTEST( wall_flags_0(k,j,nxlg:nxl), 1 ) )
-             ENDDO
-          ENDDO
-
-          DO  j = nys, nyn
-             DO  k = nzb+1, nzt-1
-                w(k,j,nxlg:nxl-1) = nest_offl%w_left(0,k,j) + ddt_lsf * t_ref *&
-                       ( nest_offl%w_left(1,k,j) - nest_offl%w_left(0,k,j) )  *&
-                         MERGE( 1.0_wp, 0.0_wp,                                &
-                                BTEST( wall_flags_0(k,j,nxlg:nxl-1), 3 ) )
-             ENDDO
-          ENDDO
-
-          DO  j = nysv, nyn
-             DO  k = nzb+1, nzt
-                v(k,j,nxlg:nxl-1) = nest_offl%v_left(0,k,j) + ddt_lsf * t_ref *&
-                       ( nest_offl%v_left(1,k,j) - nest_offl%v_left(0,k,j) )  *&
-                         MERGE( 1.0_wp, 0.0_wp,                                &
-                                BTEST( wall_flags_0(k,j,nxlg:nxl-1), 2 ) )
-             ENDDO
-          ENDDO
-
-          IF ( .NOT. neutral )  THEN
-             DO  j = nys, nyn
-                DO  k = nzb+1, nzt
-                   pt(k,j,nxlg:nxl-1) = nest_offl%pt_left(0,k,j) + ddt_lsf *   &
-                                                                   t_ref   *   &
-                       ( nest_offl%pt_left(1,k,j) - nest_offl%pt_left(0,k,j) )
- 
-                ENDDO
-             ENDDO
-          ENDIF
-
-          IF ( humidity )  THEN
-             DO  j = nys, nyn
-                DO  k = nzb+1, nzt
-                   q(k,j,nxlg:nxl-1) = nest_offl%q_left(0,k,j) + ddt_lsf   *   &
-                                                                 t_ref     *   &
-                       ( nest_offl%q_left(1,k,j) - nest_offl%q_left(0,k,j) )
- 
-                ENDDO
-             ENDDO
-          ENDIF
-
-       ENDIF
-
-       IF ( bc_dirichlet_r )  THEN
-
-          DO  j = nys, nyn
-             DO  k = nzb+1, nzt
-                u(k,j,nxr+1:nxrg) = nest_offl%u_right(0,k,j) + ddt_lsf * t_ref *&
-                      ( nest_offl%u_right(1,k,j) - nest_offl%u_right(0,k,j) )  *&
-                         MERGE( 1.0_wp, 0.0_wp,                                 &
-                                BTEST( wall_flags_0(k,j,nxr+1:nxrg), 1 ) )
-
-             ENDDO
-          ENDDO
-          DO  j = nys, nyn
-             DO  k = nzb+1, nzt-1
-                w(k,j,nxr+1:nxrg) = nest_offl%w_right(0,k,j) + ddt_lsf * t_ref *&
-                      ( nest_offl%w_right(1,k,j) - nest_offl%w_right(0,k,j) )  *&
-                         MERGE( 1.0_wp, 0.0_wp,                                 &
-                                BTEST( wall_flags_0(k,j,nxr+1:nxrg), 3 ) )
-             ENDDO
-          ENDDO
-
-          DO  j = nysv, nyn
-             DO  k = nzb+1, nzt
-                v(k,j,nxr+1:nxrg) = nest_offl%v_right(0,k,j) + ddt_lsf * t_ref *&
-                      ( nest_offl%v_right(1,k,j) - nest_offl%v_right(0,k,j) )  *&
-                         MERGE( 1.0_wp, 0.0_wp,                                 &
-                                BTEST( wall_flags_0(k,j,nxr+1:nxrg), 2 ) )
-             ENDDO
-          ENDDO
-
-          IF ( .NOT. neutral )  THEN
-             DO  j = nys, nyn
-                DO  k = nzb+1, nzt
-                   pt(k,j,nxr+1:nxrg) = nest_offl%pt_right(0,k,j) + ddt_lsf *  &
-                                                                    t_ref   *  &
-                     ( nest_offl%pt_right(1,k,j) - nest_offl%pt_right(0,k,j) )
- 
-                ENDDO
-             ENDDO
-          ENDIF
-
-          IF ( humidity )  THEN
-             DO  j = nys, nyn
-                DO  k = nzb+1, nzt
-                   q(k,j,nxr+1:nxrg) = nest_offl%q_right(0,k,j) + ddt_lsf   *  &
-                                                                    t_ref   *  &
-                       ( nest_offl%q_right(1,k,j) - nest_offl%q_right(0,k,j) )
- 
-                ENDDO
-             ENDDO
-          ENDIF
-
-       ENDIF
-
-       IF ( bc_dirichlet_s )  THEN
-
-          DO  i = nxl, nxr
-             DO  k = nzb+1, nzt
-                v(k,nysg:nys,i)   = nest_offl%v_south(0,k,i) + ddt_lsf * t_ref *&
-                      ( nest_offl%v_south(1,k,i) - nest_offl%v_south(0,k,i) )  *&
-                         MERGE( 1.0_wp, 0.0_wp,                                 &
-                                BTEST( wall_flags_0(k,nysg:nys,i), 2 ) )
-             ENDDO
-          ENDDO
-
-          DO  i = nxl, nxr
-             DO  k = nzb+1, nzt-1
-                w(k,nysg:nys-1,i) = nest_offl%w_south(0,k,i) + ddt_lsf * t_ref  *&
-                        ( nest_offl%w_south(1,k,i) - nest_offl%w_south(0,k,i) ) *&
-                           MERGE( 1.0_wp, 0.0_wp,                                &
-                                  BTEST( wall_flags_0(k,nysg:nys-1,i), 3 ) )
-             ENDDO
-          ENDDO
-
-          DO  i = nxlu, nxr
-             DO  k = nzb+1, nzt
-                u(k,nysg:nys-1,i) = nest_offl%u_south(0,k,i) + ddt_lsf * t_ref  *&
-                        ( nest_offl%u_south(1,k,i) - nest_offl%u_south(0,k,i) ) *&
-                           MERGE( 1.0_wp, 0.0_wp,                                &
-                                  BTEST( wall_flags_0(k,nysg:nys-1,i), 1 ) )
-             ENDDO
-          ENDDO
-
-          IF ( .NOT. neutral )  THEN
-             DO  i = nxl, nxr
-                DO  k = nzb+1, nzt
-                   pt(k,nysg:nys-1,i) = nest_offl%pt_south(0,k,i) + ddt_lsf *  &
-                                                                    t_ref   *  &
-                     ( nest_offl%pt_south(1,k,i) - nest_offl%pt_south(0,k,i) )
- 
-                ENDDO
-             ENDDO
-          ENDIF
-
-          IF ( humidity )  THEN
-             DO  i = nxl, nxr
-                DO  k = nzb+1, nzt
-                   q(k,nysg:nys-1,i) = nest_offl%q_south(0,k,i) + ddt_lsf   *  &
-                                                                    t_ref   *  &
-                       ( nest_offl%q_south(1,k,i) - nest_offl%q_south(0,k,i) )
- 
-                ENDDO
-             ENDDO
-          ENDIF
-
-       ENDIF
-
-       IF ( bc_dirichlet_n )  THEN
-
-          DO  i = nxl, nxr
-             DO  k = nzb+1, nzt
-                v(k,nyn+1:nyng,i)   = nest_offl%v_north(0,k,i) + ddt_lsf * t_ref *&
-                        ( nest_offl%v_north(1,k,i) - nest_offl%v_north(0,k,i) )  *&
-                           MERGE( 1.0_wp, 0.0_wp,                                 &
-                                  BTEST( wall_flags_0(k,nyn+1:nyng,i), 2 ) )
-             ENDDO
-          ENDDO
-          DO  i = nxl, nxr
-             DO  k = nzb+1, nzt-1
-                w(k,nyn+1:nyng,i) = nest_offl%w_north(0,k,i) + ddt_lsf * t_ref  *&
-                        ( nest_offl%w_north(1,k,i) - nest_offl%w_north(0,k,i) ) *&
-                           MERGE( 1.0_wp, 0.0_wp,                                &
-                                  BTEST( wall_flags_0(k,nyn+1:nyng,i), 3 ) )
-             ENDDO
-          ENDDO
-
-          DO  i = nxlu, nxr
-             DO  k = nzb+1, nzt
-                u(k,nyn+1:nyng,i) = nest_offl%u_north(0,k,i) + ddt_lsf * t_ref  *&
-                        ( nest_offl%u_north(1,k,i) - nest_offl%u_north(0,k,i) ) *&
-                           MERGE( 1.0_wp, 0.0_wp,                                &
-                                  BTEST( wall_flags_0(k,nyn+1:nyng,i), 1 ) )
-
-             ENDDO
-          ENDDO
-
-          IF ( .NOT. neutral )  THEN
-             DO  i = nxl, nxr
-                DO  k = nzb+1, nzt
-                   pt(k,nyn+1:nyng,i) = nest_offl%pt_north(0,k,i) + ddt_lsf *  &
-                                                                    t_ref   *  &
-                     ( nest_offl%pt_north(1,k,i) - nest_offl%pt_north(0,k,i) )
- 
-                ENDDO
-             ENDDO
-          ENDIF
-
-          IF ( humidity )  THEN
-             DO  i = nxl, nxr
-                DO  k = nzb+1, nzt
-                   q(k,nyn+1:nyng,i) = nest_offl%q_north(0,k,i) + ddt_lsf   *  &
-                                                                    t_ref   *  &
-                       ( nest_offl%q_north(1,k,i) - nest_offl%q_north(0,k,i) )
- 
-                ENDDO
-             ENDDO
-          ENDIF
-
-       ENDIF
-!
-!--    Top boundary. 
-       DO  i = nxlu, nxr
-          DO  j = nys, nyn
-             u(nzt+1,j,i) = nest_offl%u_top(0,j,i) + ddt_lsf * t_ref *         &
-                        ( nest_offl%u_top(1,j,i) - nest_offl%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) = nest_offl%v_top(0,j,i) + ddt_lsf * t_ref *         &
-                        ( nest_offl%v_top(1,j,i) - nest_offl%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) = nest_offl%w_top(0,j,i) + ddt_lsf * t_ref *     &
-                        ( nest_offl%w_top(1,j,i) - nest_offl%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) = nest_offl%pt_top(0,j,i) + ddt_lsf * t_ref *    &
-                        ( nest_offl%pt_top(1,j,i) - nest_offl%pt_top(0,j,i) )
-             ENDDO
-          ENDDO
-       ENDIF
-
-       IF ( humidity )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                q(nzt+1,j,i) = nest_offl%q_top(0,j,i) + ddt_lsf * t_ref *      &
-                        ( nest_offl%q_top(1,j,i) - nest_offl%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 ( bc_dirichlet_l  .AND.  bc_dirichlet_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 ( bc_dirichlet_l  .AND.  bc_dirichlet_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 ( bc_dirichlet_r  .AND.  bc_dirichlet_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 ( bc_dirichlet_r  .AND.  bc_dirichlet_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
-!
-!--    Moreover, set Neumann boundary condition for subgrid-scale TKE, 
-!--    passive scalar, dissipation, and chemical species if required
-       IF ( rans_mode  .AND.  rans_tke_e )  THEN
-          IF (  bc_dirichlet_l )  diss(:,:,nxl-1) = diss(:,:,nxl)
-          IF (  bc_dirichlet_r )  diss(:,:,nxr+1) = diss(:,:,nxr)
-          IF (  bc_dirichlet_s )  diss(:,nys-1,:) = diss(:,nys,:)
-          IF (  bc_dirichlet_n )  diss(:,nyn+1,:) = diss(:,nyn,:)
-       ENDIF
-       IF ( .NOT. constant_diffusion )  THEN
-          IF (  bc_dirichlet_l )  e(:,:,nxl-1) = e(:,:,nxl)
-          IF (  bc_dirichlet_r )  e(:,:,nxr+1) = e(:,:,nxr)
-          IF (  bc_dirichlet_s )  e(:,nys-1,:) = e(:,nys,:)
-          IF (  bc_dirichlet_n )  e(:,nyn+1,:) = e(:,nyn,:)
-          e(nzt+1,:,:) = e(nzt,:,:)
-       ENDIF
-       IF ( passive_scalar )  THEN
-          IF (  bc_dirichlet_l )  s(:,:,nxl-1) = s(:,:,nxl)
-          IF (  bc_dirichlet_r )  s(:,:,nxr+1) = s(:,:,nxr)
-          IF (  bc_dirichlet_s )  s(:,nys-1,:) = s(:,nys,:)
-          IF (  bc_dirichlet_n )  s(:,nyn+1,:) = s(:,nyn,:)
-       ENDIF
-
-
-
-       CALL exchange_horiz( u, nbgp )
-       CALL exchange_horiz( v, nbgp )
-       CALL exchange_horiz( w, nbgp )
-       IF ( .NOT. neutral )  CALL exchange_horiz( pt, nbgp )
-       IF ( humidity      )  CALL exchange_horiz( q,  nbgp )
-
-!
-!--    Set surface pressure. Please note, time-dependent surface
-!--    pressure would require changes in anelastic approximation and 
-!--    treatment of fluxes. 
-!--    For the moment, comment this out! 
-!      surface_pressure = nest_offl%surface_pressure(nest_offl%tind) +                 &
-!                                                      ddt_lsf * t_ref *       &
-!                                    ( nest_offl%surface_pressure(nest_offl%tind_p)    &
-!                                    - nest_offl%surface_pressure(nest_offl%tind) )
-
-    END SUBROUTINE lsf_nesting_offline
 
 !------------------------------------------------------------------------------!
@@ -935,7 +456,4 @@
            ONLY:  bc_lr_cyc, bc_ns_cyc
 
-       USE netcdf_data_input_mod,                                              &
-           ONLY:  netcdf_data_input_lsf  
-
        IMPLICIT NONE
 
@@ -967,107 +485,105 @@
        REAL(wp) ::  r_dummy           !<
 
-       IF ( nesting_offline )  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( nest_offl%ug(0:1,nzb:nzt+1) )
-             ALLOCATE( nest_offl%vg(0:1,nzb:nzt+1) )
+       ALLOCATE( p_surf(0:nlsf), pt_surf(0:nlsf), q_surf(0:nlsf),              &
+                 qsws_surf(0:nlsf), shf_surf(0:nlsf),                          &
+                 td_lsa_lpt(nzb:nzt+1,0:nlsf), td_lsa_q(nzb:nzt+1,0:nlsf),     &
+                 td_sub_lpt(nzb:nzt+1,0:nlsf), td_sub_q(nzb:nzt+1,0:nlsf),     &
+                 time_vert(0:nlsf), time_surf(0:nlsf),                         &
+                 ug_vert(nzb:nzt+1,0:nlsf), vg_vert(nzb:nzt+1,0:nlsf),         &
+                 wsubs_vert(nzb:nzt+1,0:nlsf) )
+
+       p_surf = 0.0_wp; pt_surf = 0.0_wp; q_surf = 0.0_wp; qsws_surf = 0.0_wp
+       shf_surf = 0.0_wp; time_vert = 0.0_wp; td_lsa_lpt = 0.0_wp
+       td_lsa_q = 0.0_wp; td_sub_lpt = 0.0_wp; td_sub_q = 0.0_wp
+       time_surf = 0.0_wp; ug_vert = 0.0_wp; vg_vert = 0.0_wp
+       wsubs_vert = 0.0_wp
+
+!
+!--    Array for storing large scale forcing and nudging tendencies at each 
+!--    timestep for data output
+       ALLOCATE( sums_ls_l(nzb:nzt+1,0:7) )
+       sums_ls_l = 0.0_wp
+
+       ngp_sums_ls = (nz+2)*6
+
+       OPEN ( finput, FILE='LSF_DATA', STATUS='OLD', &
+              FORM='FORMATTED', IOSTAT=ierrn )
+
+       IF ( ierrn /= 0 )  THEN
+          message_string = 'file LSF_DATA does not exist'
+          CALL message( 'ls_forcing', 'PA0368', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       ierrn = 0
+!
+!--    First three lines of LSF_DATA contain header
+       READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
+       READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
+       READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
+
+       IF ( ierrn /= 0 )  THEN
+          message_string = 'errors in file LSF_DATA'
+          CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
+       ENDIF
+
+!
+!--    Surface values are read in
+       nt     = 0
+       ierrn = 0
+
+       DO WHILE ( time_surf(nt) < end_time )
+          nt = nt + 1
+          READ ( finput, *, IOSTAT = ierrn ) time_surf(nt), shf_surf(nt),      &
+                                             qsws_surf(nt), pt_surf(nt),       &
+                                             q_surf(nt), p_surf(nt)
+
+          IF ( ierrn /= 0 )  THEN
+            WRITE ( message_string, * ) 'No time dependent surface ' //        &
+                              'variables in & LSF_DATA for end of run found'
+
+             CALL message( 'ls_forcing', 'PA0363', 1, 2, 0, 6, 0 )
           ENDIF
-!
-!--       Allocate arrays for reading boundary values. Arrays will incorporate 2 
-!--       time levels in order to interpolate in between. 
-          IF ( bc_dirichlet_l )  THEN
-             ALLOCATE( nest_offl%u_left(0:1,nzb+1:nzt,nys:nyn)  )
-             ALLOCATE( nest_offl%v_left(0:1,nzb+1:nzt,nysv:nyn) )
-             ALLOCATE( nest_offl%w_left(0:1,nzb+1:nzt-1,nys:nyn) )
-             IF ( humidity )      ALLOCATE( nest_offl%q_left(0:1,nzb+1:nzt,nys:nyn)  )
-             IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_left(0:1,nzb+1:nzt,nys:nyn) )
-          ENDIF
-          IF ( bc_dirichlet_r )  THEN
-             ALLOCATE( nest_offl%u_right(0:1,nzb+1:nzt,nys:nyn)  )
-             ALLOCATE( nest_offl%v_right(0:1,nzb+1:nzt,nysv:nyn) )
-             ALLOCATE( nest_offl%w_right(0:1,nzb+1:nzt-1,nys:nyn) )
-             IF ( humidity )      ALLOCATE( nest_offl%q_right(0:1,nzb+1:nzt,nys:nyn)  )
-             IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_right(0:1,nzb+1:nzt,nys:nyn) )
-          ENDIF
-          IF ( bc_dirichlet_n )  THEN
-             ALLOCATE( nest_offl%u_north(0:1,nzb+1:nzt,nxlu:nxr) )
-             ALLOCATE( nest_offl%v_north(0:1,nzb+1:nzt,nxl:nxr)  )
-             ALLOCATE( nest_offl%w_north(0:1,nzb+1:nzt-1,nxl:nxr) )
-             IF ( humidity )      ALLOCATE( nest_offl%q_north(0:1,nzb+1:nzt,nxl:nxr)  )
-             IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_north(0:1,nzb+1:nzt,nxl:nxr) )
-          ENDIF
-          IF ( bc_dirichlet_s )  THEN
-             ALLOCATE( nest_offl%u_south(0:1,nzb+1:nzt,nxlu:nxr) )
-             ALLOCATE( nest_offl%v_south(0:1,nzb+1:nzt,nxl:nxr)  )
-             ALLOCATE( nest_offl%w_south(0:1,nzb+1:nzt-1,nxl:nxr)    )
-             IF ( humidity )      ALLOCATE( nest_offl%q_south(0:1,nzb+1:nzt,nxl:nxr)  )
-             IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_south(0:1,nzb+1:nzt,nxl:nxr) )
-          ENDIF
-          
-          ALLOCATE( nest_offl%u_top(0:1,nys:nyn,nxlu:nxr) )
-          ALLOCATE( nest_offl%v_top(0:1,nysv:nyn,nxl:nxr) )
-          ALLOCATE( nest_offl%w_top(0:1,nys:nyn,nxl:nxr)  )
-          IF ( humidity )      ALLOCATE( nest_offl%q_top(0:1,nys:nyn,nxl:nxr)  )
-          IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_top(0:1,nys:nyn,nxl:nxr) )
-
-!
-!--       Read COSMO data at lateral and top boundaries
-          CALL netcdf_data_input_lsf
-!
-!--       Write COSMO data at lateral and top boundaries
-          CALL lsf_nesting_offline
-!
-!--       After 3D data is initialized, ensure mass conservation
-          CALL lsf_nesting_offline_mass_conservation
-!
-!--       Initialize surface pressure. Please note, time-dependent surface
-!--       pressure would require changes in anelastic approximation and 
-!--       treatment of fluxes. 
-!--       For the moment, comment this out! 
-!         surface_pressure = nest_offl%surface_pressure(0)
-
-       ELSE
-
-          ALLOCATE( p_surf(0:nlsf), pt_surf(0:nlsf), q_surf(0:nlsf),           &
-                    qsws_surf(0:nlsf), shf_surf(0:nlsf),                       &
-                    td_lsa_lpt(nzb:nzt+1,0:nlsf), td_lsa_q(nzb:nzt+1,0:nlsf),  &
-                    td_sub_lpt(nzb:nzt+1,0:nlsf), td_sub_q(nzb:nzt+1,0:nlsf),  &
-                    time_vert(0:nlsf), time_surf(0:nlsf),                      &
-                    ug_vert(nzb:nzt+1,0:nlsf), vg_vert(nzb:nzt+1,0:nlsf),      &
-                    wsubs_vert(nzb:nzt+1,0:nlsf) )
-
-          p_surf = 0.0_wp; pt_surf = 0.0_wp; q_surf = 0.0_wp; qsws_surf = 0.0_wp
-          shf_surf = 0.0_wp; time_vert = 0.0_wp; td_lsa_lpt = 0.0_wp
-          td_lsa_q = 0.0_wp; td_sub_lpt = 0.0_wp; td_sub_q = 0.0_wp
-          time_surf = 0.0_wp; ug_vert = 0.0_wp; vg_vert = 0.0_wp
-          wsubs_vert = 0.0_wp
-
-!
-!--       Array for storing large scale forcing and nudging tendencies at each 
-!--       timestep for data output
-          ALLOCATE( sums_ls_l(nzb:nzt+1,0:7) )
-          sums_ls_l = 0.0_wp
-
-          ngp_sums_ls = (nz+2)*6
-
-          OPEN ( finput, FILE='LSF_DATA', STATUS='OLD', &
-                 FORM='FORMATTED', IOSTAT=ierrn )
-
-          IF ( ierrn /= 0 )  THEN
-             message_string = 'file LSF_DATA does not exist'
-             CALL message( 'ls_forcing', 'PA0368', 1, 2, 0, 6, 0 )
-          ENDIF
-
-          ierrn = 0
-!
-!--       First three lines of LSF_DATA contain header
-          READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
-          READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
-          READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
-
+       ENDDO
+
+       IF ( time_surf(1) > end_time )  THEN
+          WRITE ( message_string, * ) 'Time dependent surface variables in ' //&
+                                      '&LSF_DATA set in after end of ' ,       &
+                                      'simulation - lsf_surf is set to FALSE'
+          CALL message( 'ls_forcing', 'PA0371', 0, 0, 0, 6, 0 )
+          lsf_surf = .FALSE.
+       ENDIF
+
+!
+!--    Go to the end of the list with surface variables
+       DO WHILE ( ierrn == 0 )
+          READ ( finput, *, IOSTAT = ierrn ) r_dummy
+       ENDDO
+
+!
+!--    Profiles of ug, vg and w_subs are read in (large scale forcing)
+
+       nt = 0
+       DO WHILE ( time_vert(nt) < end_time )
+          nt = nt + 1
+          hash = "#"
+          ierrn = 1 ! not zero
+!
+!--       Search for the next line consisting of "# time", 
+!--       from there onwards the profiles will be read
+          DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) ) 
+             READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(nt)
+             IF ( ierrn < 0 )  THEN 
+                WRITE( message_string, * ) 'No time dependent vertical profiles',&
+                                 ' in & LSF_DATA for end of run found'
+                CALL message( 'ls_forcing', 'PA0372', 1, 2, 0, 6, 0 )
+             ENDIF
+          ENDDO
+
+          IF ( nt == 1 .AND. time_vert(nt) > end_time ) EXIT
+
+          READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert,  &
+                                           lowwsubs_vert, low_td_lsa_lpt,      &
+                                           low_td_lsa_q, low_td_sub_lpt,       &
+                                           low_td_sub_q
           IF ( ierrn /= 0 )  THEN
              message_string = 'errors in file LSF_DATA'
@@ -1075,146 +591,80 @@
           ENDIF
 
-!
-!--       Surface values are read in
-          nt     = 0
-          ierrn = 0
-
-          DO WHILE ( time_surf(nt) < end_time )
-             nt = nt + 1
-             READ ( finput, *, IOSTAT = ierrn ) time_surf(nt), shf_surf(nt),   &
-                                                qsws_surf(nt), pt_surf(nt),    &
-                                                q_surf(nt), p_surf(nt)
-
-             IF ( ierrn /= 0 )  THEN
-               WRITE ( message_string, * ) 'No time dependent surface ' //     &
-                                 'variables in & LSF_DATA for end of run found'
-
-                CALL message( 'ls_forcing', 'PA0363', 1, 2, 0, 6, 0 )
+          READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,            &
+                                           highvg_vert, highwsubs_vert,        &
+                                           high_td_lsa_lpt, high_td_lsa_q,     &
+                                           high_td_sub_lpt, high_td_sub_q
+       
+          IF ( ierrn /= 0 )  THEN
+             message_string = 'errors in file LSF_DATA'
+             CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
+          ENDIF
+
+
+          DO  k = nzb, nzt+1
+             IF ( highheight < zu(k) )  THEN
+                lowheight      = highheight
+                lowug_vert     = highug_vert
+                lowvg_vert     = highvg_vert
+                lowwsubs_vert  = highwsubs_vert
+                low_td_lsa_lpt = high_td_lsa_lpt
+                low_td_lsa_q   = high_td_lsa_q
+                low_td_sub_lpt = high_td_sub_lpt
+                low_td_sub_q   = high_td_sub_q
+
+                ierrn = 0
+                READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,      &
+                                                 highvg_vert, highwsubs_vert,  &
+                                                 high_td_lsa_lpt,              &
+                                                 high_td_lsa_q,                &
+                                                 high_td_sub_lpt, high_td_sub_q
+
+                IF ( ierrn /= 0 )  THEN
+                   WRITE( message_string, * ) 'zu(',k,') = ', zu(k), 'm ',     &
+                        'is higher than the maximum height in LSF_DATA ',      &
+                        'which is ', lowheight, 'm. Interpolation on PALM ',   &
+                        'grid is not possible.'
+                   CALL message( 'ls_forcing', 'PA0395', 1, 2, 0, 6, 0 )
+                ENDIF
+
              ENDIF
+
+!
+!--          Interpolation of prescribed profiles in space 
+             fac = (highheight-zu(k))/(highheight - lowheight)
+
+             ug_vert(k,nt)    = fac * lowug_vert                               &
+                                + ( 1.0_wp - fac ) * highug_vert
+             vg_vert(k,nt)    = fac * lowvg_vert                               &
+                                + ( 1.0_wp - fac ) * highvg_vert
+             wsubs_vert(k,nt) = fac * lowwsubs_vert                            &
+                                + ( 1.0_wp - fac ) * highwsubs_vert
+
+             td_lsa_lpt(k,nt) = fac * low_td_lsa_lpt                           &
+                                + ( 1.0_wp - fac ) * high_td_lsa_lpt
+             td_lsa_q(k,nt)   = fac * low_td_lsa_q                             &
+                                + ( 1.0_wp - fac ) * high_td_lsa_q
+             td_sub_lpt(k,nt) = fac * low_td_sub_lpt                           &
+                                + ( 1.0_wp - fac ) * high_td_sub_lpt
+             td_sub_q(k,nt)   = fac * low_td_sub_q                             &
+                                + ( 1.0_wp - fac ) * high_td_sub_q
+
           ENDDO
 
-          IF ( time_surf(1) > end_time )  THEN
-             WRITE ( message_string, * ) 'Time dependent surface variables in ' // &
-                                         '&LSF_DATA set in after end of ' ,        &
-                                         'simulation - lsf_surf is set to FALSE'
-             CALL message( 'ls_forcing', 'PA0371', 0, 0, 0, 6, 0 )
-             lsf_surf = .FALSE.
-          ENDIF
-
-!
-!--       Go to the end of the list with surface variables
-          DO WHILE ( ierrn == 0 )
-             READ ( finput, *, IOSTAT = ierrn ) r_dummy
-          ENDDO
-
-!
-!--       Profiles of ug, vg and w_subs are read in (large scale forcing)
-
-          nt = 0
-          DO WHILE ( time_vert(nt) < end_time )
-             nt = nt + 1
-             hash = "#"
-             ierrn = 1 ! not zero
-!
-!--          Search for the next line consisting of "# time", 
-!--          from there onwards the profiles will be read
-             DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) ) 
-                READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(nt)
-                IF ( ierrn < 0 )  THEN 
-                   WRITE( message_string, * ) 'No time dependent vertical profiles',&
-                                    ' in & LSF_DATA for end of run found'
-                   CALL message( 'ls_forcing', 'PA0372', 1, 2, 0, 6, 0 )
-                ENDIF
-             ENDDO
-
-             IF ( nt == 1 .AND. time_vert(nt) > end_time ) EXIT
-
-             READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert,&
-                                              lowwsubs_vert, low_td_lsa_lpt,    &
-                                              low_td_lsa_q, low_td_sub_lpt,     &
-                                              low_td_sub_q
-             IF ( ierrn /= 0 )  THEN
-                message_string = 'errors in file LSF_DATA'
-                CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
-             ENDIF
-
-             READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,         &
-                                              highvg_vert, highwsubs_vert,     &
-                                              high_td_lsa_lpt, high_td_lsa_q,  &
-                                              high_td_sub_lpt, high_td_sub_q
-         
-             IF ( ierrn /= 0 )  THEN
-                message_string = 'errors in file LSF_DATA'
-                CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
-             ENDIF
-
-
-             DO  k = nzb, nzt+1
-                IF ( highheight < zu(k) )  THEN
-                   lowheight      = highheight
-                   lowug_vert     = highug_vert
-                   lowvg_vert     = highvg_vert
-                   lowwsubs_vert  = highwsubs_vert
-                   low_td_lsa_lpt = high_td_lsa_lpt
-                   low_td_lsa_q   = high_td_lsa_q
-                   low_td_sub_lpt = high_td_sub_lpt
-                   low_td_sub_q   = high_td_sub_q
-
-                   ierrn = 0
-                   READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,    &
-                                                    highvg_vert, highwsubs_vert,&
-                                                    high_td_lsa_lpt,            &
-                                                    high_td_lsa_q,              &
-                                                    high_td_sub_lpt, high_td_sub_q
-
-                   IF ( ierrn /= 0 )  THEN
-                      WRITE( message_string, * ) 'zu(',k,') = ', zu(k), 'm ',  &
-                           'is higher than the maximum height in LSF_DATA ',   &
-                           'which is ', lowheight, 'm. Interpolation on PALM ',&
-                           'grid is not possible.'
-                      CALL message( 'ls_forcing', 'PA0395', 1, 2, 0, 6, 0 )
-                   ENDIF
-
-                ENDIF
-
-!
-!--             Interpolation of prescribed profiles in space 
-                fac = (highheight-zu(k))/(highheight - lowheight)
-
-                ug_vert(k,nt)    = fac * lowug_vert                            &
-                                   + ( 1.0_wp - fac ) * highug_vert
-                vg_vert(k,nt)    = fac * lowvg_vert                            &
-                                   + ( 1.0_wp - fac ) * highvg_vert
-                wsubs_vert(k,nt) = fac * lowwsubs_vert                         &
-                                   + ( 1.0_wp - fac ) * highwsubs_vert
-
-                td_lsa_lpt(k,nt) = fac * low_td_lsa_lpt                        &
-                                   + ( 1.0_wp - fac ) * high_td_lsa_lpt
-                td_lsa_q(k,nt)   = fac * low_td_lsa_q                          &
-                                   + ( 1.0_wp - fac ) * high_td_lsa_q
-                td_sub_lpt(k,nt) = fac * low_td_sub_lpt                        &
-                                   + ( 1.0_wp - fac ) * high_td_sub_lpt
-                td_sub_q(k,nt)   = fac * low_td_sub_q                          &
-                                   + ( 1.0_wp - fac ) * high_td_sub_q
-
-             ENDDO
-
-          ENDDO 
-
-!
-!--       Large scale vertical velocity has to be zero at the surface
-          wsubs_vert(nzb,:) = 0.0_wp
+       ENDDO 
+
+!
+!--    Large scale vertical velocity has to be zero at the surface
+       wsubs_vert(nzb,:) = 0.0_wp
     
-          IF ( time_vert(1) > end_time )  THEN
-             WRITE ( message_string, * ) 'Time dependent large scale profile ',&
-                                'forcing from&LSF_DATA sets in after end of ' ,&
-                                'simulation - lsf_vert is set to FALSE'
-             CALL message( 'ls_forcing', 'PA0373', 0, 0, 0, 6, 0 )
-             lsf_vert = .FALSE.
-          ENDIF
-
-          CLOSE( finput )
-
-       ENDIF
+       IF ( time_vert(1) > end_time )  THEN
+          WRITE ( message_string, * ) 'Time dependent large scale profile ',   &
+                             'forcing from&LSF_DATA sets in after end of ' ,   &
+                             'simulation - lsf_vert is set to FALSE'
+          CALL message( 'ls_forcing', 'PA0373', 0, 0, 0, 6, 0 )
+          lsf_vert = .FALSE.
+       ENDIF
+
+       CLOSE( finput )
 
     END SUBROUTINE lsf_init
Index: /palm/trunk/SOURCE/nesting_offl_mod.f90
===================================================================
--- /palm/trunk/SOURCE/nesting_offl_mod.f90	(revision 3347)
+++ /palm/trunk/SOURCE/nesting_offl_mod.f90	(revision 3347)
@@ -0,0 +1,1090 @@
+!> @file nesting_offl_mod.f90
+!------------------------------------------------------------------------------!
+! This file is part of the PALM model system.
+!
+! PALM is free software: you can redistribute it and/or modify it under the
+! terms of the GNU General Public License as published by the Free Software
+! Foundation, either version 3 of the License, or (at your option) any later
+! version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with
+! PALM. If not, see <http://www.gnu.org/licenses/>.
+!
+! Copyright 1997-2018 Leibniz Universitaet Hannover
+!------------------------------------------------------------------------------!
+!
+! Current revisions:
+! ------------------
+! 
+! 
+! Former revisions:
+! -----------------
+!
+! Initial Revision: 
+! - separate offline nesting from large_scale_nudging_mod
+! - revise offline nesting, adjust for usage of synthetic turbulence generator
+! - adjust Rayleigh damping depending on the time-depending atmospheric
+!   conditions
+! 
+! 
+! Description:
+! ------------
+!> Offline nesting in larger-scale models. Boundary conditions for the simulation
+!> are read from NetCDF file and are prescribed onto the respective arrays.
+!> Further, a mass-flux correction is performed to maintain the mass balance. 
+!--------------------------------------------------------------------------------!
+ MODULE nesting_offl_mod
+
+    USE arrays_3d,                                                             &
+        ONLY:  dzw, e, diss, pt, pt_init, q, q_init, s, u, u_init, ug, v,      &
+               v_init, vg, w, zu, zw          
+
+    USE control_parameters,                                                    &
+        ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, &
+               dt_3d, dz, constant_diffusion, humidity, message_string,        &
+               nesting_offline, neutral, passive_scalar, rans_mode, rans_tke_e,&
+               time_since_reference_point, volume_flow
+               
+    USE cpulog,                                                                &
+        ONLY:  cpu_log, log_point
+
+    USE grid_variables
+
+    USE indices,                                                               &
+        ONLY:  nbgp, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nys,                  &
+               nysv, nysg, nyn, nyng, nzb, nz, nzt, wall_flags_0
+
+    USE kinds
+
+    USE netcdf_data_input_mod,                                                 &
+        ONLY:  nest_offl
+        
+    USE pegrid
+
+    REAL(wp) ::  zi_ribulk = 0.0_wp  !< boundary-layer depth according to bulk Richardson criterion, i.e. the height where Ri_bulk exceeds the critical
+                                     !< bulk Richardson number of 0.25
+    
+    SAVE
+    PRIVATE
+!
+!-- Public subroutines
+    PUBLIC nesting_offl_bc, nesting_offl_header, nesting_offl_init,            &
+           nesting_offl_mass_conservation, nesting_offl_parin 
+!
+!-- Public variables
+    PUBLIC zi_ribulk   
+
+    INTERFACE nesting_offl_bc
+       MODULE PROCEDURE nesting_offl_bc
+    END INTERFACE nesting_offl_bc
+    
+    INTERFACE nesting_offl_header
+       MODULE PROCEDURE nesting_offl_header
+    END INTERFACE nesting_offl_header
+    
+    INTERFACE nesting_offl_init
+       MODULE PROCEDURE nesting_offl_init
+    END INTERFACE nesting_offl_init
+           
+    INTERFACE nesting_offl_mass_conservation
+       MODULE PROCEDURE nesting_offl_mass_conservation
+    END INTERFACE nesting_offl_mass_conservation
+    
+    INTERFACE nesting_offl_parin
+       MODULE PROCEDURE nesting_offl_parin
+    END INTERFACE nesting_offl_parin
+
+ CONTAINS
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> In this subroutine a constant mass within the model domain is guaranteed. 
+!> Larger-scale models may be based on a compressible equation system, which is
+!> not consistent with PALMs incompressible equation system. In order to avoid
+!> a decrease or increase of mass during the simulation, non-divergent flow
+!> through the lateral and top boundaries is compensated by the vertical wind 
+!> component at the top boundary. 
+!------------------------------------------------------------------------------!
+    SUBROUTINE nesting_offl_mass_conservation
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i !< grid index in x-direction
+       INTEGER(iwp) ::  j !< grid index in y-direction
+       INTEGER(iwp) ::  k !< grid index in z-direction
+
+       REAL(wp) ::  d_area_t                        !< inverse of the total area of the horizontal model domain
+       REAL(wp) ::  w_correct                       !< vertical velocity increment required to compensate non-divergent flow through the boundaries
+       REAL(wp), DIMENSION(1:3) ::  volume_flow_l   !< local volume flow
+
+       
+       CALL  cpu_log( log_point(58), 'offline nesting', 'start' )
+       
+       volume_flow   = 0.0_wp
+       volume_flow_l = 0.0_wp
+
+       d_area_t = 1.0_wp / ( ( nx + 1 ) * dx * ( ny + 1 ) * dy )
+
+       IF ( bc_dirichlet_l )  THEN
+          i = nxl
+          DO  j = nys, nyn
+             DO  k = nzb+1, nzt
+                volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) * dy   &
+                                     * MERGE( 1.0_wp, 0.0_wp,                  &
+                                              BTEST( wall_flags_0(k,j,i), 1 ) )
+             ENDDO
+          ENDDO
+       ENDIF
+       IF ( bc_dirichlet_r )  THEN
+          i = nxr+1
+          DO  j = nys, nyn
+             DO  k = nzb+1, nzt
+                volume_flow_l(1) = volume_flow_l(1) - u(k,j,i) * dzw(k) * dy   &
+                                     * MERGE( 1.0_wp, 0.0_wp,                  &
+                                              BTEST( wall_flags_0(k,j,i), 1 ) )
+             ENDDO
+          ENDDO
+       ENDIF
+       IF ( bc_dirichlet_s )  THEN
+          j = nys
+          DO  i = nxl, nxr
+             DO  k = nzb+1, nzt
+                volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) * dx   &
+                                     * MERGE( 1.0_wp, 0.0_wp,                  &
+                                              BTEST( wall_flags_0(k,j,i), 2 ) )
+             ENDDO
+          ENDDO
+       ENDIF
+       IF ( bc_dirichlet_n )  THEN
+          j = nyn+1
+          DO  i = nxl, nxr
+             DO  k = nzb+1, nzt
+                volume_flow_l(2) = volume_flow_l(2) - v(k,j,i) * dzw(k) * dx   &
+                                     * MERGE( 1.0_wp, 0.0_wp,                  &
+                                              BTEST( wall_flags_0(k,j,i), 2 ) )
+             ENDDO
+          ENDDO
+       ENDIF
+!
+!--    Top boundary
+       k = nzt
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dx * dy
+          ENDDO
+       ENDDO
+
+#if defined( __parallel )
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, 3, MPI_REAL, MPI_SUM,      &
+                           comm2d, ierr )
+#else
+       volume_flow = volume_flow_l
+#endif
+
+       w_correct = SUM( volume_flow ) * d_area_t
+
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzt, nzt + 1
+                w(k,j,i) = w(k,j,i) + w_correct
+             ENDDO
+          ENDDO
+       ENDDO
+       
+       CALL  cpu_log( log_point(58), 'offline nesting', 'stop' )
+
+    END SUBROUTINE nesting_offl_mass_conservation
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Set the lateral and top boundary conditions in case the PALM domain is 
+!> nested offline in a mesoscale model. Further, average boundary data and
+!> determine mean profiles, further used for correct damping in the sponge
+!> layer. 
+!------------------------------------------------------------------------------!
+    SUBROUTINE nesting_offl_bc                     
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i !< running index x-direction
+       INTEGER(iwp) ::  j !< running index y-direction
+       INTEGER(iwp) ::  k !< running index z-direction
+
+       REAL(wp) ::  fac_dt   !< interpolation factor
+       
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  pt_ref   !< reference profile for potential temperature
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  pt_ref_l !< reference profile for potential temperature on subdomain
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  q_ref    !< reference profile for mixing ratio on subdomain
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  q_ref_l  !< reference profile for mixing ratio on subdomain
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  u_ref    !< reference profile for u-component on subdomain
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  u_ref_l  !< reference profile for u-component on subdomain
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  v_ref    !< reference profile for v-component on subdomain
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  v_ref_l  !< reference profile for v-component on subdomain
+       
+       CALL  cpu_log( log_point(58), 'offline nesting', 'start' )      
+!
+!--    Set mean profiles, derived from boundary data, to zero
+       pt_ref   = 0.0_wp
+       q_ref    = 0.0_wp
+       u_ref    = 0.0_wp
+       v_ref    = 0.0_wp
+       
+       pt_ref_l = 0.0_wp
+       q_ref_l  = 0.0_wp
+       u_ref_l  = 0.0_wp
+       v_ref_l  = 0.0_wp
+!
+!--    Determine interpolation factor and limit it to 1. This is because
+!--    t+dt can slightly exceed time(tind_p) before boundary data is updated 
+!--    again. 
+       fac_dt = ( time_since_reference_point - nest_offl%time(nest_offl%tind)  &
+                + dt_3d ) /                                                    &
+           ( nest_offl%time(nest_offl%tind_p) - nest_offl%time(nest_offl%tind) )
+       fac_dt = MIN( 1.0_wp, fac_dt )
+!
+!--    Set boundary conditions of u-, v-, w-component, as well as q, and pt.
+!--    Note, boundary values at the left boundary: i=-1 (v,w,pt,q) and 
+!--    i=0 (u), at the right boundary: i=nxr+1 (all), at the south boundary:
+!--    j=-1 (u,w,pt,q) and j=0 (v), at the north boundary: j=nyn+1 (all).
+!--    Please note, at the left (for u) and south (for v) boundary, values 
+!--    for u and v are set also at i/j=-1, since these values are used in 
+!--    boundary_conditions() to restore prognostic values.
+!--    Further, sum up data to calculate mean profiles from boundary data, 
+!--    used for Rayleigh damping. 
+       IF ( bc_dirichlet_l )  THEN
+
+          DO  j = nys, nyn
+             DO  k = nzb+1, nzt
+                u(k,j,0) = interpolate_in_time( nest_offl%u_left(0,k,j),       &
+                                                nest_offl%u_left(1,k,j),       &
+                                                fac_dt ) *                     &
+                             MERGE( 1.0_wp, 0.0_wp,                            &
+                                    BTEST( wall_flags_0(k,j,0), 1 ) )
+                u(k,j,-1) = u(k,j,0)
+             ENDDO
+             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,j,0)
+          ENDDO
+
+          DO  j = nys, nyn
+             DO  k = nzb+1, nzt-1
+                w(k,j,-1) = interpolate_in_time( nest_offl%w_left(0,k,j),      &
+                                                 nest_offl%w_left(1,k,j),      &
+                                                 fac_dt ) *                    &
+                            MERGE( 1.0_wp, 0.0_wp,                             &
+                                   BTEST( wall_flags_0(k,j,-1), 3 ) )
+             ENDDO
+          ENDDO
+
+          DO  j = nysv, nyn
+             DO  k = nzb+1, nzt
+                v(k,j,-1) = interpolate_in_time( nest_offl%v_left(0,k,j),      &
+                                                 nest_offl%v_left(1,k,j),      &
+                                                 fac_dt ) *                    &
+                               MERGE( 1.0_wp, 0.0_wp,                          &
+                                      BTEST( wall_flags_0(k,j,-1), 2 ) )
+             ENDDO
+             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,j,-1)
+          ENDDO
+
+          IF ( .NOT. neutral )  THEN
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt
+                   pt(k,j,-1) = interpolate_in_time( nest_offl%pt_left(0,k,j), &
+                                                     nest_offl%pt_left(1,k,j), &
+                                                     fac_dt )
+ 
+                ENDDO
+                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,j,-1)
+             ENDDO
+          ENDIF
+
+          IF ( humidity )  THEN
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt
+                   q(k,j,-1) = interpolate_in_time( nest_offl%q_left(0,k,j),   &
+                                                    nest_offl%q_left(1,k,j),   &
+                                                    fac_dt )
+ 
+                ENDDO
+                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,j,-1)
+             ENDDO
+          ENDIF
+
+       ENDIF
+
+       IF ( bc_dirichlet_r )  THEN
+
+          DO  j = nys, nyn
+             DO  k = nzb+1, nzt
+                u(k,j,nxr+1) = interpolate_in_time( nest_offl%u_right(0,k,j),  &
+                                                    nest_offl%u_right(1,k,j),  &
+                                                    fac_dt ) *                 &
+                             MERGE( 1.0_wp, 0.0_wp,                            &
+                                    BTEST( wall_flags_0(k,j,nxr+1), 1 ) )
+             ENDDO
+             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,j,nxr+1)
+          ENDDO
+          DO  j = nys, nyn
+             DO  k = nzb+1, nzt-1
+                w(k,j,nxr+1) = interpolate_in_time( nest_offl%w_right(0,k,j),  &
+                                                    nest_offl%w_right(1,k,j),  &
+                                                    fac_dt ) *                 &
+                             MERGE( 1.0_wp, 0.0_wp,                            &
+                                    BTEST( wall_flags_0(k,j,nxr+1), 3 ) )
+             ENDDO
+          ENDDO
+
+          DO  j = nysv, nyn
+             DO  k = nzb+1, nzt
+                v(k,j,nxr+1) = interpolate_in_time( nest_offl%v_right(0,k,j),  &
+                                                    nest_offl%v_right(1,k,j),  &
+                                                    fac_dt ) *                 &
+                             MERGE( 1.0_wp, 0.0_wp,                            &
+                                    BTEST( wall_flags_0(k,j,nxr+1), 2 ) )
+             ENDDO
+             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,j,nxr+1)
+          ENDDO
+
+          IF ( .NOT. neutral )  THEN
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt
+                   pt(k,j,nxr+1) = interpolate_in_time(                        &
+                                                  nest_offl%pt_right(0,k,j),   &
+                                                  nest_offl%pt_right(1,k,j),   &
+                                                  fac_dt )
+                ENDDO
+                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,j,nxr+1)
+             ENDDO
+          ENDIF
+
+          IF ( humidity )  THEN
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt
+                   q(k,j,nxr+1) = interpolate_in_time(                         &
+                                                  nest_offl%q_right(0,k,j),    &
+                                                  nest_offl%q_right(1,k,j),    &
+                                                  fac_dt ) 
+ 
+                ENDDO
+                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,j,nxr+1)
+             ENDDO
+          ENDIF
+
+       ENDIF
+
+       IF ( bc_dirichlet_s )  THEN
+
+          DO  i = nxl, nxr
+             DO  k = nzb+1, nzt
+                v(k,0,i) = interpolate_in_time( nest_offl%v_south(0,k,i),      &
+                                                nest_offl%v_south(1,k,i),      &
+                                                fac_dt ) *                     &
+                           MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(k,0,i), 2 ) )
+                v(k,-1,i) = v(k,0,i) 
+             ENDDO
+             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,0,i)
+          ENDDO
+
+          DO  i = nxl, nxr
+             DO  k = nzb+1, nzt-1
+                w(k,-1,i) = interpolate_in_time( nest_offl%w_south(0,k,i),     &
+                                                 nest_offl%w_south(1,k,i),     &
+                                                 fac_dt ) *                    &
+                           MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(k,-1,i), 3 ) )
+             ENDDO
+          ENDDO
+
+          DO  i = nxlu, nxr
+             DO  k = nzb+1, nzt
+                u(k,-1,i) = interpolate_in_time( nest_offl%u_south(0,k,i),     &
+                                                 nest_offl%u_south(1,k,i),     &
+                                                 fac_dt ) *                    &
+                           MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(k,-1,i), 1 ) )
+             ENDDO
+             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,-1,i)
+          ENDDO
+
+          IF ( .NOT. neutral )  THEN
+             DO  i = nxl, nxr
+                DO  k = nzb+1, nzt
+                   pt(k,-1,i) = interpolate_in_time(                           &
+                                                 nest_offl%pt_south(0,k,i),    &
+                                                 nest_offl%pt_south(1,k,i),    &
+                                                 fac_dt )
+ 
+                ENDDO
+                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,-1,i)
+             ENDDO
+          ENDIF
+
+          IF ( humidity )  THEN
+             DO  i = nxl, nxr
+                DO  k = nzb+1, nzt
+                   q(k,-1,i) = interpolate_in_time(                            &
+                                                 nest_offl%q_south(0,k,i),     &
+                                                 nest_offl%q_south(1,k,i),     &
+                                                 fac_dt )
+ 
+                ENDDO
+                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,-1,i)
+             ENDDO
+          ENDIF
+
+       ENDIF
+
+       IF ( bc_dirichlet_n )  THEN
+
+          DO  i = nxl, nxr
+             DO  k = nzb+1, nzt
+                v(k,nyn+1,i) = interpolate_in_time( nest_offl%v_north(0,k,i),  &
+                                                    nest_offl%v_north(1,k,i),  &
+                                                    fac_dt ) *                 &
+                               MERGE( 1.0_wp, 0.0_wp,                          &
+                                      BTEST( wall_flags_0(k,nyn+1,i), 2 ) )
+             ENDDO
+             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,nyn+1,i)
+          ENDDO
+          DO  i = nxl, nxr
+             DO  k = nzb+1, nzt-1
+                w(k,nyn+1,i) = interpolate_in_time( nest_offl%w_north(0,k,i),  &
+                                                    nest_offl%w_north(1,k,i),  &
+                                                    fac_dt ) *                 &
+                               MERGE( 1.0_wp, 0.0_wp,                          &
+                                      BTEST( wall_flags_0(k,nyn+1,i), 3 ) )
+             ENDDO
+          ENDDO
+
+          DO  i = nxlu, nxr
+             DO  k = nzb+1, nzt
+                u(k,nyn+1,i) = interpolate_in_time( nest_offl%u_north(0,k,i),  &
+                                                    nest_offl%u_north(1,k,i),  &
+                                                    fac_dt ) *                 &
+                               MERGE( 1.0_wp, 0.0_wp,                          &
+                                      BTEST( wall_flags_0(k,nyn+1,i), 1 ) )
+
+             ENDDO
+             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,nyn+1,i)
+          ENDDO
+
+          IF ( .NOT. neutral )  THEN
+             DO  i = nxl, nxr
+                DO  k = nzb+1, nzt
+                   pt(k,nyn+1,i) = interpolate_in_time(                        &
+                                                    nest_offl%pt_north(0,k,i), &
+                                                    nest_offl%pt_north(1,k,i), &
+                                                    fac_dt )
+ 
+                ENDDO
+                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,nyn+1,i)
+             ENDDO
+          ENDIF
+
+          IF ( humidity )  THEN
+             DO  i = nxl, nxr
+                DO  k = nzb+1, nzt
+                   q(k,nyn+1,i) = interpolate_in_time(                         &
+                                                    nest_offl%q_north(0,k,i),  &
+                                                    nest_offl%q_north(1,k,i),  &
+                                                    fac_dt )
+ 
+                ENDDO
+                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,nyn+1,i)
+             ENDDO
+          ENDIF
+
+       ENDIF
+!
+!--    Top boundary
+       DO  i = nxlu, nxr
+          DO  j = nys, nyn
+             u(nzt+1,j,i) = interpolate_in_time( nest_offl%u_top(0,j,i),       &
+                                                 nest_offl%u_top(1,j,i),       &
+                                                 fac_dt ) *                    &
+                           MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
+             u_ref_l(nzt+1) = u_ref_l(nzt+1) + u(nzt+1,j,i)
+          ENDDO
+       ENDDO
+
+       DO  i = nxl, nxr
+          DO  j = nysv, nyn
+             v(nzt+1,j,i) = interpolate_in_time( nest_offl%v_top(0,j,i),       &
+                                                 nest_offl%v_top(1,j,i),       &
+                                                 fac_dt ) *                    &
+                           MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
+             v_ref_l(nzt+1) = v_ref_l(nzt+1) + v(nzt+1,j,i)
+          ENDDO
+       ENDDO
+
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             w(nzt,j,i) = interpolate_in_time( nest_offl%w_top(0,j,i),         &
+                                               nest_offl%w_top(1,j,i),         &
+                                               fac_dt ) *                      &
+                           MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(nzt,j,i), 3 ) )
+             w(nzt+1,j,i) = w(nzt,j,i)
+          ENDDO
+       ENDDO
+
+
+       IF ( .NOT. neutral )  THEN
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                pt(nzt+1,j,i) = interpolate_in_time( nest_offl%pt_top(0,j,i),  &
+                                                     nest_offl%pt_top(1,j,i),  &
+                                                     fac_dt )
+                pt_ref_l(nzt+1) = pt_ref_l(nzt+1) + pt(nzt+1,j,i)
+             ENDDO
+          ENDDO
+       ENDIF
+
+       IF ( humidity )  THEN
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                q(nzt+1,j,i) = interpolate_in_time( nest_offl%q_top(0,j,i),    &
+                                                    nest_offl%q_top(1,j,i),    &
+                                                    fac_dt )
+                q_ref_l(nzt+1) = q_ref_l(nzt+1) + q(nzt+1,j,i)
+             ENDDO
+          ENDDO
+       ENDIF
+!
+!--    Moreover, set Neumann boundary condition for subgrid-scale TKE, 
+!--    passive scalar, dissipation, and chemical species if required
+       IF ( rans_mode  .AND.  rans_tke_e )  THEN
+          IF (  bc_dirichlet_l )  diss(:,:,nxl-1) = diss(:,:,nxl)
+          IF (  bc_dirichlet_r )  diss(:,:,nxr+1) = diss(:,:,nxr)
+          IF (  bc_dirichlet_s )  diss(:,nys-1,:) = diss(:,nys,:)
+          IF (  bc_dirichlet_n )  diss(:,nyn+1,:) = diss(:,nyn,:)
+       ENDIF
+       IF ( .NOT. constant_diffusion )  THEN
+          IF (  bc_dirichlet_l )  e(:,:,nxl-1) = e(:,:,nxl)
+          IF (  bc_dirichlet_r )  e(:,:,nxr+1) = e(:,:,nxr)
+          IF (  bc_dirichlet_s )  e(:,nys-1,:) = e(:,nys,:)
+          IF (  bc_dirichlet_n )  e(:,nyn+1,:) = e(:,nyn,:)
+          e(nzt+1,:,:) = e(nzt,:,:)
+       ENDIF
+       IF ( passive_scalar )  THEN
+          IF (  bc_dirichlet_l )  s(:,:,nxl-1) = s(:,:,nxl)
+          IF (  bc_dirichlet_r )  s(:,:,nxr+1) = s(:,:,nxr)
+          IF (  bc_dirichlet_s )  s(:,nys-1,:) = s(:,nys,:)
+          IF (  bc_dirichlet_n )  s(:,nyn+1,:) = s(:,nyn,:)
+       ENDIF
+
+       CALL exchange_horiz( u, nbgp )
+       CALL exchange_horiz( v, nbgp )
+       CALL exchange_horiz( w, nbgp )
+       IF ( .NOT. neutral )  CALL exchange_horiz( pt, nbgp )
+       IF ( humidity      )  CALL exchange_horiz( q,  nbgp )
+!
+!--    In case of Rayleigh damping, where the profiles u_init, v_init
+!--    q_init and pt_init are still used, update these profiles from the
+!--    averaged boundary data. 
+!--    But first, average these data.
+#if defined( __parallel )
+       CALL MPI_ALLREDUCE( u_ref_l, u_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,     &
+                           comm2d, ierr )
+       CALL MPI_ALLREDUCE( v_ref_l, v_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,     &
+                           comm2d, ierr )
+       IF ( humidity )  THEN
+          CALL MPI_ALLREDUCE( q_ref_l, q_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,  &
+                              comm2d, ierr )
+       ENDIF
+       IF ( .NOT. neutral )  THEN
+          CALL MPI_ALLREDUCE( pt_ref_l, pt_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,&
+                              comm2d, ierr )
+       ENDIF
+#else
+    u_ref  = u_ref_l
+    v_ref  = v_ref_l
+    IF ( humidity )       q_ref  = q_ref_l
+    IF ( .NOT. neutral )  pt_ref = pt_ref_l
+#endif
+!
+!-- Average data. Note, reference profiles up to nzt are derived from lateral
+!-- boundaries, at the model top it is derived from the top boundary. Thus,
+!-- number of input data is different from nzb:nzt compared to nzt+1. 
+!-- Derived from lateral boundaries. 
+    u_ref(nzb:nzt) = u_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx     ),      &
+                                            KIND = wp )  
+    v_ref(nzb:nzt) = v_ref(nzb:nzt) / REAL( 2.0_wp * ( ny   + nx + 1   ),      &
+                                            KIND = wp )
+    IF ( humidity )                                                            &
+       q_ref(nzb:nzt) = q_ref(nzb:nzt)   / REAL( 2.0_wp * ( ny + 1 + nx + 1 ), &
+                                                 KIND = wp )
+    IF ( .NOT. neutral )                                                       &
+       pt_ref(nzb:nzt) = pt_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx + 1 ), &
+                                           KIND = wp )
+!
+!-- Derived from top boundary.    
+    u_ref(nzt+1) = u_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx     ), KIND = wp ) 
+    v_ref(nzt+1) = v_ref(nzt+1) / REAL( ( ny     ) * ( nx + 1 ), KIND = wp )
+    IF ( humidity )                                                            &
+       q_ref(nzt+1) = q_ref(nzt+1)   / REAL( ( ny + 1 ) * ( nx + 1 ), KIND = wp )
+    IF ( .NOT. neutral )                                                       &
+       pt_ref(nzt+1) = pt_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ), KIND = wp )
+!
+!-- Write onto init profiles, which are used for damping
+    u_init = u_ref
+    v_init = v_ref
+    IF ( humidity      )  q_init  = q_ref
+    IF ( .NOT. neutral )  pt_init = pt_ref
+!
+!-- Set bottom boundary condition
+    IF ( humidity      )  q_init(nzb)  = q_init(nzb+1)
+    IF ( .NOT. neutral )  pt_init(nzb) = pt_init(nzb+1)
+!
+!-- Further, adjust Rayleigh damping height in case of time-changing conditions.
+!-- Therefore, calculate boundary-layer depth first. 
+    CALL calc_zi
+    CALL adjust_sponge_layer 
+    
+!
+!-- Update geostrophic wind components from dynamic input file. 
+!-- Note, at the moment this is not available. As a workaround, use reference
+!-- profiles derived from the boundary values. 
+!-- This will be revised as soon as geostrophic wind components are provided by
+!-- dynamic input file.
+    ug = u_init
+    vg = v_init
+    
+    CALL  cpu_log( log_point(58), 'offline nesting', 'stop' )
+
+    END SUBROUTINE nesting_offl_bc
+
+!------------------------------------------------------------------------------!
+! Description:
+!------------------------------------------------------------------------------!
+!> Calculates the boundary-layer depth from the boundary data, according to 
+!> bulk-Richardson criterion. 
+!------------------------------------------------------------------------------!
+    SUBROUTINE calc_zi
+       
+       USE basic_constants_and_equations_mod,                                  &
+           ONLY:  g
+       
+       USE kinds
+       
+       USE surface_mod,                                                        &
+           ONLY:  get_topography_top_index, get_topography_top_index_ji
+
+       IMPLICIT NONE 
+
+       INTEGER(iwp) :: i            !< loop index in x-direction
+       INTEGER(iwp) :: j            !< loop index in y-direction
+       INTEGER(iwp) :: k            !< loop index in z-direction
+       INTEGER(iwp) :: k_surface    !< topography top index in z-direction
+    
+       REAL(wp) ::  ri_bulk                 !< bulk Richardson number
+       REAL(wp) ::  ri_bulk_crit = 0.25_wp  !< critical bulk Richardson number
+       REAL(wp) ::  topo_max                !< maximum topography level in model domain
+       REAL(wp) ::  topo_max_l              !< maximum topography level in subdomain
+       REAL(wp) ::  u_comp                  !< u-component
+       REAL(wp) ::  v_comp                  !< v-component
+       REAL(wp) ::  vpt_surface             !< near-surface virtual potential temperature
+       REAL(wp) ::  zi_l                    !< mean boundary-layer depth on subdomain   
+       REAL(wp) ::  zi_local                !< local boundary-layer depth  
+    
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  vpt_col !< vertical profile of virtual potential temperature at (j,i)-grid point
+
+       
+!
+!--    Calculate mean boundary-layer height from boundary data.
+!--    Start with the left and right boundaries. 
+       zi_l      = 0.0_wp
+       IF ( bc_dirichlet_l  .OR.  bc_dirichlet_r )  THEN
+!
+!--       Determine index along x. Please note, index indicates boundary
+!--       grid point for scalars. 
+          i = MERGE( -1, nxr + 1, bc_dirichlet_l )
+    
+          DO  j = nys, nyn
+!
+!--          Determine topography top index at current (j,i) index
+             k_surface = get_topography_top_index_ji( j, i, 's' )
+!
+!--          Pre-compute surface virtual temperature. Therefore, use 2nd 
+!--          prognostic level according to Heinze et al. (2017).
+             IF ( humidity )  THEN
+                vpt_surface = pt(k_surface+2,j,i) *                            &
+                            ( 1.0_wp + 0.61_wp * q(k_surface+2,j,i) )
+                vpt_col     = pt(:,j,i) * ( 1.0_wp + 0.61_wp * q(:,j,i) )
+             ELSE
+                vpt_surface = pt(k_surface+2,j,i)
+                vpt_col     = pt(:,j,i)
+             ENDIF
+!
+!--          Calculate local boundary layer height from bulk Richardson number,
+!--          i.e. the height where the bulk Richardson number exceeds its
+!--          critical value of 0.25 (according to Heinze et al., 2017).
+!--          Note, no interpolation of u- and v-component is made, as both 
+!--          are mainly mean inflow profiles with very small spatial variation. 
+             zi_local = 0.0_wp
+             DO  k = k_surface+1, nzt
+                u_comp = MERGE( u(k,j,i+1), u(k,j,i), bc_dirichlet_l )
+                v_comp = v(k,j,i)
+                ri_bulk = zu(k) * g / vpt_surface *                            &
+                          ( vpt_col(k) - vpt_surface ) /                       &
+                          ( u_comp * u_comp + v_comp * v_comp ) 
+                       
+                IF ( zi_local == 0.0_wp  .AND.  ri_bulk > ri_bulk_crit )       &
+                   zi_local = zu(k)           
+             ENDDO
+!
+!--          Assure that the minimum local boundary-layer depth is at least at 
+!--          the second vertical grid level.
+             zi_l = zi_l + MAX( zi_local, zu(k_surface+2) )   
+             
+          ENDDO
+       
+       ENDIF
+!
+!--    Do the same at the north and south boundaries.
+       IF ( bc_dirichlet_s  .OR.  bc_dirichlet_n )  THEN
+       
+          j = MERGE( -1, nyn + 1, bc_dirichlet_s )
+       
+          DO  i = nxl, nxr
+             k_surface = get_topography_top_index_ji( j, i, 's' )
+ 
+             IF ( humidity )  THEN
+                vpt_surface = pt(k_surface+2,j,i) *                            &
+                            ( 1.0_wp + 0.61_wp * q(k_surface+2,j,i) )
+                vpt_col     = pt(:,j,i) * ( 1.0_wp + 0.61_wp * q(:,j,i) )
+             ELSE
+                vpt_surface = pt(k_surface+2,j,i)
+                vpt_col  = pt(:,j,i)
+             ENDIF
+          
+             zi_local = 0.0_wp
+             DO  k = k_surface+1, nzt                
+                u_comp = u(k,j,i)
+                v_comp = MERGE( v(k,j+1,i), v(k,j,i), bc_dirichlet_s )
+                ri_bulk = zu(k) * g / vpt_surface *                            &
+                       ( vpt_col(k) - vpt_surface ) /                          &
+                       ( u_comp * u_comp + v_comp * v_comp ) 
+                       
+                IF ( zi_local == 0.0_wp  .AND.  ri_bulk > 0.25_wp )            &
+                   zi_local = zu(k)    
+             ENDDO
+             zi_l = zi_l + MAX( zi_local, zu(k_surface+2) )   
+          
+          ENDDO
+          
+       ENDIF
+    
+#if defined( __parallel )
+       CALL MPI_ALLREDUCE( zi_l, zi_ribulk, 1, MPI_REAL, MPI_SUM,              &
+                           comm2d, ierr )
+#else
+       zi_ribulk = zi_l
+#endif
+       zi_ribulk = zi_ribulk / REAL( 2 * nx + 2 * ny, KIND = wp )
+!
+!--    Finally, check if boundary layer depth is not below the any topography. 
+!--    zi_ribulk will be used to adjust rayleigh damping height, i.e. the 
+!--    lower level of the sponge layer, as well as to adjust the synthetic
+!--    turbulence generator accordingly. If Rayleigh damping would be applied
+!--    near buildings, etc., this would spoil the simulation results. 
+       topo_max_l = zw(MAXVAL( get_topography_top_index( 's' )))
+       
+#if defined( __parallel )
+       CALL MPI_ALLREDUCE( topo_max_l, topo_max, 1, MPI_REAL, MPI_MAX,            &
+                           comm2d, ierr )
+#else
+       topo_max     = topo_max_l
+#endif
+
+       zi_ribulk = MAX( zi_ribulk, topo_max )
+
+    END SUBROUTINE calc_zi
+    
+    
+!------------------------------------------------------------------------------!
+! Description:
+!------------------------------------------------------------------------------!
+!> Adjust the height where the rayleigh damping starts, i.e. the lower level
+!> of the sponge layer. 
+!------------------------------------------------------------------------------!
+    SUBROUTINE adjust_sponge_layer
+       
+       USE arrays_3d,                                                          &
+           ONLY:  rdf, rdf_sc, zu
+       
+       USE basic_constants_and_equations_mod,                                  &
+           ONLY:  pi
+           
+       USE control_parameters,                                                 &
+           ONLY:  rayleigh_damping_height, rayleigh_damping_factor
+       
+       USE kinds
+
+       IMPLICIT NONE 
+
+       INTEGER(iwp) :: k   !< loop index in z-direction
+    
+       REAL(wp) ::  rdh    !< updated Rayleigh damping height
+ 
+    
+       IF ( rayleigh_damping_height > 0.0_wp  .AND.                            &
+            rayleigh_damping_factor > 0.0_wp )  THEN
+!
+!--       Update Rayleigh-damping height and re-calculate height-depending 
+!--       damping coefficients. 
+!--       Assure that rayleigh damping starts well above the boundary layer.   
+          rdh = MIN( MAX( zi_ribulk * 1.3_wp, 10.0_wp * dz(1) ),               & 
+                     0.8_wp * zu(nzt), rayleigh_damping_height )
+!
+!--       Update Rayleigh damping factor
+          DO  k = nzb+1, nzt
+             IF ( zu(k) >= rdh )  THEN
+                rdf(k) = rayleigh_damping_factor *                             &
+                                          ( SIN( pi * 0.5_wp * ( zu(k) - rdh ) &
+                                        / ( zu(nzt) - rdh ) )                  &
+                                          )**2
+             ENDIF
+          ENDDO
+          rdf_sc = rdf
+       
+       ENDIF
+
+    END SUBROUTINE adjust_sponge_layer
+    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Performs consistency checks
+!------------------------------------------------------------------------------!
+    SUBROUTINE nesting_offl_check_parameters 
+
+       IMPLICIT NONE
+!
+!--    Perform checks
+
+
+    END SUBROUTINE nesting_offl_check_parameters
+    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Reads the parameter list nesting_offl_parameters
+!------------------------------------------------------------------------------!
+    SUBROUTINE nesting_offl_parin 
+
+       IMPLICIT NONE
+       
+       CHARACTER (LEN=80) ::  line   !< dummy string that contains the current line of the parameter file
+
+
+       NAMELIST /nesting_offl_parameters/   nesting_offline
+
+       line = ' '
+
+!
+!--    Try to find stg package
+       REWIND ( 11 )
+       line = ' '
+       DO WHILE ( INDEX( line, '&nesting_offl_parameters' ) == 0 )
+          READ ( 11, '(A)', END=20 )  line
+       ENDDO
+       BACKSPACE ( 11 )
+
+!
+!--    Read namelist
+       READ ( 11, nesting_offl_parameters, ERR = 10, END = 20 )
+
+       GOTO 20
+
+ 10    BACKSPACE( 11 )
+       READ( 11 , '(A)') line
+       CALL parin_fail_message( 'nesting_offl_parameters', line )
+
+ 20    CONTINUE
+
+
+    END SUBROUTINE nesting_offl_parin
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Writes information about offline nesting into HEADER file
+!------------------------------------------------------------------------------!
+    SUBROUTINE nesting_offl_header ( io )
+
+       IMPLICIT NONE
+
+       INTEGER(iwp), INTENT(IN) ::  io !< Unit of the output file
+
+       WRITE ( io, 1 )
+       IF ( nesting_offline )  THEN
+          WRITE ( io, 3 )
+       ELSE
+          WRITE ( io, 2 )
+       ENDIF
+
+1 FORMAT (//' Offline nesting in COSMO model:'/                                &
+              ' -------------------------------'/)
+2 FORMAT (' --> No offlince nesting is used (default) ')
+3 FORMAT (' --> Offlince nesting is used. Boundary data is read from dynamic input file ')
+
+    END SUBROUTINE nesting_offl_header 
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Allocate arrays used to read boundary data from NetCDF file and initialize
+!> boundary data. 
+!------------------------------------------------------------------------------!
+    SUBROUTINE nesting_offl_init
+    
+       USE netcdf_data_input_mod,                                              &
+           ONLY:  netcdf_data_input_offline_nesting  
+
+       IMPLICIT NONE
+
+
+!--    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( nest_offl%ug(0:1,nzb:nzt+1) )
+!           ALLOCATE( nest_offl%vg(0:1,nzb:nzt+1) )
+!        ENDIF
+!
+!--    Allocate arrays for reading boundary values. Arrays will incorporate 2 
+!--    time levels in order to interpolate in between. 
+       IF ( bc_dirichlet_l )  THEN
+          ALLOCATE( nest_offl%u_left(0:1,nzb+1:nzt,nys:nyn)  )
+          ALLOCATE( nest_offl%v_left(0:1,nzb+1:nzt,nysv:nyn) )
+          ALLOCATE( nest_offl%w_left(0:1,nzb+1:nzt-1,nys:nyn) )
+          IF ( humidity )       ALLOCATE( nest_offl%q_left(0:1,nzb+1:nzt,nys:nyn)  )
+          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_left(0:1,nzb+1:nzt,nys:nyn) )
+       ENDIF
+       IF ( bc_dirichlet_r )  THEN
+          ALLOCATE( nest_offl%u_right(0:1,nzb+1:nzt,nys:nyn)  )
+          ALLOCATE( nest_offl%v_right(0:1,nzb+1:nzt,nysv:nyn) )
+          ALLOCATE( nest_offl%w_right(0:1,nzb+1:nzt-1,nys:nyn) )
+          IF ( humidity )       ALLOCATE( nest_offl%q_right(0:1,nzb+1:nzt,nys:nyn)  )
+          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_right(0:1,nzb+1:nzt,nys:nyn) )
+       ENDIF
+       IF ( bc_dirichlet_n )  THEN
+          ALLOCATE( nest_offl%u_north(0:1,nzb+1:nzt,nxlu:nxr) )
+          ALLOCATE( nest_offl%v_north(0:1,nzb+1:nzt,nxl:nxr)  )
+          ALLOCATE( nest_offl%w_north(0:1,nzb+1:nzt-1,nxl:nxr) )
+          IF ( humidity )       ALLOCATE( nest_offl%q_north(0:1,nzb+1:nzt,nxl:nxr)  )
+          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_north(0:1,nzb+1:nzt,nxl:nxr) )
+       ENDIF
+       IF ( bc_dirichlet_s )  THEN
+          ALLOCATE( nest_offl%u_south(0:1,nzb+1:nzt,nxlu:nxr) )
+          ALLOCATE( nest_offl%v_south(0:1,nzb+1:nzt,nxl:nxr)  )
+          ALLOCATE( nest_offl%w_south(0:1,nzb+1:nzt-1,nxl:nxr)    )
+          IF ( humidity )       ALLOCATE( nest_offl%q_south(0:1,nzb+1:nzt,nxl:nxr)  )
+          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_south(0:1,nzb+1:nzt,nxl:nxr) )
+       ENDIF
+       
+       ALLOCATE( nest_offl%u_top(0:1,nys:nyn,nxlu:nxr) )
+       ALLOCATE( nest_offl%v_top(0:1,nysv:nyn,nxl:nxr) )
+       ALLOCATE( nest_offl%w_top(0:1,nys:nyn,nxl:nxr)  )
+       IF ( humidity )       ALLOCATE( nest_offl%q_top(0:1,nys:nyn,nxl:nxr)  )
+       IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_top(0:1,nys:nyn,nxl:nxr) )
+
+!
+!--    Read COSMO data at lateral and top boundaries
+       CALL netcdf_data_input_offline_nesting
+!
+!--    Initialize boundary data 
+       IF ( bc_dirichlet_l )  THEN
+          u(nzb+1:nzt,nys:nyn,0)    = nest_offl%u_left(0,nzb+1:nzt,nys:nyn)
+          v(nzb+1:nzt,nysv:nyn,-1)  = nest_offl%v_left(0,nzb+1:nzt,nysv:nyn) 
+          w(nzb+1:nzt-1,nys:nyn,-1) = nest_offl%w_left(0,nzb+1:nzt-1,nys:nyn)
+          IF ( .NOT. neutral )  pt(nzb+1:nzt,nys:nyn,-1) =                     &
+                                   nest_offl%pt_left(0,nzb+1:nzt,nys:nyn)
+          IF ( humidity      )  q(nzb+1:nzt,nys:nyn,-1)  =                     &
+                                   nest_offl%q_left(0,nzb+1:nzt,nys:nyn) 
+       ENDIF
+       IF ( bc_dirichlet_r )  THEN
+          u(nzb+1:nzt,nys:nyn,nxr+1)   = nest_offl%u_right(0,nzb+1:nzt,nys:nyn) 
+          v(nzb+1:nzt,nysv:nyn,nxr+1)  = nest_offl%v_right(0,nzb+1:nzt,nysv:nyn)
+          w(nzb+1:nzt-1,nys:nyn,nxr+1) = nest_offl%w_right(0,nzb+1:nzt-1,nys:nyn)
+          IF ( .NOT. neutral )  pt(nzb+1:nzt,nys:nyn,nxr+1) =                  &
+                                   nest_offl%pt_right(0,nzb+1:nzt,nys:nyn)
+          IF ( humidity      )  q(nzb+1:nzt,nys:nyn,nxr+1)  =                  &
+                                   nest_offl%q_right(0,nzb+1:nzt,nys:nyn)
+       ENDIF
+       IF ( bc_dirichlet_s )  THEN
+          u(nzb+1:nzt,-1,nxlu:nxr)  = nest_offl%u_south(0,nzb+1:nzt,nxlu:nxr)
+          v(nzb+1:nzt,0,nxl:nxr)    = nest_offl%v_south(0,nzb+1:nzt,nxl:nxr) 
+          w(nzb+1:nzt-1,-1,nxl:nxr) = nest_offl%w_south(0,nzb+1:nzt-1,nxl:nxr) 
+          IF ( .NOT. neutral )  pt(nzb+1:nzt,-1,nxl:nxr) =                     &
+                                   nest_offl%pt_south(0,nzb+1:nzt,nxl:nxr)
+          IF ( humidity      )  q(nzb+1:nzt,-1,nxl:nxr)  =                     &
+                                   nest_offl%q_south(0,nzb+1:nzt,nxl:nxr) 
+       ENDIF
+       IF ( bc_dirichlet_n )  THEN
+          u(nzb+1:nzt,nyn+1,nxlu:nxr)  = nest_offl%u_north(0,nzb+1:nzt,nxlu:nxr)
+          v(nzb+1:nzt,nyn+1,nxl:nxr)   = nest_offl%v_north(0,nzb+1:nzt,nxl:nxr) 
+          w(nzb+1:nzt-1,nyn+1,nxl:nxr) = nest_offl%w_north(0,nzb+1:nzt-1,nxl:nxr)
+          IF ( .NOT. neutral )  pt(nzb+1:nzt,nyn+1,nxl:nxr) =                  &
+                                   nest_offl%pt_north(0,nzb+1:nzt,nxl:nxr) 
+          IF ( humidity      )  q(nzb+1:nzt,nyn+1,nxl:nxr)  =                  &
+                                   nest_offl%q_north(0,nzb+1:nzt,nxl:nxr)
+       ENDIF
+!
+!--    Initial calculation of the boundary layer depth from the prescribed 
+!--    boundary data. This is requiered for initialize the synthetic turbulence
+!--    generator correctly. 
+       CALL calc_zi
+       
+!
+!--    After boundary data is initialized, mask topography at the 
+!--    boundaries for the velocity components.
+       u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
+       v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )
+       w = MERGE( w, 0.0_wp, BTEST( wall_flags_0, 3 ) )
+       
+       CALL calc_zi
+       
+!
+!--    After boundary data is initialized, ensure mass conservation
+       CALL nesting_offl_mass_conservation
+
+    END SUBROUTINE nesting_offl_init
+    
+!------------------------------------------------------------------------------!
+! Description:
+!------------------------------------------------------------------------------!
+!> Interpolation function, used to interpolate boundary data in time. 
+!------------------------------------------------------------------------------!
+    FUNCTION interpolate_in_time( var_t1, var_t2, fac  ) 
+       
+       USE kinds
+
+       IMPLICIT NONE 
+
+       REAL(wp)            :: interpolate_in_time !< time-interpolated boundary value
+       REAL(wp)            :: var_t1              !< boundary value at t1
+       REAL(wp)            :: var_t2              !< boundary value at t2
+       REAL(wp)            :: fac                 !< interpolation factor
+
+       interpolate_in_time = ( 1.0_wp - fac ) * var_t1 + fac * var_t2      
+
+    END FUNCTION interpolate_in_time
+
+
+
+ END MODULE nesting_offl_mod
Index: /palm/trunk/SOURCE/netcdf_data_input_mod.f90
===================================================================
--- /palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 3346)
+++ /palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 3347)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Subroutine renamed
+! 
+! 3257 2018-09-17 17:11:46Z suehring
 ! (from branch resler)
 ! Formatting
@@ -678,7 +681,7 @@
     END INTERFACE netcdf_data_input_init_lsm
 
-    INTERFACE netcdf_data_input_lsf
-       MODULE PROCEDURE netcdf_data_input_lsf
-    END INTERFACE netcdf_data_input_lsf
+    INTERFACE netcdf_data_input_offline_nesting
+       MODULE PROCEDURE netcdf_data_input_offline_nesting
+    END INTERFACE netcdf_data_input_offline_nesting
 
     INTERFACE netcdf_data_input_surface_data
@@ -735,5 +738,5 @@
            netcdf_data_input_init, netcdf_data_input_init_lsm,                 &
            netcdf_data_input_init_3d,                                          &
-           netcdf_data_input_interpolate, netcdf_data_input_lsf,               &
+           netcdf_data_input_interpolate, netcdf_data_input_offline_nesting,   &
            netcdf_data_input_surface_data, netcdf_data_input_topo
 
@@ -3078,5 +3081,5 @@
 !> (COSMO) by Inifor.
 !------------------------------------------------------------------------------!
-    SUBROUTINE netcdf_data_input_lsf
+    SUBROUTINE netcdf_data_input_offline_nesting
 
        USE control_parameters,                                                 &
@@ -3157,5 +3160,5 @@
 !
 !--    Obtain time index for current input starting at 0.
-!--    @todo: At the moment time, in INIFOR and simulated time correspond
+!--    @todo: At the moment INIFOR and simulated time correspond
 !--           to each other. If required, adjust to daytime.
        nest_offl%tind = MINLOC( ABS( nest_offl%time -                          &
@@ -3348,5 +3351,5 @@
        CALL cpu_log( log_point_s(86), 'NetCDF input forcing', 'stop' )
 
-    END SUBROUTINE netcdf_data_input_lsf
+    END SUBROUTINE netcdf_data_input_offline_nesting
 
 
Index: /palm/trunk/SOURCE/parin.f90
===================================================================
--- /palm/trunk/SOURCE/parin.f90	(revision 3346)
+++ /palm/trunk/SOURCE/parin.f90	(revision 3347)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! - offline nesting separated from large-scale forcing module
+! - top boundary condition for pressure in offline nesting changed
+! 
+! 3343 2018-10-15 10:38:52Z suehring
 ! Introduced reading of date_init to inipar.(Russo)
 ! Add extra profiles for chemistry (basit)
@@ -477,4 +481,7 @@
         ONLY:  mas_parin
 
+    USE nesting_offl_mod,                                                      &
+        ONLY:  nesting_offl_parin
+        
     USE netcdf_interface,                                                      &
         ONLY:  netcdf_data_format, netcdf_deflate, netcdf_precision
@@ -570,5 +577,5 @@
              loop_optimization, lsf_exception, masking_method, mg_cycles,      &
              mg_switch_to_pe0_level, mixing_length_1d, momentum_advec,         &
-             most_method, nesting_offline,                  &
+             most_method,                                                      &
              netcdf_precision, neutral, ngsrb,                                 &
              nsor, nsor_ini, nudging, nx, ny, nz, ocean_mode, omega,           &
@@ -643,5 +650,5 @@
              loop_optimization, lsf_exception, masking_method, mg_cycles,      &
              mg_switch_to_pe0_level, mixing_length_1d, momentum_advec,         &
-             most_method, nesting_offline,                                     &
+             most_method,                                                      &
              netcdf_precision, neutral, ngsrb,                                 &
              nsor, nsor_ini, nudging, nx, ny, nz, ocean_mode, omega,           &
@@ -878,4 +885,5 @@
           CALL gust_parin
           CALL mas_parin
+          CALL nesting_offl_parin
           CALL ocean_parin
           CALL pcm_parin
@@ -979,5 +987,12 @@
            !  bc_s_t   = 'nesting_offline'  ! scalar boundary condition is not clear
            !  bc_cs_t  = 'nesting_offline'  ! same for chemical species
-             bc_p_t   = 'neumann'
+!
+!--          For the pressure set Dirichlet conditions, in contrast to the 
+!--          self nesting. This gives less oscilations within the 
+!--          free atmosphere since the pressure solver has more degrees of 
+!--          freedom. In constrast to the self nesting, this might be 
+!--          justified since the top boundary is far away from the domain
+!--          of interest. 
+             bc_p_t   = 'dirichlet' !'neumann'
           ENDIF
 
Index: /palm/trunk/SOURCE/pres.f90
===================================================================
--- /palm/trunk/SOURCE/pres.f90	(revision 3346)
+++ /palm/trunk/SOURCE/pres.f90	(revision 3347)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Bugfixes in offline nesting.
+! Add comment. 
+! 
+! 3341 2018-10-15 10:31:27Z suehring
 ! Rename variables for boundary flags and nesting
 ! 
@@ -169,5 +173,6 @@
                dt_3d, gathered_size, ibc_p_b, ibc_p_t,                         &
                intermediate_timestep_count, intermediate_timestep_count_max,   &
-               mg_switch_to_pe0_level, psolver, subdomain_size,                &
+               mg_switch_to_pe0_level, nesting_offline,                        &
+               psolver, subdomain_size,                                        &
                topography, volume_flow, volume_flow_area, volume_flow_initial
 
@@ -386,5 +391,5 @@
     ENDIF
 
-    IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1  .AND.                               &
+    IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1  .AND.  .NOT. nesting_offline  .AND. &
          .NOT. child_domain_nvn  .AND. intermediate_timestep_count /= 0 )       &
     THEN
@@ -719,5 +724,11 @@
        volume_flow_l(2) = 0.0_wp
     ENDIF
-
+!
+!-- Add pressure gradients to the velocity components. Note, the loops are
+!-- running over the entire model domain, even though, in case of non-cyclic 
+!-- boundaries u- and v-component are not prognostic at i=0 and j=0, 
+!-- respectiveley. However, in case of Dirichlet boundary conditions for the 
+!-- velocities, zero-gradient conditions for the pressure are set, so that
+!-- no modification is imposed at the boundaries. 
     !$OMP PARALLEL PRIVATE (i,j,k)
     !$OMP DO
Index: /palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90
===================================================================
--- /palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90	(revision 3346)
+++ /palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90	(revision 3347)
@@ -25,4 +25,10 @@
 ! -----------------
 ! $Id$
+! - Revise structure of init routine
+! - introduce new parameters to skip STG for some timesteps
+! - introduce time-dependent parametrization of Reynolds-stress tensor
+! - Bugfix in allocation of mean_inflow_profiles
+! 
+! 3341 2018-10-15 10:31:27Z suehring
 ! Introduce module parameter for number of inflow profiles
 ! 
@@ -122,5 +128,6 @@
 ! Authors:
 ! --------
-! @author Tobias Gronemeier, Atsushi Inagaki, Micha Gryschka, Christoph Knigge
+! @author Tobias Gronemeier, Matthias Suehring, Atsushi Inagaki, Micha Gryschka, Christoph Knigge
+! 
 !
 !
@@ -133,5 +140,6 @@
 !> a time scale for each velocity component. The profiles of length and time
 !> scales, mean u, v, w, e and pt, and all components of the Reynolds stress
-!> tensor are read from file STG_PROFILES.
+!> tensor can be either read from file STG_PROFILES, or will be parametrized
+!> within the boundary layer. 
 !>
 !> @todo test restart
@@ -149,8 +157,8 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  mean_inflow_profiles, u, v, w
+        ONLY:  mean_inflow_profiles, q, pt, u, v, w, zu, zw
 
     USE basic_constants_and_equations_mod,                                     &
-        ONLY:  pi
+        ONLY:  g, kappa, pi
 
     USE control_parameters,                                                    &
@@ -159,5 +167,5 @@
 
     USE cpulog,                                                                &
-        ONLY:  cpu_log, log_point, log_point_s
+        ONLY:  cpu_log, log_point
 
     USE indices,                                                               &
@@ -170,4 +178,7 @@
 #endif
 
+    USE nesting_offl_mod,                                                      &
+        ONLY:  zi_ribulk
+
     USE pegrid,                                                                &
         ONLY:  comm1dx, comm1dy, comm2d, ierr, myidx, myidy, pdims
@@ -184,6 +195,7 @@
 
 
-    LOGICAL :: velocity_seed_initialized = .FALSE.  !< true after first call of stg_main
-    LOGICAL :: use_syn_turb_gen = .FALSE.           !< switch to use synthetic turbulence generator
+    LOGICAL ::  velocity_seed_initialized = .FALSE.     !< true after first call of stg_main
+    LOGICAL ::  parametrize_inflow_turbulence = .FALSE. !< flag indicating that inflow turbulence is either read from file (.FALSE.) or if it parametrized
+    LOGICAL ::  use_syn_turb_gen = .FALSE.              !< switch to use synthetic turbulence generator
 
     INTEGER(iwp) ::  id_stg_left        !< left lateral boundary core id in case of turbulence generator
@@ -202,57 +214,79 @@
     INTEGER(iwp) ::  nzt_y_stg          !< upper bound of z coordinate (required for transposing z on PEs along y)
 
-    REAL(wp) :: mc_factor    !< mass flux correction factor
-
-    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
-
-    REAL(wp), DIMENSION(:), ALLOCATABLE :: a11             !< coefficient for Lund rotation
-    REAL(wp), DIMENSION(:), ALLOCATABLE :: a21             !< coefficient for Lund rotation
-    REAL(wp), DIMENSION(:), ALLOCATABLE :: a22             !< coefficient for Lund rotation
-    REAL(wp), DIMENSION(:), ALLOCATABLE :: a31             !< coefficient for Lund rotation
-    REAL(wp), DIMENSION(:), ALLOCATABLE :: a32             !< coefficient for Lund rotation
-    REAL(wp), DIMENSION(:), ALLOCATABLE :: a33             !< coefficient for Lund rotation
-    REAL(wp), DIMENSION(:), ALLOCATABLE :: tu              !< Lagrangian time scale of u
-    REAL(wp), DIMENSION(:), ALLOCATABLE :: tv              !< Lagrangian time scale of v
-    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_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
-
+    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
+
+    REAL(wp) ::  cosmo_ref = 10.0_wp      !< height of first vertical grid level in mesoscale model, used for calculation of scaling parameters
+    REAL(wp) ::  dt_stg_adjust = 300.0_wp !< time interval for adjusting turbulence statistics
+    REAL(wp) ::  dt_stg_call = 5.0_wp     !< time interval for calling synthetic turbulence generator
+    REAL(wp) ::  mc_factor                !< mass flux correction factor
+    REAL(wp) ::  scale_l                  !< scaling parameter used for turbulence parametrization - Obukhov length
+    REAL(wp) ::  scale_us                 !< scaling parameter used for turbulence parametrization - friction velocity
+    REAL(wp) ::  scale_wm                 !< scaling parameter used for turbulence parametrization - momentum scale  
+    REAL(wp) ::  time_stg_adjust = 0.0_wp !< time counter for adjusting turbulence information   
+    REAL(wp) ::  time_stg_call = 0.0_wp   !< time counter for calling generator   
+    
+    
+    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r11              !< Reynolds parameter
+    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r21              !< Reynolds parameter
+    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r22              !< Reynolds parameter
+    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r31              !< Reynolds parameter
+    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r32              !< Reynolds parameter
+    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r33              !< Reynolds parameter
+    
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a11             !< coefficient for Lund rotation
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a21             !< coefficient for Lund rotation
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a22             !< coefficient for Lund rotation
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a31             !< coefficient for Lund rotation
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a32             !< coefficient for Lund rotation
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a33             !< coefficient for Lund rotation
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  tu              !< Lagrangian time scale of u
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  tv              !< Lagrangian time scale of v
+    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_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
+    
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist_xz     !< imposed disturbances at north/south boundary
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist_yz     !< imposed disturbances at north/south boundary
 
 !
 !-- PALM interfaces:
+!-- Adjust time and lenght scales, Reynolds stress, and filter functions
+    INTERFACE stg_adjust
+       MODULE PROCEDURE stg_adjust
+    END INTERFACE stg_adjust
+!
 !-- Input parameter checks to be done in check_parameters
     INTERFACE stg_check_parameters
@@ -319,11 +353,12 @@
 !
 !-- Public interfaces
-    PUBLIC  stg_check_parameters, stg_header, stg_init, stg_main, stg_parin,   &
-            stg_wrd_global, stg_rrd_global
+    PUBLIC  stg_adjust, stg_check_parameters, stg_header, stg_init, stg_main,  &
+            stg_parin, stg_rrd_global, stg_wrd_global
 
 !
 !-- Public variables
-    PUBLIC  id_stg_left, id_stg_north, id_stg_right, id_stg_south,             &
-            use_syn_turb_gen
+    PUBLIC  dt_stg_call, dt_stg_adjust, id_stg_left, id_stg_north,             &
+            id_stg_right, id_stg_south, parametrize_inflow_turbulence,         &
+            time_stg_adjust, time_stg_call, use_syn_turb_gen
 
 
@@ -376,5 +411,6 @@
              message_string = 'Using synthetic turbulence generator ' //       &
                               'requires %initializing_actions = '         //   &
-                              '"set_constant_profiles" or "read_restart_data"'
+                              '"set_constant_profiles" or "read_restart_data"' //&
+                              ', if not offline nesting is applied.'
              CALL message( 'stg_check_parameters', 'PA0015', 1, 2, 0, 6, 0 )
           ENDIF
@@ -382,10 +418,12 @@
           IF ( bc_lr /= 'dirichlet/radiation' )  THEN
              message_string = 'Using synthetic turbulence generator ' //       &
-                              'requires &bc_lr = "dirichlet/radiation"'
+                              'requires &bc_lr = "dirichlet/radiation", ' //   &
+                              'if not offline nesting is applied.'
              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"'
+                              'requires &bc_ns = "cyclic", ' //               &
+                              'if not offline nesting is applied.'
              CALL message( 'stg_check_parameters', 'PA0037', 1, 2, 0, 6, 0 )
           ENDIF
@@ -425,9 +463,18 @@
        WRITE( io, 3 )
     ENDIF
+    
+    IF ( parametrize_inflow_turbulence )  THEN
+       WRITE( io, 4 ) dt_stg_adjust
+    ELSE
+       WRITE( io, 5 )
+    ENDIF
 
 1   FORMAT (//' Synthetic turbulence generator information:'/                  &
               ' ------------------------------------------'/)
-2   FORMAT ('    synthetic turbulence generator switched on')
-3   FORMAT ('    synthetic turbulence generator switched off')
+2   FORMAT ('    synthetic turbulence generator is switched on')
+3   FORMAT ('    synthetic turbulence generator is switched off')
+4   FORMAT ('    imposed turbulence statistics are parametrized and ' \\       &
+            'ajdusted to boundary-layer development each ', F8.2, ' s' )
+5   FORMAT ('    imposed turbulence is read from file' )
 
  END SUBROUTINE stg_header
@@ -443,5 +490,5 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  dzw, ddzw, u_init, v_init, zu, zw
+        ONLY:  dzw, ddzw, u_init, v_init, zu
 
     USE control_parameters,                                                    &
@@ -460,5 +507,5 @@
     IMPLICIT NONE
 
-    LOGICAL ::  file_stg_exist = .FALSE. !< flag indication whether parameter file for Reynolds stress and length scales exist
+    LOGICAL ::  file_stg_exist = .FALSE. !< flag indicating whether parameter file for Reynolds stress and length scales exist
 
 #if defined( __parallel )
@@ -484,5 +531,4 @@
     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
@@ -494,16 +540,4 @@
     REAL(wp) :: zz      !< height
 
-    REAL(wp) :: length_scale_surface !< typical length scale
-    REAL(wp) :: r_ii_0               !< correction factor
-    REAL(wp) :: time_scale           !< typical time scale
-    REAL(wp) :: length_scale_z       !< typical length scale
-
-    REAL(wp),DIMENSION(nzb:nzt+1) :: r11  !< Reynolds parameter
-    REAL(wp),DIMENSION(nzb:nzt+1) :: r21  !< Reynolds parameter
-    REAL(wp),DIMENSION(nzb:nzt+1) :: r22  !< Reynolds parameter
-    REAL(wp),DIMENSION(nzb:nzt+1) :: r31  !< Reynolds parameter
-    REAL(wp),DIMENSION(nzb:nzt+1) :: r32  !< Reynolds parameter
-    REAL(wp),DIMENSION(nzb:nzt+1) :: r33  !< Reynolds parameter
-
 
 #if defined( __parallel )
@@ -512,14 +546,4 @@
 
     CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
-
-    IF ( .NOT. ALLOCATED( mean_inflow_profiles ) )                             &
-       ALLOCATE( mean_inflow_profiles(nzb:nzt+1,1:num_mean_inflow_profiles) )
-
-    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),                 &
-               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 )
@@ -615,5 +639,26 @@
     ENDDO
     CALL RANDOM_SEED(put=seed)
-
+!
+!-- Allocate required arrays
+!-- mean_inflow profiles must not be allocated in offline nesting
+    IF ( .NOT. nesting_offline  .AND.  .NOT.  child_domain )  THEN 
+       IF ( .NOT. ALLOCATED( mean_inflow_profiles ) )                          &
+          ALLOCATE( mean_inflow_profiles(nzb:nzt+1,1:num_mean_inflow_profiles) )
+    ENDIF
+
+    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),                 &
+               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),                 &
+               r11(nzb:nzt+1), r21(nzb:nzt+1), r22(nzb:nzt+1),                 &
+               r31(nzb:nzt+1), r32(nzb:nzt+1), r33(nzb:nzt+1),                 &
+               tu(nzb:nzt+1),  tv(nzb:nzt+1),  tw(nzb:nzt+1)   )
+               
+    ALLOCATE ( dist_xz(nzb:nzt+1,nxlg:nxrg,3) )
+    ALLOCATE ( dist_yz(nzb:nzt+1,nysg:nyng,3) )
+    dist_xz = 0.0_wp
+    dist_yz = 0.0_wp
+!
 !-- Read inflow profile
 !-- Height levels of profiles in input profile is as follows:
@@ -669,80 +714,34 @@
        
        CLOSE( 90 )
-
+!
+!--    Calculate coefficient matrix from Reynolds stress tensor  
+!--    (Lund rotation)
+       CALL calc_coeff_matrix
+!
+!-- No information about turbulence and its length scales are available.
+!-- Instead, parametrize turbulence which is imposed at the boundaries. 
+!-- Set flag which indicates that turbulence is parametrized, which is done 
+!-- later when energy-balance models are already initialized. This is 
+!-- because the STG needs information about surface properties such as 
+!-- roughness to build 'realistic' turbulence profiles. 
     ELSE
 !
-!--    Set-up default turbulent length scales. From the numerical point of
-!--    view the imposed perturbations should not be immediately dissipated
-!--    by the numerics. The numerical dissipation, however, acts on scales
-!--    up to 8 x the grid spacing. For this reason, set the turbulence 
-!--    length scale to 8 time the grid spacing. Further, above the boundary
-!--    layer height, set turbulence lenght scales to zero (equivalent to not
-!--    imposing any perturbations) in order to save computational costs. 
-!--    Typical length (time) scales of 100 m (s) should be a good 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 = 8.0_wp * MIN( dx, dy, dzw(k) )
-          
-!           IF ( zu(k) > zi_richardson )  
-
-          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
-
+!--    Set flag indicating that turbulence is parametrized
+       parametrize_inflow_turbulence = .TRUE.
+!
+!--    Determine boundary-layer depth, which is used to initialize lenght 
+!--    scales
+       CALL calc_scaling_variables
+!
+!--    Initialize lenght and time scales, which in turn are used
+!--    to initialize the filter functions.
+       CALL calc_length_and_time_scale
+            
     ENDIF
 
 !
-!-- Assign initial profiles
+!-- Assign initial profiles. Note, this is only required if turbulent
+!-- inflow from the left is desired, not in case of any of the 
+!-- nesting (offline or self nesting) approaches.
     IF ( .NOT. nesting_offline  .AND.  .NOT.  child_domain )  THEN 
        u_init = mean_inflow_profiles(:,1)
@@ -751,56 +750,5 @@
        e_init = MAXVAL( mean_inflow_profiles(:,5) )
     ENDIF
-!
-!-- Calculate coefficient matrix from Reynolds stress tensor (Lund rotation)
-    DO  k = nzb, nzt+1
-       IF ( r11(k) > 0.0_wp )  THEN
-          a11(k) = SQRT( r11(k) )
-          a21(k) = r21(k) / a11(k)
-       ELSE
-          a11(k) = 0.0_wp
-          a21(k) = 0.0_wp
-       ENDIF
-
-       a22(k) = r22(k) - a21(k)**2
-       IF ( a22(k) > 0.0_wp )  THEN
-          a22(k) = SQRT( a22(k) )
-       ELSE
-          a22(k) = 0.0_wp
-       ENDIF
-
-!
-!--    a31, a32, a33 must be calculated with interpolated a11, a21, a22 (d11,
-!--    d21, d22) because of different vertical grid
-       IF ( k .le. nzt )  THEN
-          d11 = 0.5_wp * ( r11(k) + r11(k+1) )
-          IF ( d11 > 0.0_wp )  THEN
-             d11 = SQRT( d11 )
-             d21 = ( 0.5_wp * ( r21(k) + r21(k+1) ) ) / d11
-             a31(k) = r31(k) / d11
-          ELSE
-             d21 = 0.0_wp
-             a31(k) = 0.0_wp
-          ENDIF
-
-          d22 = 0.5_wp * ( r22(k) + r22(k+1) ) - d21 ** 2
-          IF ( d22 > 0.0_wp )  THEN
-             a32(k) = ( r32(k) - d21 * a31(k) ) / SQRT( d22 )
-          ELSE
-             a32(k) = 0.0_wp
-          ENDIF
-
-          a33(k) = r33(k) - a31(k) ** 2 - a32(k) ** 2
-          IF ( a33(k) > 0.0_wp )  THEN
-             a33(k) = SQRT( a33(k) )
-          ELSE
-             a33(k) = 0.0_wp
-          ENDIF
-       ELSE
-          a31(k) = a31(k-1)
-          a32(k) = a32(k-1)
-          a33(k) = a33(k-1)
-       ENDIF
-
-    ENDDO
+    
 !
 !-- Define the size of the filter functions and allocate them.
@@ -961,5 +909,5 @@
     ENDDO
 
-END SUBROUTINE stg_filter_func
+ END SUBROUTINE stg_filter_func
 
 
@@ -977,5 +925,5 @@
 
 
-    NAMELIST /stg_par/   use_syn_turb_gen
+    NAMELIST /stg_par/  dt_stg_adjust, dt_stg_call, use_syn_turb_gen
 
     line = ' '
@@ -1108,14 +1056,9 @@
     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
-
-
     CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
 
 !
 !-- Calculate time step which is needed for filter functions
-    dt_stg = dt_3d * weight_substep(intermediate_timestep_count)
-
+    dt_stg = MAX( dt_3d, dt_stg_call ) !* weight_substep(intermediate_timestep_count)
 !
 !-- Initial value of fu, fv, fw
@@ -1145,4 +1088,5 @@
        velocity_seed_initialized = .TRUE.
     ENDIF
+    
 !
 !-- New set of fu, fv, fw
@@ -1150,27 +1094,28 @@
     CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_left )
     CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_left )
-
+    
     IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent         &
-                             .AND.  .NOT.  rans_mode ) )  THEN
-!
-!--       Generate turbulence at right boundary
-          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_right )
-          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_right )
-          CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_right )
-!
-!--       Generate turbulence at north boundary
-          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_north )
-          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_north )
-          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_north )
-!
-!--       Generate turbulence at south boundary
-          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_south )
-          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_south )
-          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_south )
+                          .AND.  .NOT.  rans_mode ) )  THEN
+!
+!--    Generate turbulence at right boundary
+       CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_right )
+       CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_right )
+       CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_right )
+!
+!--    Generate turbulence at north boundary
+       CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_north )
+       CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_north )
+       CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_north )
+!
+!--    Generate turbulence at south boundary
+       CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_south )
+       CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_south )
+       CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_south )
     ENDIF
+    
 !
 !-- Turbulence generation at left and or right boundary
     IF ( myidx == id_stg_left  .OR.  myidx == id_stg_right )  THEN
-
+    
        DO  j = nysg, nyng
           DO  k = nzb, nzt + 1
@@ -1180,5 +1125,5 @@
                 fu_yz(k,j) = fuo_yz(k,j)
              ELSE
-                fu_yz(k,j) = fu_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) +     &
+                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
@@ -1187,5 +1132,5 @@
                 fv_yz(k,j) = fvo_yz(k,j)
              ELSE
-                fv_yz(k,j) = fv_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) +     &
+                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
@@ -1194,5 +1139,5 @@
                 fw_yz(k,j) = fwo_yz(k,j)
              ELSE
-                fw_yz(k,j) = fw_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) +     &
+                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
@@ -1204,21 +1149,16 @@
                 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
-                dist_yz(k,j,1) = a11(k) * fu_yz(k,j)
+                dist_yz(k,j,1) = MIN( a11(k) * fu_yz(k,j), 3.0_wp )
                 !experimental test of 1.2
-                dist_yz(k,j,2) = ( SQRT( a22(k) / MAXVAL(a22) )                &
-                                         * 1.2_wp )                            &
-                                       * (   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_yz(k,j)               &
-                                           + a32(k) * fv_yz(k,j)               &
-                                           + a33(k) * fw_yz(k,j) )
-                ! Calculation for pt and e not yet implemented
-!                 dist_yz(k,j,4) = 0.0_wp
-!                 dist_yz(k,j,5) = 0.0_wp
+                dist_yz(k,j,2) = MIN( ( SQRT( a22(k) / MAXVAL(a22) )           &
+                                      * 1.2_wp )                               &
+                                     * (   a21(k) * fu_yz(k,j)                 &
+                                         + a22(k) * fv_yz(k,j) ), 3.0_wp )
+                dist_yz(k,j,3) = MIN( ( SQRT(a33(k) / MAXVAL(a33) )            &
+                                      * 1.3_wp )                               &
+                                    * (   a31(k) * fu_yz(k,j)                  &
+                                        + a32(k) * fv_yz(k,j)                  &
+                                        + a33(k) * fw_yz(k,j) ), 3.0_wp )
              ENDIF
 
@@ -1241,6 +1181,6 @@
 
 #if defined( __parallel )
-          CALL MPI_ALLREDUCE( mc_factor_l, mc_factor,  &
-                              1, MPI_REAL, MPI_SUM, comm1dy, ierr )
+          CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, 1, MPI_REAL, MPI_SUM,    &
+                              comm1dy, ierr )
 #else
           mc_factor = mc_factor_l
@@ -1298,10 +1238,10 @@
 !--       Add disturbances
           IF ( myidx == id_stg_left  )  THEN
-
              DO  j = nys, nyn
-                DO  k = nzb, nzt
+                DO  k = nzb+1, nzt
                    u(k,j,0) = ( u(k,j,0) + dist_yz(k,j,1) )                    &
                        * mc_factor * MERGE( 1.0_wp, 0.0_wp,                    &
                                             BTEST( wall_flags_0(k,j,0), 1 ) )
+                   u(k,j,-1) = u(k,j,0)
                    v(k,j,-1)  = ( v(k,j,-1)  + dist_yz(k,j,2)  )               &
                        * mc_factor * MERGE( 1.0_wp, 0.0_wp,                    &
@@ -1314,7 +1254,6 @@
           ENDIF
           IF ( myidx == id_stg_right  )  THEN
-
              DO  j = nys, nyn
-                DO  k = nzb, nzt
+                DO  k = nzb+1, nzt
                    u(k,j,nxr+1) = ( u(k,j,nxr+1) + dist_yz(k,j,1) )            &
                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
@@ -1335,5 +1274,5 @@
 !-- 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
@@ -1351,12 +1290,12 @@
              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) ) )
+                      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) ) )
+                      fwo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
              ENDIF
 !
@@ -1369,19 +1308,20 @@
 
              ELSE
-                dist_xz(k,i,1) = a11(k) * fu_xz(k,i)
+                dist_xz(k,i,1) = MIN( a11(k) * fu_xz(k,i), 3.0_wp )
                 !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) )
+                dist_xz(k,i,2) = MIN( ( SQRT( a22(k) / MAXVAL(a22) )           &
+                                      * 1.2_wp )                               &
+                                     * (   a21(k) * fu_xz(k,i)                 &
+                                         + a22(k) * fv_xz(k,i) ), 3.0_wp )
+                dist_xz(k,i,3) = MIN( ( SQRT(a33(k) / MAXVAL(a33) )            &
+                                      * 1.3_wp )                               &
+                                    * (   a31(k) * fu_xz(k,i)                  &
+                                        + a32(k) * fv_xz(k,i)                  &
+                                        + a33(k) * fw_xz(k,i) ), 3.0_wp )
              ENDIF
 
           ENDDO
        ENDDO
+
 !
 !--    Mass flux correction following Kim et al. (2013)
@@ -1425,5 +1365,5 @@
 
           DO  i = nxl, nxr
-             DO  k = nzb, nzt
+             DO  k = nzb+1, nzt
                 u(k,-1,i) = ( u(k,-1,i) + dist_xz(k,i,1) )                     &
                       * mc_factor * MERGE( 1.0_wp, 0.0_wp,                     &
@@ -1432,4 +1372,5 @@
                       * mc_factor * MERGE( 1.0_wp, 0.0_wp,                     &
                                            BTEST( wall_flags_0(k,0,i), 2 ) )
+                v(k,-1,i) = v(k,0,i)
                 w(k,-1,i) = ( w(k,-1,i) + dist_xz(k,i,3)  )                    &
                       * mc_factor * MERGE( 1.0_wp, 0.0_wp,                     &
@@ -1441,13 +1382,13 @@
 
           DO  i = nxl, nxr
-             DO  k = nzb, nzt
-                u(k,nyn+1,i) = ( u(k,nyn+1,i) + dist_xz(k,i,1) )                &
-                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                       &
+             DO  k = nzb+1, nzt
+                u(k,nyn+1,i) = ( u(k,nyn+1,i) + dist_xz(k,i,1) )               &
+                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
                                           BTEST( wall_flags_0(k,nyn+1,i), 1 ) )
-                v(k,nyn+1,i) = ( v(k,nyn+1,i) + dist_xz(k,i,1) )                &
-                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                       &
+                v(k,nyn+1,i) = ( v(k,nyn+1,i) + dist_xz(k,i,2) )               &
+                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
                                           BTEST( wall_flags_0(k,nyn+1,i), 2 ) )
-                w(k,nyn+1,i) = ( w(k,nyn+1,i) + dist_xz(k,i,1) )                &
-                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                       &
+                w(k,nyn+1,i) = ( w(k,nyn+1,i) + dist_xz(k,i,3) )               &
+                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
                                           BTEST( wall_flags_0(k,nyn+1,i), 3 ) )
              ENDDO
@@ -1455,4 +1396,7 @@
        ENDIF
     ENDIF
+!
+!-- Finally, set time counter for calling STG to zero
+    time_stg_call = 0.0_wp
 
     CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
@@ -1493,6 +1437,6 @@
     REAL(wp) :: rand_sigma_inv  !< inverse of stdev of random number
 
-    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(-merg:merg,nzb:nzt+1)    :: b_y     !< filter function in y-direction
+    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter function in z-direction
     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
@@ -1613,5 +1557,5 @@
     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_x2   !< n_x*2
     INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z2   !< n_z*2
 
@@ -1620,6 +1564,6 @@
     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(-merg:merg,nzb:nzt+1)    :: b_x     !< filter function in x-direction
+    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter function in z-direction
     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
@@ -1712,3 +1656,453 @@
  END SUBROUTINE stg_generate_seed_xz
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Parametrization of the Reynolds stress tensor, following the parametrization
+!> described in Rotach et al. (1996), which is applied in state-of-the-art
+!> dispserion modelling. Please note, the parametrization does not distinguish
+!> between along-wind and cross-wind turbulence. 
+!------------------------------------------------------------------------------!
+ SUBROUTINE parametrize_reynolds_stress
+
+    USE arrays_3d,                                                             &
+        ONLY:  zu
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) :: k            !< loop index in z-direction
+    
+    REAL(wp)     ::  zzi         !< ratio of z/zi
+    
+!
+!-- 
+    DO  k = nzb+1, nzt+1
+          
+       IF ( zu(k) <= zi_ribulk )  THEN
+!
+!--       Determine normalized height coordinate
+          zzi = zu(k) / zi_ribulk
+!
+!--       u'u' and v'v'. Assume isotropy. Note, add a small negative number
+!--       to the denominator, else the merge-function can crash if scale_l is 
+!--       zero. 
+          r11(k) = scale_us**2 * (                                             &
+                      MERGE( 0.35_wp * (                                       &
+                           - zi_ribulk / ( kappa * scale_l - 10E-4_wp )        &
+                                       )**( 2.0_wp / 3.0_wp ),                 &
+                             0.0_wp,                                           &
+                             scale_l < 0.0_wp )                                &
+                    + 5.0_wp - 4.0_wp * zzi                                    &
+                                 )
+                                 
+          r22(k) = r11(k)
+!
+!--       w'w'
+          r33(k) = scale_wm**2 * (                                             &
+                      1.5_wp * zzi**( 2.0_wp / 3.0_wp ) * EXP( -2.0_wp * zzi ) &
+                    + ( 1.7_wp - zzi ) * ( scale_us / scale_wm )**2            &                     
+                                 )
+!
+!--       u'w' and v'w'. Assume isotropy.
+          r31(k) = - scale_us**2 * (                                           &
+                         1.0_wp - EXP( 3.0_wp * ( zzi - 1.0_wp ) )             &
+                                   )
+                              
+          r32(k) = r31(k)
+!
+!--       For u'v' no parametrization exist so far - ?. For simplicity assume
+!--       a similar profile as for u'w'. 
+          r21(k) = r31(k)
+!
+!--    Above the boundary layer, assmume laminar flow conditions.
+       ELSE
+          r11(k) = 10E-8_wp
+          r22(k) = 10E-8_wp
+          r33(k) = 10E-8_wp
+          r21(k) = 10E-8_wp
+          r31(k) = 10E-8_wp
+          r32(k) = 10E-8_wp
+       ENDIF
+!        write(9,*) zu(k), r11(k), r33(k), r31(k), zi_ribulk, scale_us, scale_wm, scale_l
+    ENDDO
+
+!
+!-- Set bottom boundary condition    
+    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)
+    
+
+ END SUBROUTINE parametrize_reynolds_stress 
+ 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculate the coefficient matrix from the Lund rotation. 
+!------------------------------------------------------------------------------!
+ SUBROUTINE calc_coeff_matrix
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) :: k   !< loop index in z-direction
+    
+!
+!-- Calculate coefficient matrix. Split loops to allow for loop vectorization.
+    DO  k = nzb+1, nzt+1
+       IF ( r11(k) > 0.0_wp )  THEN
+          a11(k) = SQRT( r11(k) ) 
+          a21(k) = r21(k) / a11(k)
+          a31(k) = r31(k) / a11(k)
+       ELSE
+          a11(k) = 10E-8_wp
+          a21(k) = 10E-8_wp
+          a31(k) = 10E-8_wp
+       ENDIF
+    ENDDO
+    DO  k = nzb+1, nzt+1
+       a22(k) = r22(k) - a21(k)**2 
+       IF ( a22(k) > 0.0_wp )  THEN
+          a22(k) = SQRT( a22(k) )
+          a32(k) = r32(k) - a21(k) * a31(k) / a22(k)
+       ELSE 
+          a22(k) = 10E-8_wp
+          a32(k) = 10E-8_wp
+       ENDIF
+    ENDDO
+    DO  k = nzb+1, nzt+1
+       a33(k) = r33(k) - a31(k)**2 - a32(k)**2
+       IF ( a33(k) > 0.0_wp )  THEN
+          a33(k) =  SQRT( a33(k) )
+       ELSE
+          a33(k) = 10E-8_wp
+       ENDIF
+    ENDDO   
+!
+!-- Set bottom boundary condition
+    a11(nzb) = a11(nzb+1)
+    a22(nzb) = a22(nzb+1)
+    a21(nzb) = a21(nzb+1)
+    a33(nzb) = a33(nzb+1)    
+    a31(nzb) = a31(nzb+1)
+    a32(nzb) = a32(nzb+1)    
+
+ END SUBROUTINE calc_coeff_matrix
+ 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This routine controls the re-adjustment of the turbulence statistics used
+!> for generating turbulence at the lateral boundaries.  
+!------------------------------------------------------------------------------!
+ SUBROUTINE stg_adjust
+
+    IMPLICIT NONE
+    
+    INTEGER(iwp) ::  k
+
+    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
+!
+!-- Compute mean boundary layer height according to Richardson-Bulk
+!-- criterion using the inflow profiles. Further velocity scale as well as
+!-- mean friction velocity are calculated.
+    CALL calc_scaling_variables
+!
+!-- Set length and time scales depending on boundary-layer height
+    CALL calc_length_and_time_scale
+!
+!-- Parametrize Reynolds-stress tensor, diagonal elements as well
+!-- as r21 (v'u'), r31 (w'u'), r32 (w'v'). Parametrization follows
+!-- Rotach et al. (1996) and is based on boundary-layer depth, 
+!-- friction velocity and velocity scale. 
+    CALL parametrize_reynolds_stress
+!
+!-- Calculate coefficient matrix from Reynolds stress tensor  
+!-- (Lund rotation)
+    CALL calc_coeff_matrix
+!
+!-- Determine filter functions on basis of updated length scales
+    CALL stg_filter_func( nux, bux ) !filter ux
+    CALL stg_filter_func( nuy, buy ) !filter uy
+    CALL stg_filter_func( nuz, buz ) !filter uz
+    CALL stg_filter_func( nvx, bvx ) !filter vx
+    CALL stg_filter_func( nvy, bvy ) !filter vy
+    CALL stg_filter_func( nvz, bvz ) !filter vz
+    CALL stg_filter_func( nwx, bwx ) !filter wx
+    CALL stg_filter_func( nwy, bwy ) !filter wy
+    CALL stg_filter_func( nwz, bwz ) !filter wz
+!
+!-- Reset time counter for controlling next adjustment to zero
+    time_stg_adjust = 0.0_wp
+    
+    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
+    
+ END SUBROUTINE stg_adjust 
+ 
+ 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculates turbuluent length and time scales if these are not available
+!> from measurements. 
+!------------------------------------------------------------------------------!
+ SUBROUTINE calc_length_and_time_scale
+
+    USE arrays_3d,                                                             &
+        ONLY:  dzw, ddzw, u_init, v_init, zu, zw
+
+    USE grid_variables,                                                        &
+        ONLY:  ddx, ddy, dx, dy
+        
+    USE statistics,                                                            &
+        ONLY:  hom
+
+    IMPLICIT NONE
+
+
+    INTEGER(iwp) :: k            !< loop index in z-direction
+
+    REAL(wp) ::  length_scale    !< typical length scale
+    
+!
+!-- In initial call the boundary-layer depth can be zero. This case, set 
+!-- minimum value for boundary-layer depth, to setup length scales correctly.
+    zi_ribulk = MAX( zi_ribulk, zw(nzb+2) )
+!
+!-- Set-up default turbulent length scales. From the numerical point of
+!-- view the imposed perturbations should not be immediately dissipated
+!-- by the numerics. The numerical dissipation, however, acts on scales
+!-- up to 8 x the grid spacing. For this reason, set the turbulence 
+!-- length scale to 8 time the grid spacing. Further, above the boundary
+!-- layer height, set turbulence lenght scales to zero (equivalent to not
+!-- imposing any perturbations) in order to save computational costs. 
+!-- Typical time scales are derived by assuming Taylors's hypothesis, 
+!-- using the length scales and the mean profiles of u- and v-component.
+    DO  k = nzb+1, nzt+1
+    
+       length_scale = 8.0_wp * MIN( dx, dy, dzw(k) )
+          
+       IF ( zu(k) <= zi_ribulk )  THEN  
+!
+!--       Assume isotropic turbulence length scales
+          nux(k) = MAX( INT( length_scale * ddx     ), 1 )
+          nuy(k) = MAX( INT( length_scale * ddy     ), 1 )
+          nuz(k) = MAX( INT( length_scale * ddzw(k) ), 1 )
+          nvx(k) = MAX( INT( length_scale * ddx     ), 1 )
+          nvy(k) = MAX( INT( length_scale * ddy     ), 1 )
+          nvz(k) = MAX( INT( length_scale * ddzw(k) ), 1 )
+          nwx(k) = MAX( INT( length_scale * ddx     ), 1 )
+          nwy(k) = MAX( INT( length_scale * ddy     ), 1 )
+          nwz(k) = MAX( INT( length_scale * ddzw(k) ), 1 )
+!
+!--       Limit time scales, else they become very larger for low wind speed, 
+!--       imposing long-living inflow perturbations which in turn propagate
+!--       further into the model domain. Use u_init and v_init to calculate
+!--       the time scales, which will be equal to the inflow profiles, both, 
+!--       in offline nesting mode or in dirichlet/radiation mode. 
+          tu(k)  = MIN( dt_stg_adjust, length_scale /                          &
+                                  ( ABS( u_init(k) ) + 0.1_wp ) )
+          tv(k)  = MIN( dt_stg_adjust, length_scale /                          &
+                                  ( ABS( v_init(k) ) + 0.1_wp ) )
+!
+!--       Time scale of w-component is a mixture from u- and v-component. 
+          tw(k)  = SQRT( tu(k)**2 + tv(k)**2 )
+!
+!--    Above the boundary layer length scales are zero, i.e. imposed turbulence
+!--    is not correlated in space and time, just white noise. This saves
+!--    computations power. 
+       ELSE
+          nux(k) = 0.0_wp
+          nuy(k) = 0.0_wp
+          nuz(k) = 0.0_wp
+          nvx(k) = 0.0_wp
+          nvy(k) = 0.0_wp
+          nvz(k) = 0.0_wp
+          nwx(k) = 0.0_wp
+          nwy(k) = 0.0_wp
+          nwz(k) = 0.0_wp
+          
+          tu(k)  = 0.0_wp
+          tv(k)  = 0.0_wp
+          tw(k)  = 0.0_wp
+       ENDIF
+    ENDDO
+!
+!-- Set bottom boundary condition for the length and time scales 
+    nux(nzb) = nux(nzb+1)
+    nuy(nzb) = nuy(nzb+1)
+    nuz(nzb) = nuz(nzb+1)
+    nvx(nzb) = nvx(nzb+1)
+    nvy(nzb) = nvy(nzb+1)
+    nvz(nzb) = nvz(nzb+1)
+    nwx(nzb) = nwx(nzb+1)
+    nwy(nzb) = nwy(nzb+1)
+    nwz(nzb) = nwz(nzb+1)
+
+    tu(nzb)  = tu(nzb+1)
+    tv(nzb)  = tv(nzb+1)
+    tw(nzb)  = tw(nzb+1)
+
+
+ END SUBROUTINE calc_length_and_time_scale 
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculate scaling variables which are used for turbulence parametrization
+!> according to Rotach et al. (1996). Scaling variables are: friction velocity,
+!> boundary-layer depth, momentum velocity scale, and Obukhov length. 
+!------------------------------------------------------------------------------!
+ SUBROUTINE calc_scaling_variables
+
+
+    USE control_parameters,                                                    &
+        ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, &
+               humidity, pt_surface
+             
+    USE indices,                                                               &
+        ONLY:  nx, ny  
+        
+    USE surface_mod,                                                           &
+        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) :: i            !< loop index in x-direction
+    INTEGER(iwp) :: j            !< loop index in y-direction
+    INTEGER(iwp) :: k            !< loop index in z-direction
+    INTEGER(iwp) :: k_ref        !< index in z-direction for reference height 
+    INTEGER(iwp) :: m            !< surface element index
+
+    REAL(wp) ::  friction_vel_l         !< mean friction veloctiy on subdomain
+    REAL(wp) ::  shf_mean               !< mean surface sensible heat flux
+    REAL(wp) ::  shf_mean_l             !< mean surface sensible heat flux on subdomain
+    REAL(wp) ::  u_int                  !< u-component
+    REAL(wp) ::  v_int                  !< v-component
+    REAL(wp) ::  vpt_surface            !< near-surface virtual potential temperature
+    REAL(wp) ::  w_convective           !< convective velocity scale
+    REAL(wp) ::  z0_mean                !< mean roughness length
+    REAL(wp) ::  z0_mean_l              !< mean roughness length on subdomain
+   
+!
+!-- Mean friction velocity and velocity scale. Therefore, 
+!-- pre-calculate mean roughness length and surface sensible heat flux
+!-- in the model domain, which are further used to estimate friction 
+!-- velocity and velocity scale. Note, for z0 linear averaging is applied,
+!-- even though this is known to unestimate the effective roughness.
+!-- This need to be revised later.
+    z0_mean_l  = 0.0_wp
+    shf_mean_l = 0.0_wp
+    DO  m = 1, surf_def_h(0)%ns
+       z0_mean_l  = z0_mean_l  + surf_def_h(0)%z0(m)
+       shf_mean_l = shf_mean_l + surf_def_h(0)%shf(m)
+    ENDDO    
+    DO  m = 1, surf_lsm_h%ns
+       z0_mean_l  = z0_mean_l  + surf_lsm_h%z0(m)
+       shf_mean_l = shf_mean_l + surf_lsm_h%shf(m)
+    ENDDO
+    DO  m = 1, surf_usm_h%ns
+       z0_mean_l  = z0_mean_l  + surf_usm_h%z0(m)
+       shf_mean_l = shf_mean_l + surf_usm_h%shf(m)
+    ENDDO
+    
+#if defined( __parallel )
+    CALL MPI_ALLREDUCE( z0_mean_l, z0_mean, 1, MPI_REAL, MPI_SUM,              &
+                        comm2d, ierr )
+    CALL MPI_ALLREDUCE( shf_mean_l, shf_mean, 1, MPI_REAL, MPI_SUM,            &
+                        comm2d, ierr )
+#else
+    z0_mean  = z0_mean_l
+    shf_mean = shf_mean_l
+#endif
+    z0_mean  = z0_mean  / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
+    shf_mean = shf_mean / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
+        
+!
+!-- Note, Inifor does not use logarithmic interpolation of the 
+!-- velocity components near the ground, so that near-surface
+!-- wind speed and thus the friction velocity is overestimated. 
+!-- However, friction velocity is used for turbulence 
+!-- parametrization, so that more physically meaningful values are important.
+!-- Hence, derive friction velocity from wind speed at a reference height.
+!-- For a first guess use 20 m, which is in the range of the first 
+!-- COSMO vertical level. 
+    k_ref = MINLOC( ABS( zu - cosmo_ref ), DIM = 1 ) - 1
+    
+    friction_vel_l = 0.0_wp
+    IF ( bc_dirichlet_l  .OR.  bc_dirichlet_r )  THEN
+       
+       i = MERGE( -1, nxr + 1, bc_dirichlet_l )
+
+       DO  j = nys, nyn
+!
+!--       Note, in u- and v- component the imposed perturbations 
+!--       from the STG are already included. Check whether this 
+!--       makes any difference compared to using the pure-mean 
+!--       inflow profiles. 
+          u_int = MERGE( u(k_ref,j,i+1), u(k_ref,j,i), bc_dirichlet_l )
+          v_int = v(k_ref,j,i)
+!
+!--       Calculate friction velocity and sum-up. Therefore, assume 
+!--       neutral condtions.
+          friction_vel_l = friction_vel_l + kappa *                            &
+                            SQRT( u_int * u_int + v_int * v_int )  /           &
+                            LOG( zu(k_ref) / z0_mean )
+          
+       ENDDO
+       
+    ENDIF
+    
+    IF ( bc_dirichlet_s  .OR.  bc_dirichlet_n )  THEN
+       
+       j = MERGE( -1, nyn + 1, bc_dirichlet_s )
+    
+       DO  i = nxl, nxr
+          
+          u_int = u(k_ref,j,i)
+          v_int = MERGE( v(k_ref,j+1,i), v(k_ref,j,i), bc_dirichlet_s )
+ 
+          friction_vel_l = friction_vel_l + kappa *                            &
+                            SQRT( u_int * u_int + v_int * v_int )  /           &
+                            LOG( zu(k_ref) / z0_mean )                       
+          
+       ENDDO
+       
+    ENDIF
+    
+#if defined( __parallel )
+    CALL MPI_ALLREDUCE( friction_vel_l, scale_us, 1, MPI_REAL, MPI_SUM,        &
+                        comm2d, ierr )
+#else
+    scale_us = friction_vel_l
+#endif
+    scale_us = scale_us / REAL( 2 * nx + 2 * ny, KIND = wp )
+    
+!
+!-- Compute mean Obukhov length
+    IF ( shf_mean > 0.0_wp )  THEN
+       scale_l = - pt_surface / ( kappa * g ) * scale_us**3 / shf_mean
+    ELSE
+       scale_l = 0.0_wp
+    ENDIF
+    
+!
+!-- Compute mean convective velocity scale. Note, in case the mean heat flux 
+!-- is negative, set convective velocity scale to zero.
+    IF ( shf_mean > 0.0_wp )  THEN
+       w_convective = ( g * shf_mean * zi_ribulk / pt_surface )**( 1.0_wp / 3.0_wp )
+    ELSE
+       w_convective = 0.0_wp
+    ENDIF
+!
+!-- Finally, in order to consider also neutral or stable stratification, 
+!-- compute momentum velocity scale from u* and convective velocity scale, 
+!-- according to Rotach et al. (1996). 
+    scale_wm = ( scale_us**3 + 0.6_wp * w_convective**3 )**( 1.0_wp / 3.0_wp )
+   
+ END SUBROUTINE calc_scaling_variables
+
  END MODULE synthetic_turbulence_generator_mod
Index: /palm/trunk/SOURCE/time_integration.f90
===================================================================
--- /palm/trunk/SOURCE/time_integration.f90	(revision 3346)
+++ /palm/trunk/SOURCE/time_integration.f90	(revision 3347)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! - offline nesting separated from large-scale forcing module
+! - changes for synthetic turbulence generator
+! 
+! 3343 2018-10-15 10:38:52Z suehring
 ! - Formatting, clean-up, comments (kanani)
 ! - Added CALL to chem_emissions_setup (Russo)
@@ -459,12 +463,15 @@
 
     USE lsf_nudging_mod,                                                       &
-        ONLY:  calc_tnudge, ls_forcing_surf, ls_forcing_vert, nudge_ref,       &
-               lsf_nesting_offline, lsf_nesting_offline_mass_conservation
+        ONLY:  calc_tnudge, ls_forcing_surf, ls_forcing_vert, nudge_ref
 
     USE multi_agent_system_mod,                                                &
         ONLY:  agents_active, multi_agent_system
-
+        
+    USE nesting_offl_mod,                                                      &
+        ONLY:  nesting_offl_bc, nesting_offl_mass_conservation
+        
     USE netcdf_data_input_mod,                                                 &
-        ONLY:  chem_emis, chem_emis_att, nest_offl, netcdf_data_input_lsf
+        ONLY:  chem_emis, chem_emis_att, nest_offl,                            &
+               netcdf_data_input_offline_nesting
 
     USE ocean_mod,                                                             &
@@ -515,5 +522,7 @@
 
     USE synthetic_turbulence_generator_mod,                                    &
-        ONLY:  stg_main, use_syn_turb_gen
+        ONLY:  dt_stg_call, dt_stg_adjust, parametrize_inflow_turbulence,      &
+               stg_adjust, stg_main, time_stg_adjust, time_stg_call,           &
+               use_syn_turb_gen
 
     USE user_actions_mod,                                                      &
@@ -638,5 +647,5 @@
        IF ( nesting_offline )  THEN
           IF ( nest_offl%time(nest_offl%tind_p) <= time_since_reference_point )&
-             CALL netcdf_data_input_lsf
+             CALL netcdf_data_input_offline_nesting
        ENDIF
 
@@ -918,8 +927,4 @@
 !--       Impose a turbulent inflow using the recycling method
           IF ( turbulent_inflow )  CALL  inflow_turbulence
-
-!
-!--       Impose a turbulent inflow using synthetic generated turbulence
-          IF ( use_syn_turb_gen )  CALL  stg_main
 
 !
@@ -960,11 +965,19 @@
 !--       boundaries. 
           IF ( nesting_offline  .AND.  intermediate_timestep_count ==          &
-                                       intermediate_timestep_count_max  )  THEN
-             CALL lsf_nesting_offline
-!
-!--          Moreover, ensure mass conservation
-             CALL lsf_nesting_offline_mass_conservation
-          ENDIF
-
+                                       intermediate_timestep_count_max  )      &
+             CALL nesting_offl_bc
+!
+!--       Impose a turbulent inflow using synthetic generated turbulence, 
+!--       only once per time step.
+          IF ( use_syn_turb_gen  .AND.  time_stg_call >= dt_stg_call  .AND.    &
+             intermediate_timestep_count == intermediate_timestep_count_max )  THEN! &
+             CALL stg_main
+          ENDIF
+!
+!--       Ensure mass conservation. This need to be done after imposing
+!--       synthetic turbulence and top boundary condition for pressure is set to
+!--       Neumann conditions. 
+!--       Is this also required in case of Dirichlet?
+          IF ( nesting_offline )  CALL nesting_offl_mass_conservation
 !
 !--       Reduce the velocity divergence via the equation for perturbation
@@ -1213,4 +1226,15 @@
        time_dopr_listing          = time_dopr_listing        + dt_3d
        time_run_control   = time_run_control + dt_3d
+!
+!--    In case of synthetic turbulence generation and parametrized turbulence
+!--    information, update the time counter and if required, adjust the 
+!--    STG to new atmospheric conditions.
+       IF ( use_syn_turb_gen  )  THEN
+          IF ( parametrize_inflow_turbulence )  THEN
+             time_stg_adjust = time_stg_adjust + dt_3d
+             IF ( time_stg_adjust >= dt_stg_adjust )  CALL stg_adjust
+          ENDIF
+          time_stg_call = time_stg_call + dt_3d
+       ENDIF
 
 !
@@ -1472,5 +1496,6 @@
 
     ENDDO   ! time loop
-
+    
+!
 !-- Vertical nesting: Deallocate variables initialized for vertical nesting    
     IF ( vnest_init )  CALL vnest_deallocate
Index: /palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3346)
+++ /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3347)
@@ -28,4 +28,8 @@
 ! -----------------
 ! $Id$
+! Enable USM initialization with default building parameters in case no static
+! input file exist. 
+! 
+! 3343 2018-10-15 10:38:52Z suehring
 ! Add output variables usm_rad_pc_inlw, usm_rad_pc_insw*
 ! 
@@ -6902,4 +6906,6 @@
         INTEGER(iwp)                                          :: wtc
         REAL(wp), DIMENSION(n_surface_params)                 :: wtp
+        
+        LOGICAL                                               :: ascii_file = .FALSE.
     
         INTEGER(iwp), DIMENSION(0:17, nysg:nyng, nxlg:nxrg)   :: usm_par
@@ -6922,4 +6928,11 @@
         IF ( building_type_f%from_file  .OR.  building_pars_f%from_file )      &
            RETURN
+!
+!--     Check if ASCII input file exists. If not, return and initialize USM
+!--     with default settings.
+        INQUIRE( FILE = 'SURFACE_PARAMETERS' // coupling_char,                 &
+                 EXIST = ascii_file )
+                 
+        IF ( .NOT. ascii_file )  RETURN
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
