Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 4471)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 4472)
@@ -25,82 +25,85 @@
 ! -----------------
 ! $Id$
+! Kolmogorov length scale eta added to profile output
+!
+! 4444 2020-03-05 15:59:50Z raasch
 ! bugfix: cpp-directives for serial mode added
-! 
+!
 ! 4392 2020-01-31 16:14:57Z pavelkrc
 ! Some error numbers revised to prevent double usage
-! 
+!
 ! 11:55:33Z oliver.maas
 ! Checks for closed channel flow implemented
-! 
+!
 ! 11:55:33Z oliver.maas
 ! Move 2-m potential temperature output to diagnostic_output_quantities
-! 
+!
 ! 11:55:33Z oliver.maas
 ! removed message PA0421, concerning old parameter recycling_yshift
-! 
+!
 ! 11:55:33Z oliver.maas
 ! adjust message to the modified parameter recycling_yshift
-! 
+!
 ! 11:55:33Z oliver.maas
 ! Check if a cross section is specified if any output cross-section quantity
 ! is given
-! 
+!
 ! 11:55:33Z oliver.maas
 ! Overwrite rotation_angle from namelist by value from static driver
-! 
+!
 ! 11:55:33Z oliver.maas
-! removed conversion from recycle_absolute_quantities to raq, added check and 
+! removed conversion from recycle_absolute_quantities to raq, added check and
 ! error message for correct input of recycling_method_for_thermodynamic_quantities
-! 
+!
 ! 11:55:33Z oliver.maas
 ! Corrected "Former revisions" section
-! 
+!
 ! 11:55:33Z oliver.maas
 ! bugfix error message: replaced PA184 by PA0184
-! 
+!
 ! 11:55:33Z oliver.maas
-! added conversion from recycle_absolute_quantities to raq for recycling of 
+! added conversion from recycle_absolute_quantities to raq for recycling of
 ! absolute quantities and added error message PA184 for not implemented quantities
-! 
+!
 ! 4142 2019-08-05 12:38:31Z suehring
 ! Consider spinup in number of output timesteps for averaged 2D output (merge
 ! from branch resler).
-! 
+!
 ! 4069 2019-07-01 14:05:51Z Giersch
-! Masked output running index mid has been introduced as a local variable to 
+! Masked output running index mid has been introduced as a local variable to
 ! avoid runtime error (Loop variable has been modified) in time_integration
-! 
+!
 ! 4048 2019-06-21 21:00:21Z knoop
 ! Moved tcm_check_data_output to module_interface
-! 
+!
 ! 4039 2019-06-18 10:32:41Z suehring
 ! Modularize diagnostic output
-! 
+!
 ! 4017 2019-06-06 12:16:46Z schwenkel
 ! output of turbulence intensity added
-! 
+!
 ! 3933 2019-04-25 12:33:20Z kanani
 ! Alphabetical resorting in CASE, condense settings for theta_2m* into one IF clause
-! 
+!
 ! 3885 2019-04-11 11:29:34Z kanani
-! Changes related to global restructuring of location messages and introduction 
+! Changes related to global restructuring of location messages and introduction
 ! of additional debug messages
-! 
+!
 ! 3766 2019-02-26 16:23:41Z raasch
 ! trim added to avoid truncation compiler warnings
-! 
+!
 ! 3761 2019-02-25 15:31:42Z raasch
 ! unused variables removed
-! 
+!
 ! 3735 2019-02-12 09:52:40Z dom_dwd_user
-! Passing variable j (averaged output?) to 
+! Passing variable j (averaged output?) to
 ! module_interface.f90:chem_check_data_output.
-! 
+!
 ! 3705 2019-01-29 19:56:39Z suehring
 ! bugfix: renamed thetav_t to vtheta_t
-! 
+!
 ! 3702 2019-01-28 13:19:30Z gronemeier
 ! most_method removed
-! 
+!
 ! 3655 2019-01-07 16:51:22Z knoop
 ! Formatting
@@ -158,5 +161,5 @@
     USE particle_attributes,                                                   &
         ONLY:  particle_advection, use_sgs_for_particles
-        
+
     USE pegrid
 
@@ -242,5 +245,5 @@
        CALL message( 'check_parameters', 'PA0476', 1, 2, 0, 6, 0 )
     ENDIF
-   
+
 !
 !-- Check dt_coupling, restart_time, dt_restart, end_time, dx, dy, nx and ny
@@ -522,9 +525,9 @@
 !
 !--    Check illegal/untested parameter combinations for closed channel
-       If ( topography == 'closed_channel' ) THEN 
+       If ( topography == 'closed_channel' ) THEN
           symmetry_flag = 1
           message_string = 'Bottom and top boundary are treated equal'
           CALL message( 'check_parameters', 'PA0699', 0, 0, 0, 6, 0 )
-       
+
           IF ( dz(1) /= dz(COUNT( dz /= -1.0_wp )) .OR.                        &
                dz_stretch_level /= -9999999.9_wp) THEN
@@ -533,5 +536,5 @@
              CALL message( 'check_parameters', 'PA0700', 1, 2, 0, 6, 0 )
           ENDIF
-       
+
           IF ( constant_flux_layer ) THEN
              WRITE( message_string, * )  'A constant flux layer is not '//     &
@@ -540,5 +543,5 @@
              CALL message( 'check_parameters', 'PA0701', 1, 2, 0, 6, 0 )
           ENDIF
-       
+
           IF ( ocean_mode ) THEN
              WRITE( message_string, * )  'The ocean mode is not allowed if '// &
@@ -546,5 +549,5 @@
              CALL message( 'check_parameters', 'PA0702', 1, 2, 0, 6, 0 )
           ENDIF
-       
+
           IF ( momentum_advec /= 'ws-scheme' .OR.                              &
                scalar_advec /= 'ws-scheme' ) THEN
@@ -629,5 +632,5 @@
 
 !
-!-- When the land- or urban-surface model is used, the flux output must be 
+!-- When the land- or urban-surface model is used, the flux output must be
 !-- dynamic.
     IF ( land_surface  .OR.  urban_surface )  THEN
@@ -663,5 +666,5 @@
                                   ' its maximum of dots_max = ', dots_max,     &
                                   '&Please increase dots_max in modules.f90.'
-       CALL message( 'init_3d_model', 'PA0194', 1, 2, 0, 6, 0 )    
+       CALL message( 'init_3d_model', 'PA0194', 1, 2, 0, 6, 0 )
     ENDIF
 
@@ -705,5 +708,5 @@
 !
 !-- Advection schemes:
-    IF ( momentum_advec /= 'pw-scheme'  .AND.                                  &  
+    IF ( momentum_advec /= 'pw-scheme'  .AND.                                  &
          momentum_advec /= 'ws-scheme'  .AND.                                  &
          momentum_advec /= 'up-scheme' )                                       &
@@ -849,6 +852,6 @@
     ENDIF
 !
-!-- In case of spinup and nested run, spinup end time must be identical 
-!-- in order to have synchronously running simulations. 
+!-- In case of spinup and nested run, spinup end time must be identical
+!-- in order to have synchronously running simulations.
     IF ( nested_run )  THEN
 #if defined( __parallel )
@@ -1149,5 +1152,5 @@
 
 !
-!-- Overwrite parameters from namelist if necessary and compute Coriolis parameter. 
+!-- Overwrite parameters from namelist if necessary and compute Coriolis parameter.
 !-- @todo - move initialization of f and fs to coriolis_mod.
     IF ( input_pids_static )  THEN
@@ -2279,4 +2282,11 @@
              ENDIF
 
+           CASE ( 'eta' )
+              dopr_index(i) = 121
+              dopr_unit(i)  = 'mm'
+              hom(:,2,121,:) = SPREAD( zu, 2, statistic_regions+1 )
+
+              kolmogorov_length_scale = .TRUE.
+
           CASE DEFAULT
              unit = 'illegal'
@@ -2463,5 +2473,5 @@
 
              IF ( ( TRIM( var ) == 'r_a*' .OR.  TRIM( var ) == 'ghf*' )        &
-                 .AND.  .NOT.  land_surface  .AND.  .NOT.  urban_surface )     &         
+                 .AND.  .NOT.  land_surface  .AND.  .NOT.  urban_surface )     &
              THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
@@ -2470,5 +2480,5 @@
                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
              ENDIF
-             
+
              IF ( TRIM( var ) == 'ssws*'  .AND.  .NOT.  passive_scalar )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
@@ -2490,6 +2500,6 @@
              IF ( TRIM( var ) == 'z0h*'   )  unit = 'm'
 !
-!--          Output of surface latent and sensible heat flux will be in W/m2 
-!--          in case of natural- and urban-type surfaces, even if 
+!--          Output of surface latent and sensible heat flux will be in W/m2
+!--          in case of natural- and urban-type surfaces, even if
 !--          flux_output_mode is set to kinematic units.
              IF ( land_surface  .OR.  urban_surface )  THEN
@@ -2705,5 +2715,5 @@
        CALL check_dt_do( dt_data_output_av, 'dt_data_output_av' )
 
