Changeset 3347

Timestamp:

Oct 15, 2018 2:21:08 PM (6 years ago)

Author:

suehring

Message:

Offline nesting revised and separated from large_scale_forcing_mod; Time-dependent synthetic turbulence generator; bugfixes in USM/LSM radiation model and init_pegrid

Location:

palm/trunk/SOURCE

Files:

• Makefile (modified) (10 diffs)
• boundary_conds.f90 (modified) (3 diffs)
• check_parameters.f90 (modified) (2 diffs)
• header.f90 (modified) (3 diffs)
• init_3d_model.f90 (modified) (6 diffs)
• init_pegrid.f90 (modified) (2 diffs)
• land_surface_model_mod.f90 (modified) (2 diffs)
• large_scale_forcing_nudging_mod.f90 (modified) (10 diffs)
• nesting_offl_mod.f90 (added)
• netcdf_data_input_mod.f90 (modified) (6 diffs)
• parin.f90 (modified) (6 diffs)
• pres.f90 (modified) (4 diffs)
• synthetic_turbulence_generator_mod.f90 (modified) (43 diffs)
• time_integration.f90 (modified) (8 diffs)
• urban_surface_mod.f90 (modified) (3 diffs)

palm/trunk/SOURCE/Makefile

@@ -25 +25 @@
 # -----------------
 # $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 +572 @@
         netcdf_data_input_mod.f90 \
         netcdf_interface_mod.f90 \
+        nesting_offl_mod.f90 \
         ocean_mod.f90 \
         outflow_turbulence.f90 \
@@ -1005 +1012 @@
         modules.o \
         netcdf_interface_mod.o \
+        nesting_offl_mod.o \
         ocean_mod.o \
         plant_canopy_model_mod.o \
@@ -1036 +1044 @@
         netcdf_data_input_mod.o \
         netcdf_interface_mod.o \
+        nesting_offl_mod.o \
         ocean_mod.o \
         plant_canopy_model_mod.o \
@@ -1045 +1054 @@
         surface_layer_fluxes_mod.o \
         surface_mod.o \
+        synthetic_turbulence_generator_mod.o \
         turbulence_closure_mod.o \
         urban_surface_mod.o \
@@ -1112 +1122 @@
         mod_kinds.o \
         modules.o \
-        netcdf_data_input_mod.o \
         surface_mod.o
 local_stop.o: \
@@ -1288 +1297 @@
         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 +1350 @@
         modules.o \
         multi_agent_system_mod.o \
+        nesting_offl_mod.o \
         netcdf_interface_mod.o \
         ocean_mod.o \
@@ -1577 +1593 @@
         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 +1618 @@
         modules.o \
         multi_agent_system_mod.o \
+        nesting_offl_mod.o \
         ocean_mod.o \
         pmc_interface_mod.o \

palm/trunk/SOURCE/boundary_conds.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Bugfix in setting boundary conditions in offline nesting
+! 
+! 3341 2018-10-15 10:31:27Z suehring
 ! changes concerning modularization of ocean option
 ! 
@@ -293 +296 @@
 !
 !-- 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 +732 @@
 !-- 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

palm/trunk/SOURCE/check_parameters.f90

@@ -25 +25 @@
 ! -----------------
 ! $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 +1957 @@
        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

palm/trunk/SOURCE/header.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Header output concerning offline nesting
+! 
+! 3343 2018-10-15 10:38:52Z suehring
 ! call of ocean_header
 ! 
@@ -438 +441 @@
     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 +909 @@
 !-- Large scale forcing and nudging
     IF ( large_scale_forcing )  CALL lsf_nudging_header( io )
+    
+!
+!-- Offline nesting 
+    IF ( nesting_offline )  CALL nesting_offl_header( io )
 
 !

palm/trunk/SOURCE/init_3d_model.f90

@@ -25 +25 @@
 ! -----------------
 ! $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 +571 @@
         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 +612 @@
 
     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 +2164 @@
     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 +2175 @@
        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 +2295 @@
        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

palm/trunk/SOURCE/init_pegrid.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Bugfix in setting number of ghost layers in neutral case 
+! 
+! 3341 2018-10-15 10:31:27Z suehring
 ! unused variables removed
 ! 
@@ -723 +726 @@
 ! 
 !-- 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

palm/trunk/SOURCE/land_surface_model_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Assign real value instead of integer
+! 
+! 3341 2018-10-15 10:31:27Z suehring
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -2781 +2784 @@
                          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
                       

palm/trunk/SOURCE/large_scale_forcing_nudging_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $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 +90 @@
 
     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 +96 @@
                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 +112 @@
     USE kinds
 
-    USE netcdf_data_input_mod,                                                 &
-        ONLY:  nest_offl
-
     USE pegrid
 
@@ -120 +122 @@
     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 +156 @@
            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 +176 @@
  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 +456 @@
            ONLY:  bc_lr_cyc, bc_ns_cyc
 
-       USE netcdf_data_input_mod,                                              &
-           ONLY:  netcdf_data_input_lsf  
-
        IMPLICIT NONE
 
