Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 3936)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 3937)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! Move initialization of synthetic turbulence generator behind initialization
+! of offline nesting. Remove call for stg_adjust, as this is now already done
+! in stg_init. 
+! 
+! 3900 2019-04-16 15:17:43Z suehring
 ! Fix problem with LOD = 2 initialization
 ! 
@@ -666,6 +671,5 @@
 
     USE synthetic_turbulence_generator_mod,                                    &
-        ONLY:  parametrize_inflow_turbulence, stg_adjust, stg_init,            &
-               use_syn_turb_gen
+        ONLY:  parametrize_inflow_turbulence, stg_init, use_syn_turb_gen
 
     USE surface_layer_fluxes_mod,                                              &
@@ -1305,7 +1309,4 @@
 !--       fluxes, etc. 
           CALL init_surfaces
-!
-!--       Initialize turbulence generator
-          IF( use_syn_turb_gen )  CALL stg_init
 
           CALL location_message( 'initializing with INIFOR', 'finished' )
@@ -1375,7 +1376,4 @@
 !--       fluxes, etc. 
           CALL init_surfaces
-!
-!--       Initialize synthetic turbulence generator if required 
-          IF( use_syn_turb_gen )  CALL stg_init
 
           CALL location_message( 'initializing with 1D model profiles', 'finished' )
@@ -1446,7 +1444,4 @@
 !--       fluxes, etc. 
           CALL init_surfaces
-!
-!--       Initialize synthetic turbulence generator if required 
-          IF( use_syn_turb_gen )  CALL stg_init
           
           CALL location_message( 'initializing with constant profiles', 'finished' )
@@ -1868,8 +1863,4 @@
        ENDIF
        IF ( passive_scalar )  ts_m  = 0.0_wp
-!
-!--    Initialize synthetic turbulence generator in case of restart.
-       IF ( TRIM( initializing_actions ) == 'read_restart_data'  .AND.         &
-            use_syn_turb_gen )  CALL stg_init
 
        CALL location_message( 'initializing in case of restart / cyclic_fill', 'finished' )
@@ -2315,12 +2306,6 @@
     IF ( agents_active )  CALL mas_init
 !
-!-- 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 )  THEN
-       CALL stg_adjust
-    ENDIF
+!-- Initialization of synthetic turbulence generator
+    IF ( use_syn_turb_gen )  CALL stg_init
 !
 !-- Initializing actions for all other modules
Index: palm/trunk/SOURCE/nesting_offl_mod.f90
===================================================================
--- palm/trunk/SOURCE/nesting_offl_mod.f90	(revision 3936)
+++ palm/trunk/SOURCE/nesting_offl_mod.f90	(revision 3937)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Set boundary conditon on upper-left and upper-south grid point for the u- and
+! v-component, respectively.
+! 
+! 3891 2019-04-12 17:52:01Z suehring
 ! Bugfix, do not overwrite lateral and top boundary data in case of restart 
 ! runs. 
@@ -633,4 +637,9 @@
           ENDDO
        ENDDO
+!
+!--    For left boundary set boundary condition for u-component also at top 
+!--    grid point. 
+!--    Note, this has no effect on the numeric solution, only for data output.
+       IF ( bc_dirichlet_l )  u(nzt+1,:,nxl) = u(nzt+1,:,nxlu)
 
        DO  i = nxl, nxr
@@ -644,4 +653,9 @@
           ENDDO
        ENDDO
+!
+!--    For south boundary set boundary condition for v-component also at top 
+!--    grid point. 
+!--    Note, this has no effect on the numeric solution, only for data output.
+       IF ( bc_dirichlet_s )  v(nzt+1,nys,:) = v(nzt+1,nysv,:)
 
        DO  i = nxl, nxr
@@ -961,7 +975,6 @@
        topo_max     = topo_max_l
 #endif
-
        zi_ribulk = MAX( zi_ribulk, topo_max )
-
+       
     END SUBROUTINE calc_zi
     
Index: palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90
===================================================================
--- palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90	(revision 3936)
+++ palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90	(revision 3937)
@@ -25,4 +25,10 @@
 ! -----------------
 ! $Id$
+! Minor bugfix in case of a very early restart where mc_factor is sill not 
+! present.
+! Some modification and fixing of potential bugs in the calculation of scaling 
+! parameters used for synthetic turbulence parametrization. 
+! 
+! 3909 2019-04-17 09:13:25Z suehring
 ! Minor bugfix for last commit
 ! 
@@ -270,5 +276,5 @@
     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) ::  mc_factor = 1.0_wp       !< mass flux correction factor
     REAL(wp) ::  scale_l                  !< scaling parameter used for turbulence parametrization - Obukhov length
     REAL(wp) ::  scale_us                 !< scaling parameter used for turbulence parametrization - friction velocity
@@ -886,5 +892,4 @@
           DO  j = nysg, nyng
              DO  k = nzb, nzt+1
-
                 IF  ( a11(k) .NE. 0._wp ) THEN
                    fu_yz(k,j) = ( u(k,j,i) / mc_factor - u_init(k) ) / a11(k)