-!--    Set needed time levels (ntdim) to 
+!--    Set needed time levels (ntdim) to
 !--    saved time levels + to be saved time levels.
        ntdim_3d(0) = do3d_time_count(0) + CEILING(                                    &
@@ -2749,6 +2759,6 @@
        ENDIF
 !
-!--    Please note, for averaged 2D data skip_time_data_output_av is the relavant 
-!--    output control parameter. 
+!--    Please note, for averaged 2D data skip_time_data_output_av is the relavant
+!--    output control parameter.
        ntdim_2d_xy(1) = do2d_xy_time_count(1) + CEILING(                              &
                      ( end_time - MAX( MERGE( skip_time_data_output_av,               &
@@ -3158,5 +3168,5 @@
        CALL message( 'check_parameters', 'PA0157', 1, 2, 0, 6, 0 )
     ENDIF
-    
+
 !
 !-- Prevent empty time records in volume, cross-section and masked data in case
@@ -3197,5 +3207,5 @@
 
 !
-!-- Check if vertical grid stretching is switched off in case of complex 
+!-- Check if vertical grid stretching is switched off in case of complex
 !-- terrain simulations
     IF ( complex_terrain  .AND.                                                &
Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 4471)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 4472)
@@ -25,6 +25,9 @@
 ! -----------------
 ! $Id$
+! Calculations of the Kolmogorov lengt scale eta implemented
+!
+! 4464 2020-03-17 11:08:46Z Giersch
 ! Reset last change (r4463)
-! 
+!
 ! 4463 2020-03-17 09:27:36Z Giersch
 ! Calculate horizontally averaged profiles of all velocity components at the
@@ -33,29 +36,29 @@
 ! 4444 2020-03-05 15:59:50Z raasch
 ! bugfix: cpp-directives for serial mode added
-! 
+!
 ! 4442 2020-03-04 19:21:13Z suehring
-! Change order of dimension in surface array %frac to allow for better 
+! Change order of dimension in surface array %frac to allow for better
 ! vectorization.
-! 
+!
 ! 4441 2020-03-04 19:20:35Z suehring
 ! Introduction of wall_flags_total_0, which currently sets bits based on static
 ! topography information used in wall_flags_static_0
-! 
+!
 ! 4329 2019-12-10 15:46:36Z motisi
 ! Renamed wall_flags_0 to wall_flags_static_0
-! 
+!
 ! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
-! 
+!
 ! 4131 2019-08-02 11:06:18Z monakurppa
 ! Allow profile output for salsa variables.
-! 
+!
 ! 4039 2019-06-18 10:32:41Z suehring
-! Correct conversion to kinematic scalar fluxes in case of pw-scheme and 
+! Correct conversion to kinematic scalar fluxes in case of pw-scheme and
 ! statistic regions
-! 
+!
 ! 3828 2019-03-27 19:36:23Z raasch
 ! unused variables removed
-! 
+!
 ! 3676 2019-01-16 15:07:05Z knoop
 ! Bugfix, terminate OMP Parallel block
@@ -71,7 +74,7 @@
 !>
 !> @note For simplicity, nzb_s_inner and nzb_diff_s_inner are being used as a
-!>       lower vertical index for k-loops for all variables, although strictly 
-!>       speaking the k-loops would have to be split up according to the staggered 
-!>       grid. However, this implies no error since staggered velocity components 
+!>       lower vertical index for k-loops for all variables, although strictly
+!>       speaking the k-loops would have to be split up according to the staggered
+!>       grid. However, this implies no error since staggered velocity components
 !>       are zero at the walls and inside buildings.
 !------------------------------------------------------------------------------!
@@ -97,8 +100,8 @@
     USE control_parameters,                                                    &
         ONLY:   air_chemistry, average_count_pr, cloud_droplets, do_sum,       &
-                dt_3d, humidity, initializing_actions, land_surface,           &
-                large_scale_forcing, large_scale_subsidence, max_pr_user,      &
-                message_string, neutral, ocean_mode, passive_scalar,           &
-                simulated_time, simulated_time_at_begin,                       &
+                dt_3d, humidity, initializing_actions, kolmogorov_length_scale,&
+                land_surface, large_scale_forcing, large_scale_subsidence,     &
+                max_pr_user, message_string, neutral, ocean_mode,              &
+                passive_scalar, simulated_time, simulated_time_at_begin,       &
                 use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
                 ws_scheme_mom, ws_scheme_sca, salsa, max_pr_salsa
@@ -109,5 +112,5 @@
     USE grid_variables,                                                        &
         ONLY:   ddx, ddy
-        
+
     USE indices,                                                               &
         ONLY:   ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, nxl, nxr, nyn, &
@@ -118,7 +121,7 @@
         ONLY:  ngp_sums, ngp_sums_ls
 #endif
-        
+
     USE kinds
-    
+
     USE land_surface_model_mod,                                                &
         ONLY:   m_soil_h, nzb_soil, nzt_soil, t_soil_h
@@ -151,7 +154,7 @@
     INTEGER(iwp) ::  j                   !<
     INTEGER(iwp) ::  k                   !<
-    INTEGER(iwp) ::  ki                  !< 
+    INTEGER(iwp) ::  ki                  !<
     INTEGER(iwp) ::  k_surface_level     !<
-    INTEGER(iwp) ::  m                   !< loop variable over all horizontal wall elements 
+    INTEGER(iwp) ::  m                   !< loop variable over all horizontal wall elements
     INTEGER(iwp) ::  l                   !< loop variable over surface facing -- up- or downward-facing
     INTEGER(iwp) ::  nt                  !<
@@ -161,10 +164,30 @@
 
     LOGICAL ::  first  !<
-    
-    REAL(wp) ::  dptdz_threshold  !< 
+
+    REAL(wp) ::  dissipation      !< dissipation rate
+    REAL(wp) ::  dptdz_threshold  !<
+    REAL(wp) ::  du_dx            !< Derivative of u fluctuations with respect to x
+    REAL(wp) ::  du_dy            !< Derivative of u fluctuations with respect to y
+    REAL(wp) ::  du_dz            !< Derivative of u fluctuations with respect to z
+    REAL(wp) ::  dv_dx            !< Derivative of v fluctuations with respect to x
+    REAL(wp) ::  dv_dy            !< Derivative of v fluctuations with respect to y
+    REAL(wp) ::  dv_dz            !< Derivative of v fluctuations with respect to z
+    REAL(wp) ::  dw_dx            !< Derivative of w fluctuations with respect to x
+    REAL(wp) ::  dw_dy            !< Derivative of w fluctuations with respect to y
+    REAL(wp) ::  dw_dz            !< Derivative of w fluctuations with respect to z
+    REAL(wp) ::  eta              !< Kolmogorov length scale
     REAL(wp) ::  fac              !<
     REAL(wp) ::  flag             !<
     REAL(wp) ::  height           !<
     REAL(wp) ::  pts              !<
+    REAL(wp) ::  s11              !< fluctuating rate-of-strain tensor component 11
+    REAL(wp) ::  s21              !< fluctuating rate-of-strain tensor component 21
+    REAL(wp) ::  s31              !< fluctuating rate-of-strain tensor component 31
+    REAL(wp) ::  s12              !< fluctuating rate-of-strain tensor component 12
+    REAL(wp) ::  s22              !< fluctuating rate-of-strain tensor component 22
+    REAL(wp) ::  s32              !< fluctuating rate-of-strain tensor component 32
+    REAL(wp) ::  s13              !< fluctuating rate-of-strain tensor component 13
+    REAL(wp) ::  s23              !< fluctuating rate-of-strain tensor component 23
+    REAL(wp) ::  s33              !< fluctuating rate-of-strain tensor component 33
     REAL(wp) ::  sums_l_etot      !<
     REAL(wp) ::  ust              !<
@@ -175,5 +198,5 @@
     REAL(wp) ::  v2               !<
     REAL(wp) ::  w2               !<
-    
+
     REAL(wp) ::  dptdz(nzb+1:nzt+1)    !<
     REAL(wp) ::  sums_ll(nzb:nzt+1,2)  !<
@@ -210,5 +233,5 @@
 !--    array
        sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
-!--    WARNING: next line still has to be adjusted for OpenMP 
+!--    WARNING: next line still has to be adjusted for OpenMP
        sums_l(:,21,0) = sums_wsts_bc_l(:,sr) *                                 &
                         heatflux_output_conversion  ! heat flux from advec_s_bc
@@ -220,14 +243,14 @@
 !--    scale) vertical fluxes and velocity variances by using commonly-
 !--    applied Reynolds-based methods ( e.g. <w'pt'> = (w-<w>)*(pt-<pt>) )
-!--    in combination with the 5th order advection scheme, pronounced 
-!--    artificial kinks could be observed in the vertical profiles near the 
-!--    surface. Please note: these kinks were not related to the model truth, 
-!--    i.e. these kinks are just related to an evaluation problem.   
-!--    In order avoid these kinks, vertical fluxes and horizontal as well 
+!--    in combination with the 5th order advection scheme, pronounced
+!--    artificial kinks could be observed in the vertical profiles near the
+!--    surface. Please note: these kinks were not related to the model truth,
+!--    i.e. these kinks are just related to an evaluation problem.
+!--    In order avoid these kinks, vertical fluxes and horizontal as well
 !--    vertical velocity variances are calculated directly within the advection
-!--    routines, according to the numerical discretization, to evaluate the 
-!--    statistical quantities as they will appear within the prognostic 
+!--    routines, according to the numerical discretization, to evaluate the
+!--    statistical quantities as they will appear within the prognostic
 !--    equations.
-!--    Copy the turbulent quantities, evaluated in the advection routines to 
+!--    Copy the turbulent quantities, evaluated in the advection routines to
 !--    the local array sums_l() for further computations.
        IF ( ws_scheme_mom .AND. sr == 0 )  THEN
@@ -248,8 +271,8 @@
              sums_l(:,15,i) = sums_wsvs_ws_l(:,i)                              &
                               * momentumflux_output_conversion ! w*v*
-             sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2 
-             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2 
-             sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2 
-             sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp *                        & 
+             sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
+             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
+             sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2
+             sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp *                        &
                               ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +      &
                                 sums_ws2_ws_l(:,i) )    ! e*
@@ -270,5 +293,5 @@
 
        ENDIF
-! 
+!
 !--    Horizontally averaged profiles of horizontal velocities and temperature.
 !--    They must have been computed before, because they are already required
@@ -446,5 +469,5 @@
 
 !
-!--    Final values are obtained by division by the total number of grid points 
+!--    Final values are obtained by division by the total number of grid points
 !--    used for summation. After that store profiles.
        sums(:,1) = sums(:,1) / ngp_2dh(sr)
@@ -482,5 +505,5 @@
 !--    Passive scalar
        IF ( passive_scalar )  hom(:,1,115,sr) = sums(:,115) /                  &
-            ngp_2dh_s_inner(:,sr)                    ! s 
+            ngp_2dh_s_inner(:,sr)                    ! s
 
 !
@@ -488,6 +511,6 @@
 !--    variances, the total and the perturbation energy (single values in last
 !--    column of sums_l) and some diagnostic quantities.
-!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly 
-!--    ----  speaking the following k-loop would have to be split up and 
+!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
+!--    ----  speaking the following k-loop would have to be split up and
 !--          rearranged according to the staggered grid.
 !--          However, this implies no error since staggered velocity components
@@ -500,8 +523,12 @@
        !$OMP DO
        !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, m) &
-       !$ACC PRIVATE(sums_l_etot, flag) &
+       !$ACC PRIVATE(sums_l_etot, flag, du_dx, du_dy, du_dz) &
+       !$ACC PRIVATE(dv_dx, dv_dy, dv_dz, dw_dx, dw_dy, dw_dz) &
+       !$ACC PRIVATE(s11, s21, s31, s12, s22, s32, s13, s23, s33) &
+       !$ACC PRIVATE(dissipation, eta) &
        !$ACC PRESENT(wall_flags_total_0, rmask, momentumflux_output_conversion) &
-       !$ACC PRESENT(hom(:,1,4,sr)) &
+       !$ACC PRESENT(hom(:,1,1:2,sr), hom(:,1,4,sr)) &
        !$ACC PRESENT(e, u, v, w, km, kh, p, pt) &
+       !$ACC PRESENT(ddx, ddy, ddzu, ddzw) &
        !$ACC PRESENT(surf_def_h(0), surf_lsm_h, surf_usm_h) &
        !$ACC PRESENT(sums_l)
@@ -557,5 +584,91 @@
                                         w(k,j,i)**2 )            * rmask(j,i,sr)&
                                                                  * flag