@@ -967 +485 @@
        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 +591 @@
           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

palm/trunk/SOURCE/netcdf_data_input_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Subroutine renamed
+! 
+! 3257 2018-09-17 17:11:46Z suehring
 ! (from branch resler)
 ! Formatting
@@ -678 +681 @@
     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 +738 @@
            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 +3081 @@
 !> (COSMO) by Inifor.
 !------------------------------------------------------------------------------!
-    SUBROUTINE netcdf_data_input_lsf
+    SUBROUTINE netcdf_data_input_offline_nesting
 
        USE control_parameters,                                                 &
@@ -3157 +3160 @@
 !
 !--    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 +3351 @@
        CALL cpu_log( log_point_s(86), 'NetCDF input forcing', 'stop' )
 
-    END SUBROUTINE netcdf_data_input_lsf
+    END SUBROUTINE netcdf_data_input_offline_nesting
 
 

palm/trunk/SOURCE/parin.f90

@@ -25 +25 @@
 ! -----------------
 ! $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 +481 @@
         ONLY:  mas_parin
 
+    USE nesting_offl_mod,                                                      &
+        ONLY:  nesting_offl_parin
+        
     USE netcdf_interface,                                                      &
         ONLY:  netcdf_data_format, netcdf_deflate, netcdf_precision
@@ -570 +577 @@
              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 +650 @@
              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 +885 @@
           CALL gust_parin
           CALL mas_parin
+          CALL nesting_offl_parin
           CALL ocean_parin
           CALL pcm_parin
@@ -979 +987 @@
            !  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
 

palm/trunk/SOURCE/pres.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Bugfixes in offline nesting.
+! Add comment. 
+! 
+! 3341 2018-10-15 10:31:27Z suehring
 ! Rename variables for boundary flags and nesting
 ! 
@@ -169 +173 @@
                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 +391 @@
     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 +724 @@
        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

palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $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 +128 @@
 ! Authors:
 ! --------
-! @author Tobias Gronemeier, Atsushi Inagaki, Micha Gryschka, Christoph Knigge
+! @author Tobias Gronemeier, Matthias Suehring, Atsushi Inagaki, Micha Gryschka, Christoph Knigge
+! 
 !
 !
@@ -133 +140 @@
 !> 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 +157 @@
 
     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 +167 @@
 
     USE cpulog,                                                                &
-        ONLY:  cpu_log, log_point, log_point_s
+        ONLY:  cpu_log, log_point
 
     USE indices,                                                               &
@@ -170 +178 @@
 #endif
 
+    USE nesting_offl_mod,                                                      &
+        ONLY:  zi_ribulk
+
     USE pegrid,                                                                &
         ONLY:  comm1dx, comm1dy, comm2d, ierr, myidx, myidy, pdims
@@ -184 +195 @@
 
 
-    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 +214 @@
     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 +353 @@
 !
 !-- 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 +411 @@
              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 +418 @@
           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 +463 @@
        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 +490 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  dzw, ddzw, u_init, v_init, zu, zw
+        ONLY:  dzw, ddzw, u_init, v_init, zu
 
     USE control_parameters,                                                    &
@@ -460 +507 @@
     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 +531 @@
     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 +540 @@
     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 +546 @@
 
     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 +639 @@
     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 +714 @@
        
        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 +750 @@
        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 +909 @@
     ENDDO
 
-END SUBROUTINE stg_filter_func
+ END SUBROUTINE stg_filter_func
 
 
@@ -977 +925 @@
 
 
-    NAMELIST /stg_par/   use_syn_turb_gen
+    NAMELIST /stg_par/  dt_stg_adjust, dt_stg_call, use_syn_turb_gen
 
     line = ' '
@@ -1108 +1056 @@
     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 +1088 @@
        velocity_seed_initialized = .TRUE.
     ENDIF
+    
 !
 !-- New set of fu, fv, fw
@@ -1150 +1094 @@
     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 +1125 @@
                 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 +1132 @@
                 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 +1139 @@
                 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 +1149 @@
                 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 +1181 @@
 
 #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 +1238 @@
 !--       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 +1254 @@
           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 +1274 @@
 !-- 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 +1290 @@
              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 +1308 @@
 
              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 +1365 @@
 
           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 +1372 @@
                       * 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 +1382 @@
 
           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 +1396 @@
        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 +1437 @@
     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 +1557 @@
     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 +1564 @@
     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 +1656 @@
  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

palm/trunk/SOURCE/time_integration.f90

@@ -25 +25 @@
 ! -----------------
 ! $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 +463 @@
 
     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 +522 @@
 
     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 +647 @@
        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 +927 @@
 !--       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 +965 @@
 !--       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 +1226 @@
        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 +1496 @@
 
     ENDDO   ! time loop
-
+    
+!
 !-- Vertical nesting: Deallocate variables initialized for vertical nesting    
     IF ( vnest_init )  CALL vnest_deallocate

palm/trunk/SOURCE/urban_surface_mod.f90

@@ -28 +28 @@
 ! -----------------
 ! $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 +6906 @@
         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 +6928 @@
         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
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!