@@ -1779,7 +1784,7 @@
 !--       zero. 
           r11(k) = scale_us**2 * (                                             &
-                      MERGE( 0.35_wp * (                                       &
+                      MERGE( 0.35_wp * ABS(                                    &
                            - zi_ribulk / ( kappa * scale_l - 10E-4_wp )        &
-                                       )**( 2.0_wp / 3.0_wp ),                 &
+                                          )**( 2.0_wp / 3.0_wp ),              &
                              0.0_wp,                                           &
                              scale_l < 0.0_wp )                                &
@@ -2052,9 +2057,7 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE calc_scaling_variables
-
-
-    USE control_parameters,                                                    &
-        ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, &
-               pt_surface
+             
+    USE arrays_3d,                                                          &
+        ONLY:  drho_air         
              
     USE indices,                                                               &
@@ -2074,4 +2077,7 @@
 
     REAL(wp) ::  friction_vel_l         !< mean friction veloctiy on subdomain
+    REAL(wp) ::  pt_surf_mean           !< mean near surface temperature (at 1st grid point)
+    REAL(wp) ::  pt_surf_mean_l         !< mean near surface temperature (at 1st grid point) on subdomain
+    REAL(wp) ::  scale_l_l              !< mean Obukhov lenght on subdomain
     REAL(wp) ::  shf_mean               !< mean surface sensible heat flux
     REAL(wp) ::  shf_mean_l             !< mean surface sensible heat flux on subdomain
@@ -2079,127 +2085,65 @@
     REAL(wp) ::  v_int                  !< v-component
     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
+!-- Calculate mean friction velocity, velocity scale, heat flux and 
+!-- near-surface temperature in the model domain.  
+    pt_surf_mean_l = 0.0_wp
+    shf_mean_l     = 0.0_wp
+    scale_l_l      = 0.0_wp
+    friction_vel_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)
+       i = surf_def_h(0)%i(m)
+       j = surf_def_h(0)%j(m)
+       k = surf_def_h(0)%k(m)
+       friction_vel_l = friction_vel_l  + surf_def_h(0)%us(m)
+       shf_mean_l     = shf_mean_l      + surf_def_h(0)%shf(m) * drho_air(k)
+       scale_l_l      = scale_l_l       + surf_def_h(0)%ol(m)
+       pt_surf_mean_l = pt_surf_mean_l  + pt(k,j,i)
     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)
+       i = surf_lsm_h%i(m)
+       j = surf_lsm_h%j(m)
+       k = surf_lsm_h%k(m)
+       friction_vel_l = friction_vel_l  + surf_lsm_h%us(m)
+       shf_mean_l     = shf_mean_l      + surf_lsm_h%shf(m) * drho_air(k)
+       scale_l_l      = scale_l_l       + surf_lsm_h%ol(m)
+       pt_surf_mean_l = pt_surf_mean_l  + pt(k,j,i)
     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)
+       i = surf_usm_h%i(m)
+       j = surf_usm_h%j(m)
+       k = surf_usm_h%k(m)
+       friction_vel_l = friction_vel_l  + surf_usm_h%us(m)
+       shf_mean_l     = shf_mean_l      + surf_usm_h%shf(m) * drho_air(k)
+       scale_l_l      = scale_l_l       + surf_usm_h%ol(m)
+       pt_surf_mean_l = pt_surf_mean_l  + pt(k,j,i)
     ENDDO
     
 #if defined( __parallel )
-    CALL MPI_ALLREDUCE( z0_mean_l, z0_mean, 1, MPI_REAL, MPI_SUM,              &
+    CALL MPI_ALLREDUCE( friction_vel_l, scale_us,     1, MPI_REAL, MPI_SUM,    &
                         comm2d, ierr )
-    CALL MPI_ALLREDUCE( shf_mean_l, shf_mean, 1, MPI_REAL, MPI_SUM,            &
+    CALL MPI_ALLREDUCE( shf_mean_l, shf_mean,         1, MPI_REAL, MPI_SUM,    &
+                        comm2d, ierr )
+    CALL MPI_ALLREDUCE( scale_l_l, scale_l,           1, MPI_REAL, MPI_SUM,    &
+                        comm2d, ierr )
+    CALL MPI_ALLREDUCE( pt_surf_mean_l, pt_surf_mean, 1, MPI_REAL, MPI_SUM,    &
                         comm2d, ierr )
 #else
-    z0_mean  = z0_mean_l
-    shf_mean = shf_mean_l
+    scale_us     = friction_vel_l
+    shf_mean     = shf_mean_l
+    scale_l      = scale_l_l
+    pt_surf_mean = pt_surf_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, 
-!-- which is 10 m, according to the height of the 1st vertical grid level 
-!-- in the COSMO level. However, in case of topography that is higher than
-!-- the reference level, the k index is determined from the 1st vertical 
-!-- PALM grid level instead. 
-!-- 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
-!
-!--       Determine the k index and topography top index
-          k_topo = MAX( get_topography_top_index_ji( j, i, 'u' ),              &
-                        get_topography_top_index_ji( j, i, 'v' ) )
-          k      = MAX( k_ref, k_topo + 1 )
-!
-!--       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,j,i+1), u(k,j,i), bc_dirichlet_l )
-          v_int = v(k,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) - zu(k_topo) ) / 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
-          
-          k_topo = MAX( get_topography_top_index_ji( j, i, 'u' ),              &
-                        get_topography_top_index_ji( j, i, 'v' ) )
-          k      = MAX( k_ref, k_topo + 1 )
-          
-          u_int = u(k,j,i)
-          v_int = MERGE( v(k,j+1,i), v(k,j,i), bc_dirichlet_s )
- 
-          friction_vel_l = friction_vel_l + kappa *                            &
-                            SQRT( u_int * u_int + v_int * v_int )  /           &
-                            LOG( ( zu(k) - zu(k_topo) ) / 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
-    
+
+    scale_us     = scale_us     / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
+    shf_mean     = shf_mean     / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
+    scale_l      = scale_l      / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
+    pt_surf_mean = pt_surf_mean / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )   
 !
 !-- 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 )
+       w_convective = ( g * shf_mean * zi_ribulk / pt_surf_mean )**( 1.0_wp / 3.0_wp )
     ELSE
        w_convective = 0.0_wp