-             ENDDO
+
+!
+!--             Computation of the Kolmogorov length scale. Calculation is based
+!--             on gradients of the deviations from the horizontal mean.
+!--             Kolmogorov scale at the boundaries (k=0/z=0m and k=nzt+1) is set to zero.
+                IF ( kolmogorov_length_scale .AND. k /= nzb .AND. k /= nzt+1) THEN
+                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 22 ) )
+
+!
+!--                Calculate components of the fluctuating rate-of-strain tensor
+!--                (0.5*(del u'_i/del x_j + del u'_j/del x_i)) defined in the
+!--                center of each grid box.
+                   du_dx = ( ( u(k,j,i+1) - hom(k,1,1,sr) ) -                  &
+                             ( u(k,j,i) - hom(k,1,1,sr) ) ) * ddx
+                   du_dy = 0.25_wp * ddy *                                     &
+                           ( ( u(k,j+1,i) - hom(k,1,1,sr) ) -                  &
+                             ( u(k,j-1,i) - hom(k,1,1,sr) ) +                  &
+                             ( u(k,j+1,i+1) - hom(k,1,1,sr) ) -                &
+                             ( u(k,j-1,i+1) - hom(k,1,1,sr) ) )
+                   du_dz = 0.25_wp * ( ( ( u(k+1,j,i) - hom(k+1,1,1,sr) ) -    &
+                                         ( u(k,j,i) - hom(k,1,1,sr) ) ) *      &
+                                        ddzu(k+1) +                            &
+                                       ( ( u(k,j,i) - hom(k,1,1,sr) ) -        &
+                                         ( u(k-1,j,i) - hom(k-1,1,1,sr) ) )*   &
+                                        ddzu(k) +                              &
+                                       ( ( u(k+1,j,i+1) - hom(k+1,1,1,sr) )-   &
+                                         ( u(k,j,i+1) - hom(k,1,1,sr) ) ) *    &
+                                        ddzu(k+1) +                            &
+                                       ( ( u(k,j,i+1) - hom(k,1,1,sr) ) -      &
+                                         ( u(k-1,j,i+1) - hom(k-1,1,1,sr) ) ) *&
+                                        ddzu(k) )
+
+                   dv_dx = 0.25_wp * ddx *                                     &
+                           ( ( v(k,j,i+1) - hom(k,1,2,sr) ) -                  &
+                             ( v(k,j,i-1) - hom(k,1,2,sr) ) +                  &
+                             ( v(k,j+1,i+1) - hom(k,1,2,sr) ) -                &
+                             ( v(k,j+1,i-1) - hom(k,1,2,sr) ) )
+                   dv_dy = ( ( v(k,j+1,i) - hom(k,1,2,sr) ) -                  &
+                             ( v(k,j,i) - hom(k,1,2,sr) ) ) * ddy
+                   dv_dz = 0.25_wp * ( ( ( v(k+1,j,i) - hom(k+1,1,2,sr) ) -    &
+                                         ( v(k,j,i) - hom(k,1,2,sr) ) ) *      &
+                                        ddzu(k+1) +                            &
+                                       ( ( v(k,j,i) - hom(k,1,2,sr) ) -        &
+                                         ( v(k-1,j,i) - hom(k-1,1,2,sr) ) ) *  &
+                                        ddzu(k) +                              &
+                                       ( ( v(k+1,j+1,i) - hom(k+1,1,2,sr) ) -  &
+                                         ( v(k,j+1,i) - hom(k,1,2,sr) ) ) *    &
+                                        ddzu(k+1) +                            &
+                                       ( ( v(k,j+1,i) - hom(k,1,2,sr) ) -      &
+                                         ( v(k-1,j+1,i) - hom(k-1,1,2,sr) ) ) *&
+                                        ddzu(k) )
+
+                   dw_dx = 0.25_wp * ddx * ( w(k,j,i+1) - w(k,j,i-1) +         &
+                                             w(k-1,j,i+1) - w(k-1,j,i-1) )
+                   dw_dy = 0.25_wp * ddy * ( w(k,j+1,i) - w(k,j-1,i) +         &
+                                             w(k-1,j+1,i) - w(k-1,j-1,i) )
+                   dw_dz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
+
+                   s11 = 0.5_wp * ( du_dx + du_dx )
+                   s21 = 0.5_wp * ( dv_dx + du_dy )
+                   s31 = 0.5_wp * ( dw_dx + du_dz )
+
+                   s12 = s21
+                   s22 = 0.5 * ( dv_dy + dv_dy )
+                   s32 = 0.5 * ( dw_dy + dv_dz )
+
+                   s13 = s31
+                   s23 = s32
+                   s33 = 0.5_wp * ( dw_dz + dw_dz )
+
+!--                Calculate 3D instantaneous energy dissipation rate after
+!--                Pope (2000): Turbulent flows, p.259. It is defined in the center
+!--                of each grid volume.
+                   dissipation = 2.0_wp * km(k,j,i) *                          &
+                                ( s11*s11 + s21*s21 + s31*s31 +                &
+                                 s12*s12 + s22*s22 + s32*s32 +                 &
+                                s13*s13 + s23*s23 + s33*s33 )
+                   eta         = ( km(k,j,i)**3.0_wp / ( dissipation+1.0E-12 ) )**(1.0_wp/4.0_wp)
+
+                   !$ACC ATOMIC
+                   sums_l(k,121,tn) = sums_l(k,121,tn) + eta * rmask(j,i,sr)   &
+                                                                * flag
+
+
+                ENDIF !Kolmogorov length scale
+
+             ENDDO !k-loop
 !
 !--          Total and perturbation energy for the total domain (being
@@ -661,14 +774,14 @@
                                             rmask(j,i,sr)
              ENDIF
-          ENDDO
-       ENDDO
+          ENDDO !j-loop
+       ENDDO !i-loop
        !$ACC UPDATE &
        !$ACC HOST(sums_l(:,3,tn), sums_l(:,8,tn), sums_l(:,9,tn)) &
        !$ACC HOST(sums_l(:,10,tn), sums_l(:,40,tn), sums_l(:,33,tn)) &
-       !$ACC HOST(sums_l(:,38,tn)) &
+       !$ACC HOST(sums_l(:,38,tn), sums_l(:,121,tn)) &
        !$ACC HOST(sums_l(nzb:nzb+4,pr_palm,tn), sums_l(nzb+14:nzb+14,pr_palm,tn))
 
 !
-!--    Computation of statistics when ws-scheme is not used. Else these 
+!--    Computation of statistics when ws-scheme is not used. Else these
 !--    quantities are evaluated in the advection routines.
        IF ( .NOT. ws_scheme_mom .OR. sr /= 0 .OR. simulated_time == 0.0_wp )   &
@@ -703,9 +816,9 @@
        ENDIF
 !
-!--    Computaion of domain-averaged perturbation energy. Please note, 
-!--    to prevent that perturbation energy is larger (even if only slightly) 
+!--    Computaion of domain-averaged perturbation energy. Please note,
+!--    to prevent that perturbation energy is larger (even if only slightly)
 !--    than the total kinetic energy, calculation is based on deviations from
 !--    the horizontal mean, instead of spatial descretization of the advection
-!--    term. 
+!--    term.
        !$OMP DO
        !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k, flag, w2, ust2, vst2) &
@@ -727,4 +840,5 @@
                                  * rmask(j,i,sr)                               &
                                  * flag
+
              ENDDO
           ENDDO
@@ -746,13 +860,13 @@
           DO  j = nys, nyn
 !
-!--          Subgridscale fluxes (without Prandtl layer from k=nzb, 
+!--          Subgridscale fluxes (without Prandtl layer from k=nzb,
 !--          oterwise from k=nzb+1)
 !--          NOTE: for simplicity, nzb_diff_s_inner is used below, although
 !--          ----  strictly speaking the following k-loop would have to be
 !--                split up according to the staggered grid.
-!--                However, this implies no error since staggered velocity 
+!--                However, this implies no error since staggered velocity
 !--                components are zero at the walls and inside buildings.
 !--          Flag 23 is used to mask surface fluxes as well as model-top fluxes,
-!--          which are added further below. 
+!--          which are added further below.
              DO  k = nzb, nzt
                 flag = MERGE( 1.0_wp, 0.0_wp,                                  &
@@ -1184,5 +1298,5 @@
                                     0.0_wp * rmask(j,i,sr)        ! u"pt"
                 sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + &
-                                    0.0_wp * rmask(j,i,sr)        ! v"pt" 
+                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
 #endif
 #ifndef _OPENACC
@@ -1227,6 +1341,6 @@
 !
 !--          Resolved fluxes (can be computed for all horizontal points)
-!--          NOTE: for simplicity, nzb_s_inner is used below, although strictly 
-!--          ----  speaking the following k-loop would have to be split up and 
+!--          NOTE: for simplicity, nzb_s_inner is used below, although strictly
+!--          ----  speaking the following k-loop would have to be split up and
 !--                rearranged according to the staggered grid.
              DO  k = nzb, nzt
@@ -1303,5 +1417,5 @@
                             sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) *    &
                                                                 rmask(j,i,sr) *&
-                                                                flag 
+                                                                flag
                             sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) *    &
                                                                 rmask(j,i,sr) *&
@@ -1320,5 +1434,5 @@
                                                              flag
                       ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
-                         sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp *              & 
+                         sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp *              &
                                                hom(k,1,41,sr) ) *              &
                                              sums_l(k,17,tn) +                 &
@@ -1359,5 +1473,5 @@
        !$ACC HOST(sums_l(:,35,tn), sums_l(:,36,tn), sums_l(:,37,tn))
 !
-!--    Treat land-surface quantities according to new wall model structure. 
+!--    Treat land-surface quantities according to new wall model structure.
        IF ( land_surface )  THEN
           tn = 0
@@ -1368,7 +1482,7 @@
              i = surf_lsm_h%i(m)
              j = surf_lsm_h%j(m)
-        
+
              IF ( i >= nxl  .AND.  i <= nxr  .AND.                             &
-                  j >= nys  .AND.  j <= nyn )  THEN 
+                  j >= nys  .AND.  j <= nyn )  THEN
                 sums_l(nzb,93,tn)  = sums_l(nzb,93,tn) + surf_lsm_h%ghf(m)
                 sums_l(nzb,94,tn)  = sums_l(nzb,94,tn) + surf_lsm_h%qsws_liq(m)
@@ -1387,9 +1501,9 @@
           DO  m = 1, surf_lsm_h%ns
 
-             i = surf_lsm_h%i(m)            
+             i = surf_lsm_h%i(m)
              j = surf_lsm_h%j(m)
 
              IF ( i >= nxl  .AND.  i <= nxr  .AND.                             &
-                  j >= nys  .AND.  j <= nyn )  THEN 
+                  j >= nys  .AND.  j <= nyn )  THEN
 
                 DO  k = nzb_soil, nzt_soil
@@ -1417,6 +1531,6 @@
                 DO  k = nzb, nzt
 !
-!--                Flag 23 is used to mask surface fluxes as well as model-top 
-!--                fluxes, which are added further below. 
+!--                Flag 23 is used to mask surface fluxes as well as model-top
+!--                fluxes, which are added further below.
                    flag = MERGE( 1.0_wp, 0.0_wp,                               &
                                  BTEST( wall_flags_total_0(k,j,i), 23 ) ) *    &
@@ -1569,6 +1683,6 @@
 
 !
-!--    Horizontal heat fluxes (subgrid, resolved, total). 
-!--    Do it only, if profiles shall be plotted. 
+!--    Horizontal heat fluxes (subgrid, resolved, total).
+!--    Do it only, if profiles shall be plotted.
        IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN
 
@@ -1628,5 +1742,5 @@
        IF ( large_scale_forcing  .AND.  ( simulated_time > 0.0_wp ) )  THEN
 !
-!--       Interpolation in time of LSF_DATA 
+!--       Interpolation in time of LSF_DATA
           nt = 1
           DO WHILE ( simulated_time - dt_3d > time_vert(nt) )
@@ -1669,5 +1783,5 @@
        tn = 0
        !$OMP PARALLEL PRIVATE( i, j, k, tn )
-       !$ tn = omp_get_thread_num()       
+       !$ tn = omp_get_thread_num()
        IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
           !$OMP DO
@@ -1718,5 +1832,5 @@
 
              IF ( air_chemistry )  THEN
-                IF ( max_pr_cs > 0 )  THEN                                 
+                IF ( max_pr_cs > 0 )  THEN
                      sums_l(:,pr_palm+max_pr_user+1:pr_palm + max_pr_user+ max_pr_cs,0) =          &
                                sums_l(:,pr_palm+max_pr_user+1:pr_palm + max_pr_user+max_pr_cs,0) + &
@@ -1774,5 +1888,5 @@
 
 !
-!--    Final values are obtained by division by the total number of grid points 
+!--    Final values are obtained by division by the total number of grid points
 !--    used for summation. After that store profiles.
 !--    Check, if statistical regions do contain at least one grid point at the
@@ -1834,5 +1948,5 @@
 
        IF ( air_chemistry ) THEN
-          IF ( max_pr_cs > 0 )  THEN                  
+          IF ( max_pr_cs > 0 )  THEN
              DO k = nzb, nzt+1
                 sums(k, pr_palm+1:pr_palm+max_pr_user+max_pr_cs) = &
@@ -1840,5 +1954,5 @@
                                  ngp_2dh_s_inner(k,sr)
              ENDDO
-          ENDIF  
+          ENDIF
        ENDIF
 
@@ -1850,5 +1964,5 @@
                   / ngp_2dh_s_inner(k,sr)
              ENDDO
-          ENDIF  
+          ENDIF
        ENDIF
 
@@ -1873,5 +1987,5 @@
        hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
                                        ! profile 24 is initial profile (sa)
-                                       ! profiles 25-29 left empty for initial 
+                                       ! profiles 25-29 left empty for initial
                                        ! profiles
        hom(:,1,30,sr) = sums(:,30)     ! u*2
@@ -1887,5 +2001,5 @@
        hom(:,1,40,sr) = sums(:,40)     ! p
        hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
-       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*       
+       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*
        hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
        hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
@@ -1893,5 +2007,5 @@
        hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
        hom(:,1,51,sr) = sums(:,51)     ! w"qv"
-       hom(:,1,52,sr) = sums(:,52)     ! w*qv*       
+       hom(:,1,52,sr) = sums(:,52)     ! w*qv*
        hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
        hom(:,1,54,sr) = sums(:,54)     ! ql
@@ -1979,5 +2093,5 @@
           hom(:,1,117,sr) = sums(:,117)     ! w"s"
           hom(:,1,114,sr) = sums(:,114)     ! w*s*
-          hom(:,1,118,sr) = sums(:,117) + sums(:,114)    ! ws 
+          hom(:,1,118,sr) = sums(:,117) + sums(:,114)    ! ws
           hom(:,1,116,sr) = sums(:,116)     ! s*2
        ENDIF
@@ -1985,4 +2099,9 @@
        hom(:,1,119,sr) = rho_air       ! rho_air in Kg/m^3
        hom(:,1,120,sr) = rho_air_zw    ! rho_air_zw in Kg/m^3
+
+       IF ( kolmogorov_length_scale ) THEN
+          hom(:,1,121,sr) = sums(:,121) * 1E3_wp  ! eta in mm
+       ENDIF
+
 
        hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
@@ -1995,5 +2114,5 @@
 
        IF ( air_chemistry )  THEN
-          IF ( max_pr_cs > 0 )  THEN    ! chem_spcs profiles      
+          IF ( max_pr_cs > 0 )  THEN    ! chem_spcs profiles
              hom(:, 1, pr_palm+max_pr_user+1:pr_palm + max_pr_user+max_pr_cs, sr) = &
                                sums(:, pr_palm+max_pr_user+1:pr_palm+max_pr_user+max_pr_cs)
@@ -2013,9 +2132,9 @@
 !--    The corresponding height is assumed as the boundary layer height, if it
 !--    is less than 1.5 times the height where the heat flux becomes negative
-!--    (positive) for the first time. Attention: the resolved vertical sensible 
+!--    (positive) for the first time. Attention: the resolved vertical sensible
 !--    heat flux (hom(:,1,17,sr) = w*pt*) is not known at the beginning because
-!--    the calculation happens in advec_s_ws which is called after 
-!--    flow_statistics. Therefore z_i is directly taken from restart data at 
-!--    the beginning of restart runs.  
+!--    the calculation happens in advec_s_ws which is called after
+!--    flow_statistics. Therefore z_i is directly taken from restart data at
+!--    the beginning of restart runs.
        IF ( TRIM( initializing_actions ) /= 'read_restart_data' .OR.           &
             simulated_time_at_begin /= simulated_time ) THEN
@@ -2059,11 +2178,11 @@
 
 !
-!--       Second scheme: Gradient scheme from Sullivan et al. (1998), modified 
-!--       by Uhlenbrock(2006). The boundary layer height is the height with the 
+!--       Second scheme: Gradient scheme from Sullivan et al. (1998), modified
+!--       by Uhlenbrock(2006). The boundary layer height is the height with the
 !--       maximal local temperature gradient: starting from the second (the
-!--       last but one) vertical gridpoint, the local gradient must be at least 
+!--       last but one) vertical gridpoint, the local gradient must be at least
 !--       0.2K/100m and greater than the next four gradients.
 !--       WARNING: The threshold value of 0.2K/100m must be adjusted for the
-!--       ocean case! 
+!--       ocean case!
           z_i(2) = 0.0_wp
           DO  k = nzb+1, nzt+1
@@ -2126,9 +2245,9 @@
 
 !
-!--    Collect the time series quantities. Please note, timeseries quantities 
-!--    which are collected from horizontally averaged profiles, e.g. wpt 
-!--    or pt(zp), are treated specially. In case of elevated model surfaces, 
+!--    Collect the time series quantities. Please note, timeseries quantities
+!--    which are collected from horizontally averaged profiles, e.g. wpt
+!--    or pt(zp), are treated specially. In case of elevated model surfaces,
 !--    index nzb+1 might be within topography and data will be zero. Therefore,
-!--    take value for the first atmosphere index, which is topo_min_level+1. 
+!--    take value for the first atmosphere index, which is topo_min_level+1.
        ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)        ! E
        ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)        ! E*
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 4471)
+++ palm/trunk/SOURCE/modules.f90	(revision 4472)
@@ -25,91 +25,95 @@
 ! -----------------
 ! $Id$
+! Additional switch added to activate calculations in flow_statistics for the
+! kolmogorov length scale
+!
+! 4461 2020-03-12 16:51:59Z raasch
 ! +virtual_pe_grid, communicator_configurations
-! 
+!
 ! 4414 2020-02-19 20:16:04Z suehring
-! - nzb_diff_s_inner, nzb_diff_s_outer, nzb_inner,nzb_outer, nzb_s_inner, 
-!   nzb_s_outer, nzb_u_inner, nzb_u_outer, nzb_v_inner, nzb_v_outer, 
+! - nzb_diff_s_inner, nzb_diff_s_outer, nzb_inner,nzb_outer, nzb_s_inner,
+!   nzb_s_outer, nzb_u_inner, nzb_u_outer, nzb_v_inner, nzb_v_outer,
 !   nzb_w_inner, nzb_w_outer
 !
-! 
+!
 ! 4360 2020-01-07 11:25:50Z suehring
 ! Introduction of wall_flags_total_0, which currently sets bits based on static
 ! topography information used in wall_flags_static_0
-! 
+!
 ! 4340 2019-12-16 08:17:03Z Giersch
 ! Flag for topography closed channel flow with symmetric boundaries introduced
-! 
+!
 ! 4331 2019-12-10 18:25:02Z suehring
 ! - do_output_at_2m, pt_2m_av
-! 
+!
 ! 4329 2019-12-10 15:46:36Z motisi
 ! Renamed wall_flags_0 to wall_flags_static_0
-! 
+!
 ! 4301 2019-11-22 12:09:09Z oliver.maas
 ! removed recycling_yshift
-! 
+!
 ! 4297 2019-11-21 10:37:50Z oliver.maas
 ! changed variable type of recycling_yshift from LOGICAL to INTEGER
-! 
+!
 ! 4293 2019-11-12 14:44:01Z Giersch
 ! Add origin_date_time
-! 
+!
 ! 4146 2019-08-07 07:47:36Z gronemeier
 ! Added rotation_angle
-! 
+!
 ! 4184 2019-08-23 08:07:40Z oliver.maas
-! changed allocated length of recycling_method_for_thermodynamic_quantities 
+! changed allocated length of recycling_method_for_thermodynamic_quantities
 ! from 20 to 80 characters
-! 
+!
 ! 4183 2019-08-23 07:33:16Z oliver.maas
 ! removed recycle_absolute_quantities and raq
 ! added recycling_method_for_thermodynamic_quantities
-! 
+!
 ! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
-! 
+!
 ! 4173 2019-08-20 12:04:06Z gronemeier
 ! add vdi_internal_controls
-! 
+!
 ! 4172 2019-08-20 11:55:33Z oliver.maas
 ! added recycle_absolute_quantities and raq
-! 
+!
 ! 4168 2019-08-16 13:50:17Z suehring
 ! +topo_top_ind
-! 
+!
 ! 4131 2019-08-02 11:06:18Z monakurppa
 ! Add max_pr_salsa to control_parameters. Used in creating profile output for
 ! salsa.
-! 
+!
 ! 4110 2019-07-22 17:05:21Z suehring
 ! -advc_flags_1, advc_flags_2
 ! +advc_flags_m, advc_flags_s
-! 
+!
 ! 4109 2019-07-22 17:00:34Z suehring
 ! remove old_dt
-! 
+!
 ! 4079 2019-07-09 18:04:41Z suehring
 ! + monotonic_limiter_z
-! 
+!
 ! 4069 2019-07-01 14:05:51Z Giersch
-! Masked output running index mid has been introduced as a local variable to 
+! Masked output running index mid has been introduced as a local variable to
 ! avoid runtime error (Loop variable has been modified) in time_integration
-! 
+!
 ! 4017 2019-06-06 12:16:46Z schwenkel
 ! increase maximum number of virtual flights
-! 
+!
 ! 3987 2019-05-22 09:52:13Z kanani
 ! Introduce alternative switch for debug output during timestepping
-! 
+!
 ! 3885 2019-04-11 11:29:34Z kanani
-! Changes related to global restructuring of location messages and introduction 
+! Changes related to global restructuring of location messages and introduction
 ! of additional debug messages
-! 
+!
 ! 3871 2019-04-08 14:38:39Z knoop
 ! Initialized parameter region
-! 
+!
 ! 3746 2019-02-16 12:41:27Z gronemeier
 ! Removed most_method
-! 
+!
 ! 3648 2019-01-02 16:35:46Z suehring
 ! -surface_data_output +surface_output
@@ -128,5 +132,5 @@
 !------------------------------------------------------------------------------!
  MODULE advection
- 
+
     USE kinds
 
@@ -135,5 +139,5 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  dex  !< exponential coefficient for the Bott-Chlond advection scheme
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  eex  !< exponential coefficient for the Bott-Chlond advection scheme
-    
+
     SAVE
 
@@ -150,5 +154,5 @@
 !------------------------------------------------------------------------------!
  MODULE mas_global_attributes
- 
+
     USE kinds
 
@@ -175,20 +179,20 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  c_w_m                  !< mean phase velocity at outflow for w-component used in radiation boundary condition
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  c_w_m_l                !< mean phase velocity at outflow for w-component used in radiation boundary condition (local subdomain value)
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ddzu                   !< 1/dzu 
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ddzu_pres              !< modified ddzu for pressure solver 
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  dd2zu                  !< 1/(dzu(k)+dzu(k+1)) 
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  dzu                    !< vertical grid size (u-grid) 
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ddzw                   !< 1/dzw 
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  dzw                    !< vertical grid size (w-grid) 
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  hyp                    !< hydrostatic pressure 
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  inflow_damping_factor  !< used for turbulent inflow (non-cyclic boundary conditions) 
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ptdf_x                 !< damping factor for potential temperature in x-direction 
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ptdf_y                 !< damping factor for potential temperature in y-direction 
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_init                !< initial profile of potential temperature 
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  q_init                 !< initial profile of total water mixing ratio 
-                                                                   !< (or total water content with active cloud physics) 
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  rdf                    !< rayleigh damping factor for velocity components 
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  rdf_sc                 !< rayleigh damping factor for scalar quantities 
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ref_state              !< reference state of potential temperature 
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ddzu                   !< 1/dzu
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ddzu_pres              !< modified ddzu for pressure solver
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  dd2zu                  !< 1/(dzu(k)+dzu(k+1))
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  dzu                    !< vertical grid size (u-grid)
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ddzw                   !< 1/dzw
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  dzw                    !< vertical grid size (w-grid)
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  hyp                    !< hydrostatic pressure
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  inflow_damping_factor  !< used for turbulent inflow (non-cyclic boundary conditions)
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ptdf_x                 !< damping factor for potential temperature in x-direction
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ptdf_y                 !< damping factor for potential temperature in y-direction
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_init                !< initial profile of potential temperature
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  q_init                 !< initial profile of total water mixing ratio
+                                                                   !< (or total water content with active cloud physics)
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  rdf                    !< rayleigh damping factor for velocity components
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  rdf_sc                 !< rayleigh damping factor for scalar quantities
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ref_state              !< reference state of potential temperature
                                                                    !< (and density in case of ocean simulation)
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  s_init                 !< initial profile of passive scalar concentration
@@ -215,10 +219,10 @@
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_diss           !< artificial numerical dissipation flux at south face of grid box - TKE dissipation
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_e              !< artificial numerical dissipation flux at south face of grid box - subgrid-scale TKE
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_nc             !< artificial numerical dissipation flux at south face of grid box - clouddrop-number concentration    
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_nr             !< artificial numerical dissipation flux at south face of grid box - raindrop-number concentration    
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_pt             !< artificial numerical dissipation flux at south face of grid box - potential temperature 
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_q              !< artificial numerical dissipation flux at south face of grid box - mixing ratio 
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_qc             !< artificial numerical dissipation flux at south face of grid box - cloudwater 
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_qr             !< artificial numerical dissipation flux at south face of grid box - rainwater 
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_nc             !< artificial numerical dissipation flux at south face of grid box - clouddrop-number concentration
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_nr             !< artificial numerical dissipation flux at south face of grid box - raindrop-number concentration
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_pt             !< artificial numerical dissipation flux at south face of grid box - potential temperature
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_q              !< artificial numerical dissipation flux at south face of grid box - mixing ratio
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_qc             !< artificial numerical dissipation flux at south face of grid box - cloudwater
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_qr             !< artificial numerical dissipation flux at south face of grid box - rainwater
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_s              !< artificial numerical dissipation flux at south face of grid box - passive scalar
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_sa             !< artificial numerical dissipation flux at south face of grid box - salinity
@@ -230,8 +234,8 @@
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_diss           !< 6th-order advective flux at south face of grid box - TKE dissipation
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_e              !< 6th-order advective flux at south face of grid box - subgrid-scale TKE
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_nc             !< 6th-order advective flux at south face of grid box - clouddrop-number concentration 
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_nr             !< 6th-order advective flux at south face of grid box - raindrop-number concentration 
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_pt             !< 6th-order advective flux at south face of grid box - potential temperature 
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_q              !< 6th-order advective flux at south face of grid box - mixing ratio 
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_nc             !< 6th-order advective flux at south face of grid box - clouddrop-number concentration
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_nr             !< 6th-order advective flux at south face of grid box - raindrop-number concentration
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_pt             !< 6th-order advective flux at south face of grid box - potential temperature
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_q              !< 6th-order advective flux at south face of grid box - mixing ratio
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_qc             !< 6th-order advective flux at south face of grid box - cloudwater
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_qr             !< 6th-order advective flux at south face of grid box - rainwater
@@ -246,9 +250,9 @@
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  mean_inflow_profiles  !< used for turbulent inflow (non-cyclic boundary conditions)
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  precipitation_amount  !< precipitation amount due to gravitational settling (bulk microphysics)
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pt_slope_ref          !< potential temperature in rotated coordinate system 
-                                                                    !< (in case of sloped surface) 
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pt_slope_ref          !< potential temperature in rotated coordinate system
+                                                                    !< (in case of sloped surface)
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  total_2d_a            !< horizontal array to store the total domain data, used for atmosphere-ocean coupling (atmosphere data)
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  total_2d_o            !< horizontal array to store the total domain data, used for atmosphere-ocean coupling (ocean data)
-    
+
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  d           !< divergence
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  de_dx       !< gradient of sgs tke in x-direction (lpm)
@@ -281,5 +285,5 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_v    !< 6th-order advective flux at south face of grid box - v-component
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_w    !< 6th-order advective flux at south face of grid box - w-component
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  kh  !< eddy diffusivity for heat 
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  kh  !< eddy diffusivity for heat
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  km  !< eddy diffusivity for momentum
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  prr         !< rain rate
@@ -309,8 +313,8 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  prho_1  !< pointer for swapping of timelevels for respective quantity
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nc_1    !< pointer for swapping of timelevels for respective quantity
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nc_2    !< pointer for swapping of timelevels for respective quantity 
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nc_2    !< pointer for swapping of timelevels for respective quantity
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nc_3    !< pointer for swapping of timelevels for respective quantity
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nr_1    !< pointer for swapping of timelevels for respective quantity
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nr_2    !< pointer for swapping of timelevels for respective quantity 
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nr_2    !< pointer for swapping of timelevels for respective quantity
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nr_3    !< pointer for swapping of timelevels for respective quantity
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  pt_1    !< pointer for swapping of timelevels for respective quantity
@@ -426,5 +430,5 @@
 !------------------------------------------------------------------------------!
  MODULE averaging
- 
+
     USE kinds
 
@@ -458,6 +462,6 @@
                                                                       !< (or total water content with active cloud physics)
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qc_av         !< avg. cloud water content
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ql_av         !< avg. liquid water content 
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ql_c_av       !< avg. change in liquid water content due to 
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ql_av         !< avg. liquid water content
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ql_c_av       !< avg. change in liquid water content due to
                                                                       !< condensation/evaporation during last time step
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ql_v_av       !< avg. volume of liquid water
@@ -468,12 +472,12 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_av          !< avg. passive scalar
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  sa_av         !< avg. salinity
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  u_av          !< avg. horizontal velocity component u 
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  u_av          !< avg. horizontal velocity component u
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  v_av          !< avg. horizontal velocity component v
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vpt_av        !< avg. virtual potential temperature
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  w_av          !< avg. vertical velocity component
-  
+
  END MODULE averaging
 
- 
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -489,5 +493,5 @@
        LOGICAL ::  opened_before  !< file is currently closed, but has been openend before
     END TYPE file_status
-    
+
     INTEGER, PARAMETER      ::  mask_xyz_dimension = 100  !< limit of mask dimensions (100 points in each direction)
     INTEGER, PARAMETER      ::  max_masks = 50            !< maximum number of masks
@@ -532,5 +536,5 @@
     CHARACTER (LEN=20)   ::  bc_uv_b = 'dirichlet'                        !< namelist parameter
     CHARACTER (LEN=20)   ::  bc_uv_t = 'dirichlet'                        !< namelist parameter
-    CHARACTER (LEN=20)   ::  coupling_mode = 'uncoupled'                  !< coupling mode for atmosphere-ocean coupling  
+    CHARACTER (LEN=20)   ::  coupling_mode = 'uncoupled'                  !< coupling mode for atmosphere-ocean coupling
     CHARACTER (LEN=20)   ::  coupling_mode_remote = 'uncoupled'           !< coupling mode of the remote process in case of coupled atmosphere-ocean runs
     CHARACTER (LEN=20)   ::  dissipation_1d = 'detering'                  !< namelist parameter
@@ -539,8 +543,8 @@
     CHARACTER (LEN=20)   ::  random_generator = 'random-parallel'         !< namelist parameter
     CHARACTER (LEN=80)   ::  recycling_method_for_thermodynamic_quantities = 'turbulent_fluctuation'        !< namelist parameter
-    CHARACTER (LEN=20)   ::  reference_state = 'initial_profile'          !< namelist parameter  
-    CHARACTER (LEN=20)   ::  timestep_scheme = 'runge-kutta-3'            !< namelist parameter       
+    CHARACTER (LEN=20)   ::  reference_state = 'initial_profile'          !< namelist parameter
+    CHARACTER (LEN=20)   ::  timestep_scheme = 'runge-kutta-3'            !< namelist parameter
     CHARACTER (LEN=20)   ::  turbulence_closure = 'Moeng_Wyngaard'        !< namelist parameter
-    CHARACTER (LEN=40)   ::  topography = 'flat'                          !< namelist parameter 
+    CHARACTER (LEN=40)   ::  topography = 'flat'                          !< namelist parameter
     CHARACTER (LEN=64)   ::  host = '????'                                !< configuration identifier as given by palmrun option -c, ENVPAR namelist parameter provided by palmrun
     CHARACTER (LEN=80)   ::  log_message                                  !< user-defined message for debugging (sse data_log.f90)
@@ -554,17 +558,17 @@
     CHARACTER (LEN=varnamelength), DIMENSION(500) ::  data_output = ' '       !< namelist parameter
     CHARACTER (LEN=varnamelength), DIMENSION(500) ::  data_output_user = ' '  !< namelist parameter
-    CHARACTER (LEN=varnamelength), DIMENSION(500) ::  doav = ' '              !< label array for multi-dimensional, 
+    CHARACTER (LEN=varnamelength), DIMENSION(500) ::  doav = ' '              !< label array for multi-dimensional,
                                                                               !< averaged output quantities
-                                            
+
     CHARACTER (LEN=varnamelength), DIMENSION(max_masks,100) ::  data_output_masks = ' '       !< namelist parameter
     CHARACTER (LEN=varnamelength), DIMENSION(max_masks,100) ::  data_output_masks_user = ' '  !< namelist parameter
 
     CHARACTER (LEN=varnamelength), DIMENSION(300) ::  data_output_pr = ' '  !< namelist parameter
-    
+
     CHARACTER (LEN=varnamelength), DIMENSION(200) ::  data_output_pr_user = ' '  !< namelist parameter
-    
-    CHARACTER (LEN=varnamelength), DIMENSION(max_masks,0:1,100) ::  domask = ' ' !< label array for multi-dimensional, 
+
+    CHARACTER (LEN=varnamelength), DIMENSION(max_masks,0:1,100) ::  domask = ' ' !< label array for multi-dimensional,
                                                                                  !< masked output quantities
-    
+
     CHARACTER (LEN=varnamelength), DIMENSION(0:1,500) ::  do2d = ' '  !< label array for 2d output quantities
     CHARACTER (LEN=varnamelength), DIMENSION(0:1,500) ::  do3d = ' '  !< label array for 3d output quantities
@@ -572,5 +576,5 @@
     INTEGER(iwp), PARAMETER ::  fl_max = 500     !< maximum number of virtual-flight measurements
     INTEGER(iwp), PARAMETER ::  var_fl_max = 20  !< maximum number of different sampling variables in virtual flight measurements
-    
+
     INTEGER(iwp) ::  abort_mode = 1                    !< abort condition (nested runs)
     INTEGER(iwp) ::  agt_time_count = 0                !< number of output intervals for agent data output
@@ -642,5 +646,5 @@
     INTEGER(iwp) ::  timestep_count = 0                !< number of timesteps carried out since the beginning of the initial run
     INTEGER(iwp) ::  y_shift = 0                       !< namelist parameter
-    
+
     INTEGER(iwp) ::  dist_nxl(0:1)                               !< left boundary of disturbance region
     INTEGER(iwp) ::  dist_nxr(0:1)                               !< right boundary of disturbance region
@@ -662,5 +666,5 @@
     INTEGER(iwp) ::  pt_vertical_gradient_level_ind(10) = -9999  !< grid index values of pt_vertical_gradient_level(s)
     INTEGER(iwp) ::  q_vertical_gradient_level_ind(10) = -9999   !< grid index values of q_vertical_gradient_level(s)
-    INTEGER(iwp) ::  s_vertical_gradient_level_ind(10) = -9999   !< grid index values of s_vertical_gradient_level(s)   
+    INTEGER(iwp) ::  s_vertical_gradient_level_ind(10) = -9999   !< grid index values of s_vertical_gradient_level(s)
     INTEGER(iwp) ::  section(100,3)                              !< collective array for section_xy/xz/yz
     INTEGER(iwp) ::  section_xy(100) = -9999                     !< namelist parameter
@@ -739,4 +743,5 @@
     LOGICAL ::  humidity_remote = .FALSE.                        !< switch for receiving near-surface humidity flux (atmosphere-ocean coupling)
     LOGICAL ::  indoor_model = .FALSE.                           !< switch for indoor-climate and energy-demand model
+    LOGICAL ::  kolmogorov_length_scale = .FALSE.                !< switch to activate calculations in flow_statistics for the kolmogorov length scale
     LOGICAL ::  large_scale_forcing = .FALSE.                    !< namelist parameter
     LOGICAL ::  large_scale_subsidence = .FALSE.                 !< namelist parameter
@@ -750,5 +755,5 @@
     LOGICAL ::  mg_switch_to_pe0 = .FALSE.                       !< internal multigrid switch for steering the ghost point exchange in case that data has been collected on PE0
     LOGICAL ::  monotonic_limiter_z = .FALSE.                    !< use monotonic flux limiter for vertical scalar advection
-    LOGICAL ::  nesting_offline = .FALSE.                        !< flag controlling offline nesting in COSMO model  
+    LOGICAL ::  nesting_offline = .FALSE.                        !< flag controlling offline nesting in COSMO model
     LOGICAL ::  neutral = .FALSE.                                !< namelist parameter
     LOGICAL ::  nudging = .FALSE.                                !< namelist parameter
@@ -804,7 +809,7 @@
     LOGICAL, DIMENSION(max_masks) ::  mask_surface = .FALSE.   !< flag for surface-following masked output
 
-    REAL(wp) ::  advected_distance_x = 0.0_wp                  !< advected distance of model domain along x 
-                                                               !< (galilei transformation) 
-    REAL(wp) ::  advected_distance_y = 0.0_wp                  !< advected distance of model domain along y 
+    REAL(wp) ::  advected_distance_x = 0.0_wp                  !< advected distance of model domain along x
+                                                               !< (galilei transformation)
+    REAL(wp) ::  advected_distance_y = 0.0_wp                  !< advected distance of model domain along y
                                                                !< (galilei transformation)
     REAL(wp) ::  alpha_surface = 0.0_wp                        !< namelist parameter
@@ -830,6 +835,6 @@
     REAL(wp) ::  cos_alpha_surface                             !< cosine of alpha_surface
     REAL(wp) ::  coupling_start_time = 0.0_wp                  !< namelist parameter
-    REAL(wp) ::  days_since_reference_point = 0.0_wp           !< days after atmosphere-ocean coupling has been activated, 
-                                                               !< or after spinup phase of LSM has been finished 
+    REAL(wp) ::  days_since_reference_point = 0.0_wp           !< days after atmosphere-ocean coupling has been activated,
+                                                               !< or after spinup phase of LSM has been finished
     REAL(wp) ::  disturbance_amplitude = 0.25_wp               !< namelist parameter
     REAL(wp) ::  disturbance_energy_limit = 0.01_wp            !< namelist parameter
@@ -886,5 +891,5 @@
     REAL(wp) ::  pt_damping_width = 0.0_wp                     !< namelist parameter
     REAL(wp) ::  pt_reference = 9999999.9_wp                   !< namelist parameter
-    REAL(wp) ::  pt_slope_offset = 0.0_wp                      !< temperature difference between left and right 
+    REAL(wp) ::  pt_slope_offset = 0.0_wp                      !< temperature difference between left and right
                                                                !< boundary of total domain
     REAL(wp) ::  pt_surface = 300.0_wp                         !< namelist parameter
@@ -935,5 +940,5 @@
     REAL(wp) ::  time_do3d = 0.0_wp                            !< time since last 3d output
     REAL(wp) ::  time_do_av = 0.0_wp                           !< time since last averaged-data output
-    REAL(wp) ::  time_do_sla = 0.0_wp                          !< time since last 
+    REAL(wp) ::  time_do_sla = 0.0_wp                          !< time since last
     REAL(wp) ::  time_restart = 9999999.9_wp                   !< time at which run shall be terminated and restarted
     REAL(wp) ::  time_run_control = 0.0_wp                     !< time since last RUN_CONTROL output
@@ -949,6 +954,6 @@
     REAL(wp) ::  tunnel_width_y = 9999999.9_wp                 !< namelist parameter
     REAL(wp) ::  tunnel_wall_depth = 9999999.9_wp              !< namelist parameter
-    REAL(wp) ::  ug_surface = 0.0_wp                           !< namelist parameter 
-    REAL(wp) ::  u_bulk = 0.0_wp                               !< namelist parameter 
+    REAL(wp) ::  ug_surface = 0.0_wp                           !< namelist parameter
+    REAL(wp) ::  u_bulk = 0.0_wp                               !< namelist parameter
     REAL(wp) ::  u_gtrans = 0.0_wp                             !< transformed wind component (galilei transformation)
     REAL(wp) ::  vg_surface = 0.0_wp                           !< namelist parameter
@@ -995,6 +1000,6 @@
     REAL(wp) ::  wall_humidityflux(0:5) = 0.0_wp                   !< namelist parameter
     REAL(wp) ::  wall_salinityflux(0:5) = 0.0_wp                   !< namelist parameter
-    REAL(wp) ::  wall_scalarflux(0:5) = 0.0_wp                     !< namelist parameter 
-    REAL(wp) ::  subs_vertical_gradient(10) = 0.0_wp               !< namelist parameter 
+    REAL(wp) ::  wall_scalarflux(0:5) = 0.0_wp                     !< namelist parameter
+    REAL(wp) ::  subs_vertical_gradient(10) = 0.0_wp               !< namelist parameter
     REAL(wp) ::  subs_vertical_gradient_level(10) = -9999999.9_wp  !< namelist parameter
 
@@ -1004,9 +1009,9 @@
     REAL(wp), DIMENSION(max_masks,mask_xyz_dimension) ::  mask_y = -1.0_wp  !< namelist parameter
     REAL(wp), DIMENSION(max_masks,mask_xyz_dimension) ::  mask_z = -1.0_wp  !< namelist parameter
-    
+
     REAL(wp), DIMENSION(max_masks,3) ::  mask_x_loop = -1.0_wp  !< namelist parameter
     REAL(wp), DIMENSION(max_masks,3) ::  mask_y_loop = -1.0_wp  !< namelist parameter
     REAL(wp), DIMENSION(max_masks,3) ::  mask_z_loop = -1.0_wp  !< namelist parameter
-    
+
 !
 !--    internal mask arrays ("mask,dimension,selection")
@@ -1030,5 +1035,5 @@
     REAL(wp) ::  ddx          !< 1/dx
     REAL(wp) ::  ddx2         !< 1/dx2
-    REAL(wp) ::  dx = 1.0_wp  !< horizontal grid size (along x-direction) 
+    REAL(wp) ::  dx = 1.0_wp  !< horizontal grid size (along x-direction)
     REAL(wp) ::  dx2          !< dx*dx
     REAL(wp) ::  ddy          !< 1/dy
@@ -1040,6 +1045,6 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  ddy2_mg  !< 1/dy_l**2 (dy_l: grid spacing along y on different multigrid level)
 
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zu_s_inner  !< height of topography top on scalar grid 
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zw_w_inner  !< height of topography top on w grid 
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zu_s_inner  !< height of topography top on scalar grid
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zw_w_inner  !< height of topography top on w grid
 
     SAVE
@@ -1090,5 +1095,5 @@
     INTEGER(idp), DIMENSION(:), ALLOCATABLE ::  ngp_3d        !< number of grid points of the total domain
     INTEGER(idp), DIMENSION(:), ALLOCATABLE ::  ngp_3d_inner  !< ! need to have 64 bit for grids > 2E9
-                    
+
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ngp_2dh  !< number of grid points of a horizontal cross section through the total domain
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nxl_mg   !< left-most grid index of subdomain on different multigrid level
@@ -1105,14 +1110,14 @@
     INTEGER(iwp), DIMENSION(:,:,:), POINTER ::  flags  !< pointer to wall_flags_1-10
 
-    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE,  TARGET ::  wall_flags_1   !< topograpyh masking flag on multigrid level 1 
-    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE,  TARGET ::  wall_flags_2   !< topograpyh masking flag on multigrid level 2 
-    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE,  TARGET ::  wall_flags_3   !< topograpyh masking flag on multigrid level 3 
-    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE,  TARGET ::  wall_flags_4   !< topograpyh masking flag on multigrid level 4 
-    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE,  TARGET ::  wall_flags_5   !< topograpyh masking flag on multigrid level 5 
-    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE,  TARGET ::  wall_flags_6   !< topograpyh masking flag on multigrid level 6 
-    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE,  TARGET ::  wall_flags_7   !< topograpyh masking flag on multigrid level 7 
-    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE,  TARGET ::  wall_flags_8   !< topograpyh masking flag on multigrid level 8 
-    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE,  TARGET ::  wall_flags_9   !< topograpyh masking flag on multigrid level 9 
-    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE,  TARGET ::  wall_flags_10  !< topograpyh masking flag on multigrid level 10 
+    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE,  TARGET ::  wall_flags_1   !< topograpyh masking flag on multigrid level 1
+    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE,  TARGET ::  wall_flags_2   !< topograpyh masking flag on multigrid level 2
+    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE,  TARGET ::  wall_flags_3   !< topograpyh masking flag on multigrid level 3
+    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE,  TARGET ::  wall_flags_4   !< topograpyh masking flag on multigrid level 4
+    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE,  TARGET ::  wall_flags_5   !< topograpyh masking flag on multigrid level 5
+    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE,  TARGET ::  wall_flags_6   !< topograpyh masking flag on multigrid level 6
+    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE,  TARGET ::  wall_flags_7   !< topograpyh masking flag on multigrid level 7
+    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE,  TARGET ::  wall_flags_8   !< topograpyh masking flag on multigrid level 8
+    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE,  TARGET ::  wall_flags_9   !< topograpyh masking flag on multigrid level 9
+    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE,  TARGET ::  wall_flags_10  !< topograpyh masking flag on multigrid level 10
 
     INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  advc_flags_m            !< flags used to degrade order of advection scheme for momentum
@@ -1221,5 +1226,5 @@
     CHARACTER(LEN=2) ::  send_receive = 'al'     !<
     CHARACTER(LEN=7) ::  myid_char = ''          !< character string containing processor id number
-    
+
     INTEGER(iwp) ::  comm1dx                     !< communicator for domain decomposition along x
     INTEGER(iwp) ::  comm1dy                     !< communicator for domain decomposition along y
@@ -1256,12 +1261,12 @@
     INTEGER(iwp) ::  tasks_per_node = -9999      !< MPI tasks per compute node
     INTEGER(iwp) ::  threads_per_task = 1        !< number of OPENMP threads per MPI task
-    INTEGER(iwp) ::  type_x                      !< derived MPI datatype for 2-D ghost-point exchange - north / south 
-    INTEGER(iwp) ::  type_xy                     !< derived MPI datatype for 2-D ghost-point exchange - north / south 
-    INTEGER(iwp) ::  type_y                      !< derived MPI datatype for 2-D exchange in atmosphere-ocean coupler 
+    INTEGER(iwp) ::  type_x                      !< derived MPI datatype for 2-D ghost-point exchange - north / south
+    INTEGER(iwp) ::  type_xy                     !< derived MPI datatype for 2-D ghost-point exchange - north / south
+    INTEGER(iwp) ::  type_y                      !< derived MPI datatype for 2-D exchange in atmosphere-ocean coupler
 
     INTEGER(iwp) ::  pdims(2) = 1  !< number of processors along x-y dimension
     INTEGER(iwp) ::  req(100)      !< MPI return variable indicating if send-receive operation is finished
 
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  hor_index_bounds               !< horizontal index bounds 
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  hor_index_bounds               !< horizontal index bounds
     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  hor_index_bounds_previous_run  !< horizontal index bounds of previous run
 
@@ -1283,10 +1288,10 @@
     INTEGER(iwp) ::  pcoord(2)                !< PE coordinates along x and y
     INTEGER(iwp) ::  status(MPI_STATUS_SIZE)  !< MPI status variable used in various MPI calls
-    
+
     INTEGER(iwp), DIMENSION(MPI_STATUS_SIZE,100) ::  wait_stat  !< MPI status variable used in various MPI calls
-    
+
     INTEGER(iwp) ::  type_x_byte !< derived MPI datatype for 2-D 8-bit integer ghost-point exchange - north / south
     INTEGER(iwp) ::  type_y_byte !< derived MPI datatype for 2-D integer ghost-point exchange - left / right
-    
+
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ngp_xz      !< number of ghost points in xz-plane on different multigrid level
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ngp_xz_int  !< number of ghost points in xz-plane on different multigrid level
@@ -1297,5 +1302,5 @@
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  type_xz_int !< derived MPI datatype for 3-D integer ghost-point exchange - north / south
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  type_y_int  !< derived MPI datatype for 2-D integer ghost-point exchange - left / right
-    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  type_yz     !< derived MPI datatype for 3-D integer ghost-point exchange - left / right 
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  type_yz     !< derived MPI datatype for 3-D integer ghost-point exchange - left / right
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  type_yz_int !< derived MPI datatype for 3-D integer ghost-point exchange - left / right
 
@@ -1354,5 +1359,5 @@
     INTEGER(iwp) ::  dopr_index(300) = 0                !< index number of respective profile quantity
     INTEGER(iwp) ::  dopr_initial_index(300) = 0        !< index number of initial profiles to be output
-                
+
     SAVE
 
@@ -1370,5 +1375,5 @@
     CHARACTER (LEN=40) ::  region(0:9) =  &  !< label for statistic region
                            'total domain                            '
- 
+
     INTEGER(iwp) ::  pr_palm = 200          !< maximum number of output profiles
     INTEGER(iwp) ::  statistic_regions = 0  !< identifier for statistic regions
@@ -1377,8 +1382,8 @@
     INTEGER(iwp) ::  v_max_ijk(3) = -1  !< index values (i,j,k) of location where v_max occurs
     INTEGER(iwp) ::  w_max_ijk(3) = -1  !< index values (i,j,k) of location where w_max occurs
-    
+
     LOGICAL ::  flow_statistics_called = .FALSE.  !< flag that tells other routines if flow statistics was executed
                                                   !< (after each timestep)
-    
+
     REAL(wp) ::  u_max = 0.0_wp  !< maximum of absolute u-veloctiy in entire domain
     REAL(wp) ::  v_max = 0.0_wp  !< maximum of absolute v-veloctiy in entire domain
@@ -1386,7 +1391,7 @@
 
     REAL(wp), DIMENSION(2) ::  z_i  !< inversion height
-    
+
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  mean_surface_level_height  !< mean surface level height for the different statistic regions
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  sums_divnew_l              !< subdomain sum (_l) of divergence after pressure 
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  sums_divnew_l              !< subdomain sum (_l) of divergence after pressure
                                                                        !< solver call (new)
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  sums_divold_l              !< subdomain sum (_l) of divergence before pressure
@@ -1394,5 +1399,5 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  weight_substep             !< weighting factor for substeps in timestepping
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  weight_pres                !< substep weighting factor for pressure solver
-    
+
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums             !< global sum array for the various output quantities
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_salsa_ws_l  !< subdomain sum of vertical salsa flux w's' (5th-order advection scheme only)
@@ -1413,10 +1418,10 @@
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsss_ws_l   !< subdomain sum of vertical passive scalar flux w's' (5th-order advection scheme only)
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_ls_l        !< subdomain sum of large scale forcing and nudging tendencies
-                                              
+
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  hom_sum             !< sum array for horizontal mean
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  rmask               !< REAL flag array (0.0 or 1.0) for statistic regions
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  sums_l              !< subdomain sum (_l) gathered for various quantities
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  sums_l_l            !< subdomain sum (_l) of mixing length from diffusivities
-    
+
     REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  hom  !< horizontal mean of various quantities (profiles/timeseries)
 
@@ -1452,5 +1457,5 @@
     INTEGER(iwp) ::  nzt_y   !< internal index bound for transpositions
     INTEGER(iwp) ::  nzt_yd  !< internal index bound for transpositions
-                
+
     SAVE
 
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 4471)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 4472)
@@ -20,9 +20,12 @@
 ! Current revisions:
 ! ------------------
-! 
-! 
+!
+!
 ! Former revisions:
 ! -----------------
 ! $Id$
+! OPENACC COPYIN directive for ddx and ddy added
+!
+! 4466 2020-03-20 16:14:41Z suehring
 ! Add advection fluxes to ACC copyin
 !
@@ -289,4 +292,7 @@
         ONLY:  flight_measurement
 
+    USE grid_variables,                                                                            &
+        ONLY:  ddx, ddy
+
     USE indices,                                                                                   &
         ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, nzb, nzt
@@ -524,4 +530,9 @@
 !$ACC DATA &
 !$ACC COPYIN(tsc(1:5))
+
+!
+!-- Copy data from grid_variables
+!$ACC DATA &
+!$ACC COPYIN(ddx, ddy)
 
 !
@@ -565,4 +576,11 @@
 !$ACC COPY(sums_wsss_ws_l(nzb:nzt+1,0)) &
 !$ACC COPY(sums_salsa_ws_l(nzb:nzt+1,0))
+
+!
+!-- Next statement is to avoid compiler warnings about unused variables. Please
+!-- remove in case that you are using them. ddx and ddy need to be defined in
+!-- time_integration because of ACC COPYIN directive.
+    ddx = ddx
+    ddy = ddy
 
 #if defined( _OPENACC )
@@ -1727,4 +1745,5 @@
 !$ACC END DATA
 !$ACC END DATA
+!$ACC END DATA
 
 #if defined( __parallel )
