Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 3069)
+++ palm/trunk/SOURCE/Makefile	(revision 3083)
@@ -25,4 +25,7 @@
 # -----------------
 # $Id$
+# add turbulence_closure_mod to parin
+# 
+# 2963 2018-04-12 14:47:44Z suehring
 # Introduce index for vegetation/wall, pavement/green-wall and water/window 
 # surfaces, for clearer access of surface fraction, albedo, emissivity, etc. .
@@ -1230,4 +1233,5 @@
 	spectra_mod.o \
 	synthetic_turbulence_generator_mod.o \
+	turbulence_closure_mod.o \
 	uv_exposure_model_mod.o \
 	vertical_nesting_mod.o \
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 3069)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 3083)
@@ -24,5 +24,8 @@
 ! Former revisions:
 ! -----------------
-! $Id$
+! $Id: check_parameters.f90 2520 2017-10-05 13:50:26Z gronemeier &
+! Add inital profile output for e (TG)
+! 
+! 3065 2018-06-12 07:03:02Z Giersch
 ! dz was replaced by dz(1), error message revised
 ! 
@@ -2455,9 +2458,14 @@
              ENDIF
 
-          CASE ( 'e' )
+          CASE ( 'e', '#e' )
              dopr_index(i)  = 8
              dopr_unit(i)   = 'm2/s2'
              hom(:,2,8,:)   = SPREAD( zu, 2, statistic_regions+1 )
              hom(nzb,2,8,:) = 0.0_wp
+             IF ( data_output_pr(i)(1:1) == '#' )  THEN
+                dopr_initial_index(i) = 8
+                hom(:,2,8,:)          = SPREAD( zu, 2, statistic_regions+1 )
+                data_output_pr(i)     = data_output_pr(i)(2:)
+             ENDIF
 
           CASE ( 'km', '#km' )
Index: palm/trunk/SOURCE/data_output_mask.f90
===================================================================
--- palm/trunk/SOURCE/data_output_mask.f90	(revision 3069)
+++ palm/trunk/SOURCE/data_output_mask.f90	(revision 3083)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! 
+! 
+! 3045 2018-05-28 07:55:41Z Giersch
 ! Error messages revised
 ! 
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 3069)
+++ palm/trunk/SOURCE/header.f90	(revision 3083)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Print RANS-mode constants
+! 
+! 3065 2018-06-12 07:03:02Z Giersch
 ! Header output concerning stretching revised
 ! 
@@ -450,4 +453,7 @@
     USE synthetic_turbulence_generator_mod,                                    &
         ONLY:  stg_header
+
+    USE turbulence_closure_mod,                                                &
+        ONLY:  rans_const_c, rans_const_sigma
 
     IMPLICIT NONE
@@ -1946,5 +1952,7 @@
        IF ( wall_adjustment )  WRITE ( io, 453 )  wall_adjustment_factor
     ENDIF
-
+    IF ( rans_mode )  THEN
+       WRITE ( io, 457 )  rans_const_c, rans_const_sigma
+    ENDIF
 !
 !-- Special actions during the run
@@ -2388,4 +2396,11 @@
 456 FORMAT ('    Day of the year at model start :   day_init = ',I3             &
             /'    UTC time at model start        :   time_utc_init = ',F7.1' s')
+457 FORMAT ('    RANS-mode constants: c_0 = ',F9.5/         &
+            '                         c_1 = ',F9.5/         &
+            '                         c_2 = ',F9.5/         &
+            '                         c_3 = ',F9.5/         &
+            '                         c_4 = ',F9.5/         &
+            '                         sigma_e    = ',F9.5/  &
+            '                         sigma_diss = ',F9.5)
 470 FORMAT (//' Actions during the simulation:'/ &
               ' -----------------------------'/)
Index: palm/trunk/SOURCE/model_1d_mod.f90
===================================================================
--- palm/trunk/SOURCE/model_1d_mod.f90	(revision 3069)
+++ palm/trunk/SOURCE/model_1d_mod.f90	(revision 3083)
@@ -25,4 +25,21 @@
 ! -----------------
 ! $Id$
+! Bugfixes:
+!   - preset te_diss and te_e to avoid runtime errors
+!   - implementation of buoyancy term to dissipation
+!     according to Sogachev et al. (2012)
+!   - where diss_p < 0 set diss_p = 0.1 diss
+!   - calculate progn eq(diss) starting from nzb+1
+! Changes:
+!   - add sig_e to TKE equation
+!   - adjust prognostic equation of diss
+!   - set model constants according to Koblitz (2013)
+!   - renamed c_m to c_0
+!   - rename l_black into l1d_init
+!   - calculate l_grid within init_1d_model and save it as l1d_init
+!   - calculate l1d according to DE85 if dissipation is a prognostic value
+!   - made annotations doxygen-readable
+! 
+! 3049 2018-05-29 13:52:36Z Giersch
 ! Error messages revised
 ! 
@@ -164,5 +181,6 @@
     USE kinds
 
-    USE pegrid
+    USE pegrid,                                                                &
+        ONLY:  myid
         
 
@@ -175,10 +193,10 @@
     LOGICAL ::  stop_dt_1d = .FALSE.             !< termination flag, used in case of too small timestep (1d-model)
 
-    REAL(wp) ::  c_1 = 1.44_wp                 !< model constant
-    REAL(wp) ::  c_2 = 1.92_wp                 !< model constant
-    REAL(wp) ::  c_3 = 1.44_wp                 !< model constant
-    REAL(wp) ::  c_h = 0.0015_wp               !< model constant according to Detering and Etling (1985)
-    REAL(wp) ::  c_m = 0.4_wp                  !< model constant, 0.4 according to Detering and Etling (1985)
-    REAL(wp) ::  c_mu = 0.09_wp                !< model constant
+    REAL(wp) ::  alpha_buoyancy                !< model constant according to Koblitz (2013)
+    REAL(wp) ::  c_0 = 0.03_wp**0.25_wp        !< model constant according to Koblitz (2013)
+    REAL(wp) ::  c_1 = 1.52_wp                 !< model constant according to Koblitz (2013)
+    REAL(wp) ::  c_2 = 1.83_wp                 !< model constant according to Koblitz (2013)
+    REAL(wp) ::  c_3                           !< model constant
+    REAL(wp) ::  c_mu                          !< model constant
     REAL(wp) ::  damp_level_1d = -1.0_wp       !< namelist parameter
     REAL(wp) ::  dt_1d = 60.0_wp               !< dynamic timestep (1d-model)
@@ -187,7 +205,9 @@
     REAL(wp) ::  dt_run_control_1d = 60.0_wp   !< namelist parameter
     REAL(wp) ::  end_time_1d = 864000.0_wp     !< namelist parameter
+    REAL(wp) ::  lambda                        !< maximum mixing length
     REAL(wp) ::  qs1d                          !< characteristic humidity scale (1d-model)
     REAL(wp) ::  simulated_time_1d = 0.0_wp    !< updated simulated time (1d-model)
-    REAL(wp) ::  sig_diss = 1.3_wp             !< model constant
+    REAL(wp) ::  sig_diss = 2.95_wp            !< model constant according to Koblitz (2013)
+    REAL(wp) ::  sig_e = 2.95_wp               !< model constant according to Koblitz (2013)
     REAL(wp) ::  time_pr_1d = 0.0_wp           !< updated simulated time for profile output (1d-model)
     REAL(wp) ::  time_run_control_1d = 0.0_wp  !< updated simulated time for run-control output (1d-model)
@@ -198,5 +218,4 @@
     REAL(wp) ::  z01d                          !< roughness length for momentum (1d-model)
     REAL(wp) ::  z0h1d                         !< roughness length for scalars (1d-model)
-
 
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  diss1d   !< tke dissipation rate (1d-model)
@@ -279,6 +298,4 @@
 
     INTEGER(iwp) ::  k  !< loop index
-    
-    REAL(wp) ::  lambda !< maximum mixing length
 
 !
@@ -324,11 +341,11 @@
 !
 !--       Use the same mixing length as in 3D model (LES-mode)
-          !@todo: rename (delete?) this option
-          ! As the mixing length is different between RANS and LES mode, it
-          ! must be distinguished here between these modes. For RANS mode,
-          ! the mixing length is calculated accoding to Blackadar, which is
-          ! the other option at this point.
-          ! Maybe delete this option entirely (not appropriate in LES case)
-          ! 2018-03-20, gronemeier
+          !> @todo rename (delete?) this option
+          !>  As the mixing length is different between RANS and LES mode, it
+          !>  must be distinguished here between these modes. For RANS mode,
+          !>  the mixing length is calculated accoding to Blackadar, which is
+          !>  the other option at this point.
+          !>  Maybe delete this option entirely (not appropriate in LES case)
+          !>  2018-03-20, gronemeier
           DO  k = nzb+1, nzt
              l1d_init(k)  = ( dx * dy * dzw(k) )**0.33333333333333_wp
@@ -359,5 +376,5 @@
        us1d = 0.1_wp   ! without initial friction the flow would not change
     ELSE
-       diss1d(nzb+1) = 1.0_wp
+       diss1d(nzb+1) = 0.001_wp
        e1d(nzb+1)  = 1.0_wp
        km1d(nzb+1) = 1.0_wp
@@ -372,7 +389,8 @@
 
 !
-!-- Tendencies must be preset in order to avoid runtime errors within the
-!-- first Runge-Kutta step
+!-- Tendencies must be preset in order to avoid runtime errors
+    te_diss  = 0.0_wp
     te_dissm = 0.0_wp
+    te_e  = 0.0_wp
     te_em = 0.0_wp
     te_um = 0.0_wp
@@ -381,5 +399,6 @@
 !
 !-- Set model constant
-    IF ( dissipation_1d == 'as_in_3d_model' )  c_m = 0.1_wp
+    IF ( dissipation_1d == 'as_in_3d_model' )  c_0 = 0.1_wp
+    c_mu = c_0**4
 
 !
@@ -423,5 +442,5 @@
 !-- Determine the time step at the start of a 1D-simulation and
 !-- determine and printout quantities used for run control
-    dt_1d = 1.0_wp
+    dt_1d = 0.01_wp
     CALL run_control_1d
 
@@ -483,5 +502,5 @@
 !--             dissipation rate
                 IF ( dissipation_1d == 'detering' )  THEN
-                   diss1d(k) = c_m**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
+                   diss1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
                 ELSEIF ( dissipation_1d == 'as_in_3d_model' )  THEN
                    diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) &
@@ -497,5 +516,5 @@
                                      kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1)  &
                                    - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k)    &
-                                                 ) * ddzw(k)                   &
+                                                 ) * ddzw(k) / sig_e           &
                                    - diss1d(k)
 
@@ -503,17 +522,24 @@
 !
 !--                dissipation rate
-                   te_diss(k) = km1d(k) *                                      &
+                   IF ( rif1d(k) >= 0.0_wp )  THEN
+                      alpha_buoyancy = 1.0_wp - l1d(k) / lambda
+                   ELSE
+                      alpha_buoyancy = 1.0_wp - ( 1.0_wp + ( c_2 - 1.0_wp )    &
+                                                         / ( c_2 - c_1    ) )  &
+                                              * l1d(k) / lambda
+                   ENDIF
+                   c_3 = ( c_1 - c_2 ) * alpha_buoyancy
+                   te_diss(k) = ( km1d(k) *                                    &
                                   ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2  &
                                   + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2  &
-                                  ) * c_1 * c_mu**0.75 / c_h * f / us1d        &
-                                    * SQRT(e1d(k))                             &
+                                  ) * ( c_1 + (c_2 - c_1) * l1d(k) / lambda )  &
                                   - g / pt_0 * kh1d(k) * flux * c_3            &
-                                    * diss1d(k) / ( e1d(k) + 1.0E-20_wp )      &
-                                  + ( kmzp * ( diss1d(k+1) - diss1d(k) )       &
+                                  - c_2 * diss1d(k)                            &
+                                ) * diss1d(k) / ( e1d(k) + 1.0E-20_wp )        &
+                                + (   kmzp * ( diss1d(k+1) - diss1d(k) )       &
                                            * ddzu(k+1)                         &
                                     - kmzm * ( diss1d(k) - diss1d(k-1) )       &
                                            * ddzu(k)                           &
-                                    ) * ddzw(k) / sig_diss                     &
-                                  - c_2 * diss1d(k)**2 / ( e1d(k) + 1.0E-20_wp )
+                                  ) * ddzw(k) / sig_diss
 
                 ENDIF
@@ -546,5 +572,5 @@
 !--          dissipation rate
              IF ( dissipation_1d == 'detering' )  THEN
-                diss1d(k) = c_m**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
+                diss1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
              ELSEIF ( dissipation_1d == 'as_in_3d_model' )  THEN
                 diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) )  &
@@ -565,4 +591,8 @@
 !--          TKE
              IF ( .NOT. dissipation_1d == 'prognostic' )  THEN
+                !> @query why integrate over 2dz
+                !>   Why is it allowed to integrate over two delta-z for e
+                !>   while for u and v it is not?
+                !>   2018-04-23, gronemeier
                 te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2&
                                     + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2&
@@ -572,5 +602,5 @@
                                      kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1)  &
                                    - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k)    &
-                                                 ) * ddzw(k)                   &
+                                                 ) * ddzw(k) / sig_e           &
                                    - diss1d(k)
              ENDIF
@@ -595,10 +625,4 @@
              ENDDO
 
-             IF ( dissipation_1d == 'prognostic' )  THEN
-                DO  k = nzb_diff, nzt
-                   diss1d_p(k) = diss1d(k) + dt_1d * ( tsc(2) * te_diss(k) + &
-                                                 tsc(3) * te_dissm(k) )
-                ENDDO
-             ENDIF
 !
 !--          Eliminate negative TKE values, which can result from the
@@ -606,4 +630,12 @@
 !--          value is reduced to 10 percent of its old value.
              WHERE ( e1d_p < 0.0_wp )  e1d_p = 0.1_wp * e1d
+
+             IF ( dissipation_1d == 'prognostic' )  THEN
+                DO  k = nzb+1, nzt
+                   diss1d_p(k) = diss1d(k) + dt_1d * ( tsc(2) * te_diss(k) + &
+                                                       tsc(3) * te_dissm(k) )
+                ENDDO
+                WHERE ( diss1d_p < 0.0_wp )  diss1d_p = 0.1_wp * diss1d
+             ENDIF
           ENDIF
 
@@ -643,5 +675,6 @@
                    IF ( dissipation_1d == 'prognostic' )  THEN
                       DO k = nzb+1, nzt
-                         te_dissm(k) = -9.5625_wp * te_diss(k) + 5.3125_wp * te_dissm(k)
+                         te_dissm(k) = -9.5625_wp * te_diss(k)  &
+                                     +  5.3125_wp * te_dissm(k)
                       ENDDO
                    ENDIF
@@ -650,5 +683,4 @@
              ENDIF
           ENDIF
-
 
 !
@@ -701,5 +733,8 @@
 
              ENDIF    ! constant_flux_layer
-
+             !> @todo combine if clauses
+             !>   The previous and following if clauses can be combined into a
+             !>   single clause
+             !>   2018-04-23, gronemeier
 !
 !--          Compute the Richardson-flux numbers,
@@ -789,8 +824,8 @@
 !--             compatibility with the 3D model.
                 IF ( ibc_e_b == 2 )  THEN
-                   e1d(nzb+1) = ( us1d / c_m )**2
+                   e1d(nzb+1) = ( us1d / c_0 )**2
                 ENDIF
                 IF ( dissipation_1d == 'prognostic' )  THEN
-                   e1d(nzb+1) = us1d**2 / SQRT( c_mu )
+                   e1d(nzb+1) = ( us1d / c_0 )**2
                    diss1d(nzb+1) = us1d**3 / ( kappa * zu(nzb+1) )
                    diss1d(nzb) = diss1d(nzb+1)
@@ -829,26 +864,33 @@
 !--          in the dissipation of TKE via l1d_diss. Otherwise, km1d would be
 !--          too large.
-             IF ( mixing_length_1d == 'blackadar' )  THEN
+             IF ( dissipation_1d /= 'prognostic' )  THEN
+                IF ( mixing_length_1d == 'blackadar' )  THEN
+                   DO  k = nzb+1, nzt
+                      IF ( rif1d(k) >= 0.0_wp )  THEN
+                         l1d(k) = l1d_init(k) / ( 1.0_wp + 5.0_wp * rif1d(k) )
+                         l1d_diss(k) = l1d(k)
+                      ELSE
+                         l1d(k) = l1d_init(k)
+                         l1d_diss(k) = l1d_init(k) *                           &
+                                       SQRT( 1.0_wp - 16.0_wp * rif1d(k) )
+                      ENDIF
+                   ENDDO
+                ELSEIF ( mixing_length_1d == 'as_in_3d_model' )  THEN
+                   DO  k = nzb+1, nzt
+                      dpt_dz = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k)
+                      IF ( dpt_dz > 0.0_wp )  THEN
+                         l_stable = 0.76_wp * SQRT( e1d(k) )                   &
+                                  / SQRT( g / pt_init(k) * dpt_dz ) + 1E-5_wp
+                      ELSE
+                         l_stable = l1d_init(k)
+                      ENDIF
+                      l1d(k) = MIN( l1d_init(k), l_stable )
+                      l1d_diss(k) = l1d(k)
+                   ENDDO
+                ENDIF
+             ELSE
                 DO  k = nzb+1, nzt
-                   IF ( rif1d(k) >= 0.0_wp )  THEN
-                      l1d(k) = l1d_init(k) / ( 1.0_wp + 5.0_wp * rif1d(k) )
-                      l1d_diss(k) = l1d(k)
-                   ELSE
-                      l1d(k) = l1d_init(k)
-                      l1d_diss(k) = l1d_init(k) *                              &
-                                    SQRT( 1.0_wp - 16.0_wp * rif1d(k) )
-                   ENDIF
-                ENDDO
-             ELSEIF ( mixing_length_1d == 'as_in_3d_model' )  THEN
-                DO  k = nzb+1, nzt
-                   dpt_dz = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k)
-                   IF ( dpt_dz > 0.0_wp )  THEN
-                      l_stable = 0.76_wp * SQRT( e1d(k) ) /                    &
-                                     SQRT( g / pt_init(k) * dpt_dz ) + 1E-5_wp
-                   ELSE
-                      l_stable = l1d_init(k)
-                   ENDIF
-                   l1d(k) = MIN( l1d_init(k), l_stable )
-                   l1d_diss(k) = l1d(k)
+                   l1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) )                   &
+                          / ( diss1d(k) + 1.0E-30_wp )
                 ENDDO
              ENDIF
@@ -867,4 +909,5 @@
                 ENDIF
              ENDIF
+
              IF ( dissipation_1d == 'prognostic' )  THEN
                 DO  k = nzb_diff, nzt
@@ -873,5 +916,5 @@
              ELSE
                 DO  k = nzb_diff, nzt
-                   km1d(k) = c_m * SQRT( e1d(k) ) * l1d(k)
+                   km1d(k) = c_0 * SQRT( e1d(k) ) * l1d(k)
                 ENDDO
              ENDIF
@@ -1018,11 +1061,16 @@
 
     REAL(wp) ::  dt_diff  !< time step accorind to diffusion criterion
+    REAL(wp) ::  dt_old   !< previous time step
     REAL(wp) ::  fac      !< factor of criterion
     REAL(wp) ::  value    !< auxiliary variable
 
 !
+!-- Save previous time step
+    dt_old = dt_1d
+
+!
 !-- Compute the currently feasible time step according to the diffusion
 !-- criterion. At nzb+1 the half grid length is used.
-    fac = 0.125  !0.35_wp                                                       !### changed from 0.35
+    fac = 0.125
     dt_diff = dt_max_1d
     DO  k = nzb+2, nzt
@@ -1034,6 +1082,10 @@
 
 !
+!-- Limit the new time step to a maximum of 10 times the previous time step
+    dt_1d = MIN( dt_old * 10.0_wp, dt_1d )
+
+!
 !-- Set flag when the time step becomes too small
-    IF ( dt_1d < ( 0.00001_wp * dt_max_1d ) )  THEN
+    IF ( dt_1d < ( 1.0E-15_wp * dt_max_1d ) )  THEN
        stop_dt_1d = .TRUE.
 
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 3069)
+++ palm/trunk/SOURCE/modules.f90	(revision 3083)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! set dt_3d = 0.01
+! 
+! 3065 2018-06-12 07:03:02Z Giersch
 ! Variables concerning stretching introduced or revised
 ! 
@@ -1422,5 +1425,5 @@
     REAL(wp) ::  dt_run_control = 60.0_wp                      !< namelist parameter
     REAL(wp) ::  dt_spinup = 60.0_wp                           !< namelist parameter
-    REAL(wp) ::  dt_3d = 1.0_wp                                !< time step
+    REAL(wp) ::  dt_3d = 0.01_wp                               !< time step
     REAL(wp) ::  dz_max = 1000.0_wp                            !< namelist parameter
     REAL(wp) ::  dz_stretch_factor = 1.08_wp                   !< namelist parameter
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 3069)
+++ palm/trunk/SOURCE/parin.f90	(revision 3083)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Added rans_const_c and rans_const_sigma as input parameters (TG)
+! 
+! 3065 2018-06-12 07:03:02Z Giersch
 ! New initialization parameters added
 ! 
@@ -467,4 +470,7 @@
     USE synthetic_turbulence_generator_mod,                                    &
         ONLY:  stg_parin
+
+    USE turbulence_closure_mod,                                                &
+        ONLY:  rans_const_c, rans_const_sigma
 
     USE urban_surface_mod,                                                     &
@@ -531,5 +537,6 @@
              pt_vertical_gradient_level, q_surface, q_surface_initial_change,  &
              q_vertical_gradient, q_vertical_gradient_level,                   &
-             random_generator, random_heatflux, rans_mode,                     &
+             random_generator, random_heatflux, rans_const_c, rans_const_sigma,&
+             rans_mode,                                                        &
              rayleigh_damping_factor, rayleigh_damping_height,                 &
              recycling_width, recycling_yshift,                                &
@@ -601,5 +608,6 @@
              pt_vertical_gradient_level, q_surface, q_surface_initial_change,  &
              q_vertical_gradient, q_vertical_gradient_level,                   &
-             random_generator, random_heatflux, rans_mode,                     &
+             random_generator, random_heatflux, rans_const_c, rans_const_sigma,&
+             rans_mode,                                                        &
              rayleigh_damping_factor, rayleigh_damping_height,                 &
              recycling_width, recycling_yshift,                                &
Index: palm/trunk/SOURCE/timestep.f90
===================================================================
--- palm/trunk/SOURCE/timestep.f90	(revision 3069)
+++ palm/trunk/SOURCE/timestep.f90	(revision 3083)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! limit dt_3d to be at maximum 2*old_dt; define old_dt at beginning of routine
+! Add km/kh_max
+! 
+! 3049 2018-05-29 13:52:36Z Giersch
 ! Error messages revised
 ! 
@@ -161,4 +165,6 @@
     INTEGER(iwp) ::  j !< 
     INTEGER(iwp) ::  k !< 
+    INTEGER(iwp) ::  km_max_ijk(3) = -1  !< index values (i,j,k) of location where km_max occurs
+    INTEGER(iwp) ::  kh_max_ijk(3) = -1  !< index values (i,j,k) of location where kh_max occurs
 
     LOGICAL ::  stop_dt_local !< local switch for controlling the time stepping
@@ -173,4 +179,6 @@
     REAL(wp) ::  dt_w              !< 
     REAL(wp) ::  dt_w_l            !< 
+    REAL(wp) ::  km_max            !< maximum of Km in entire domain
+    REAL(wp) ::  kh_max            !< maximum of Kh in entire domain
     REAL(wp) ::  u_gtrans_l        !< 
     REAL(wp) ::  u_max_l           !< 
@@ -190,6 +198,8 @@
 
 
-
     CALL cpu_log( log_point(12), 'calculate_timestep', 'start' )
+!
+!--    Save former time step as reference
+       old_dt = dt_3d
 
 !
@@ -318,7 +328,8 @@
 !--    The time step is the minimum of the 3-4 components and the diffusion time
 !--    step minus a reduction (cfl_factor) to be on the safe side.
-!--    The time step must not exceed the maximum allowed value.
+!--    The time step must not exceed the maximum allowed value and must not
+!--    increase by more than a factor of 10.
        dt_3d = cfl_factor * MIN( dt_diff, dt_u, dt_v, dt_w, dt_precipitation )
-       dt_3d = MIN( dt_3d, dt_max )
+       dt_3d = MIN( dt_3d, dt_max, 2.0_wp * old_dt )
 
 !
@@ -334,4 +345,12 @@
        IF ( dt_3d < ( 0.00001_wp * dt_max ) )  THEN
           stop_dt = .TRUE.
+
+!
+!--       Determine the maxima of the diffusion coefficients, including their
+!--       grid index positions.
+          CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, km, 'abs',  &
+                               0.0_wp, km_max, km_max_ijk )
+          CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, kh, 'abs',  &
+                               0.0_wp, kh_max, kh_max_ijk )
 
           WRITE( message_string, * ) 'Time step has reached minimum limit.',   &
@@ -342,10 +361,14 @@
                '&dt_w            = ', dt_w, ' s',                              &
                '&dt_diff         = ', dt_diff, ' s',                           &
-               '&u_max           = ', u_max, ' m/s   k=', u_max_ijk(1),        &
+               '&u_max           = ', u_max, ' m/s    k=', u_max_ijk(1),       &
                '  j=', u_max_ijk(2), '  i=', u_max_ijk(3),                     &
-               '&v_max           = ', v_max, ' m/s   k=', v_max_ijk(1),        &
+               '&v_max           = ', v_max, ' m/s    k=', v_max_ijk(1),       &
                '  j=', v_max_ijk(2), '  i=', v_max_ijk(3),                     &
-               '&w_max           = ', w_max, ' m/s   k=', w_max_ijk(1),        &
-               '  j=', w_max_ijk(2), '  i=', w_max_ijk(3)
+               '&w_max           = ', w_max, ' m/s    k=', w_max_ijk(1),       &
+               '  j=', w_max_ijk(2), '  i=', w_max_ijk(3),                     &
+               '&km_max          = ', km_max, ' m2/s2  k=', km_max_ijk(1),     &
+               '  j=', km_max_ijk(2), '  i=', km_max_ijk(3),                   &
+               '&kh_max          = ', kh_max, ' m2/s2  k=', kh_max_ijk(1),     &
+                '  j=', kh_max_ijk(2), '  i=', kh_max_ijk(3)
           CALL message( 'timestep', 'PA0312', 0, 1, 0, 6, 0 )
 !
@@ -388,8 +411,4 @@
        dt_3d = NINT( dt_3d * 100.0_wp / div ) * div / 100.0_wp
 
-!
-!--    Adjust the time step
-       old_dt = dt_3d
-
     ENDIF
 
Index: palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3069)
+++ palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3083)
@@ -25,4 +25,25 @@
 ! -----------------
 ! $Id$
+! - set limits of diss at the end of prognostic equations
+! - call production_e to calculate production term of diss
+! - limit change of diss to -90% to +100%
+! - remove factor 0.5 from diffusion_diss_ij
+! - rename c_m into c_0, and c_h into c_4
+! - add rans_const_c and rans_const_sigma as namelist parameters
+! - add calculation of mixing length for profile output in case of rans_tke_e
+! - changed format of annotations to comply with doxygen standards
+! - calculate and save dissipation rate during rans_tke_l mode
+! - set bc at vertical walls for e, diss, km, kh
+! - bugfix: set l_wall = 0.0 within buildings
+! - set l_wall at bottom and top boundary (rans-mode)
+! - bugfix in production term for dissipation rate
+! - bugfix in diffusion of dissipation rate
+! - disable check for 1D model if rans_tke_e is used
+! - bugfixes for initialization (rans-mode):
+!    - correction of dissipation-rate formula
+!    - calculate km based on l_wall
+!    - initialize diss if 1D model is not used
+! 
+! 3045 2018-05-28 07:55:41Z Giersch
 ! Error message revised
 ! 
@@ -36,5 +57,5 @@
 ! Further todo's
 ! 
-! 2936 2018-03-27 14:49:27Z suehring
+! 2936 2018-03-27 14:49:27Z gronemeier
 ! - defined l_grid only within this module
 ! - Moved l_wall definition from modules.f90
@@ -81,5 +102,4 @@
 !>       add OpenMP directives whereever possible
 !>       remove debug output variables (dummy1, dummy2, dummy3)
-!> @todo Move initialization of wall-mixing length from init_grid
 !> @todo Check for random disturbances
 !> @note <Enter notes on the module>
@@ -146,15 +166,20 @@
 
 
-    REAL(wp) ::  c_1  = 1.44_wp    !< model constant for RANS mode
-    REAL(wp) ::  c_2  = 1.92_wp    !< model constant for RANS mode
-    REAL(wp) ::  c_3  = 1.44_wp    !< model constant for RANS mode
-    REAL(wp) ::  c_h  = 0.0015_wp  !< model constant for RANS mode
-    REAL(wp) ::  c_m               !< constant used for diffusion coefficient and dissipation (dependent on mode RANS/LES)
-    REAL(wp) ::  c_mu = 0.09_wp    !< model constant for RANS mode
-    REAL(wp) ::  l_max             !< maximum length scale for Blackadar mixing length
-    REAL(wp) ::  sig_e = 1.0_wp    !< factor to calculate Ke from Km
-    REAL(wp) ::  sig_diss = 1.3_wp !< factor to calculate K_diss from Km
-    INTEGER(iwp) ::  surf_e        !< end index of surface elements at given i-j position
-    INTEGER(iwp) ::  surf_s        !< start index of surface elements at given i-j position
+    REAL(wp) ::  c_0                !< constant used for diffusion coefficient and dissipation (dependent on mode RANS/LES)
+    REAL(wp) ::  c_1                !< model constant for RANS mode
+    REAL(wp) ::  c_2                !< model constant for RANS mode
+    REAL(wp) ::  c_3                !< model constant for RANS mode
+    REAL(wp) ::  c_4                !< model constant for RANS mode
+    REAL(wp) ::  l_max              !< maximum length scale for Blackadar mixing length
+    REAL(wp) ::  dsig_e = 1.0_wp    !< factor to calculate Ke from Km (1/sigma_e)
+    REAL(wp) ::  dsig_diss = 1.0_wp !< factor to calculate K_diss from Km (1/sigma_diss)
+    INTEGER(iwp) ::  surf_e         !< end index of surface elements at given i-j position
+    INTEGER(iwp) ::  surf_s         !< start index of surface elements at given i-j position
+
+    REAL(wp), DIMENSION(0:4) :: rans_const_c = &       !< model constants for RANS mode (namelist param)
+       (/ 0.55_wp, 1.44_wp, 1.92_wp, 0.0_wp, 0.0_wp /) !> default values fit for standard-tke-e closure
+
+    REAL(wp), DIMENSION(2) :: rans_const_sigma = &     !< model constants for RANS mode, sigma values (sigma_e, sigma_diss) (namelist param)
+       (/ 1.0_wp, 0.77_wp /)
 
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  l_black    !< mixing length according to Blackadar
@@ -163,10 +188,12 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  l_wall !< near-wall mixing length
 
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: dummy1 !< debug output variable
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: dummy2 !< debug output variable
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: dummy3 !< debug output variable
-
-
-    PUBLIC c_m, dummy1, dummy2, dummy3
+    !> @todo remove debug variables
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: diss_prod1, diss_adve1, diss_diff1, &
+                                               diss_prod2, diss_adve2, diss_diff2, &
+                                               diss_prod3, diss_adve3, diss_diff3, &
+                                               dummy1, dummy2, dummy3
+
+
+    PUBLIC c_0, rans_const_c, rans_const_sigma
 
 !
@@ -294,4 +321,9 @@
 ! ------------
 !> Check parameters routine for turbulence closure module.
+!> @todo remove rans_mode from initialization namelist and rework checks
+!>   The way it is implemented at the moment, the user has to set two variables
+!>   so that the RANS mode is working. It would be better if only one parameter
+!>   has to be set.
+!>   2018-06-18, gronemeier
 !------------------------------------------------------------------------------!
  SUBROUTINE tcm_check_parameters
@@ -307,5 +339,14 @@
     IF ( rans_mode )  THEN
 
-       c_m = 0.4_wp  !according to Detering and Etling (1985)
+!
+!--    Assign values to constants for RANS mode
+       dsig_e    = 1.0_wp / rans_const_sigma(1)
+       dsig_diss = 1.0_wp / rans_const_sigma(2)
+
+       c_0 = rans_const_c(0)
+       c_1 = rans_const_c(1)
+       c_2 = rans_const_c(2)
+       !c_3 = rans_const_c(3)   !> @todo clarify how to switch between different models
+       c_4 = rans_const_c(4)
 
        SELECT CASE ( TRIM( turbulence_closure ) )
@@ -316,28 +357,29 @@
           CASE ( 'TKE-e' )
              rans_tke_e = .TRUE.
-
-             IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) == 0  &
-                  .AND.  .NOT.  nest_domain )  THEN
-                message_string = 'Initializing without 1D model while ' //     &
-                                 'using TKE-e closure&' //                     &
-                                 'is not possible at the moment!'
-                CALL message( 'tcm_check_parameters', 'TG0005', 1, 2, 0, 6, 0 )
-             ENDIF
 
           CASE DEFAULT
              message_string = 'Unknown turbulence closure: ' //                &
                               TRIM( turbulence_closure )
-             CALL message( 'tcm_check_parameters', 'TG0001', 1, 2, 0, 6, 0 )
+             CALL message( 'tcm_check_parameters', 'PA0500', 1, 2, 0, 6, 0 )
 
        END SELECT
+
+       IF ( turbulent_inflow .OR. turbulent_outflow )  THEN
+          message_string = 'turbulent inflow/outflow is not yet '//            &
+                           'implemented for RANS mode'
+          CALL message( 'tcm_check_parameters', 'PA0501', 1, 2, 0, 6, 0 )
+       ENDIF
 
        message_string = 'RANS mode is still in development! ' //               &
                         '&Not all features of PALM are yet compatible '//      &
                         'with RANS mode. &Use at own risk!'
-       CALL message( 'tcm_check_parameters', 'TG0003', 0, 1, 0, 6, 0 )
+       CALL message( 'tcm_check_parameters', 'PA0502', 0, 1, 0, 6, 0 )
 
     ELSE
 
-       c_m = 0.1_wp !according to Lilly (1967) and Deardorff (1980)
+       c_0 = 0.1_wp !according to Lilly (1967) and Deardorff (1980)
+
+       dsig_e = 1.0_wp !assure to use K_m to calculate TKE instead 
+                       !of K_e which is used in RANS mode
 
        SELECT CASE ( TRIM( turbulence_closure ) )
@@ -347,19 +389,10 @@
 
           CASE DEFAULT
+             !> @todo rework this part so that only one call of this error exists
              message_string = 'Unknown turbulence closure: ' //                &
                               TRIM( turbulence_closure )
-             CALL message( 'tcm_check_parameters', 'TG0001', 1, 2, 0, 6, 0 )
+             CALL message( 'tcm_check_parameters', 'PA0500', 1, 2, 0, 6, 0 )
 
        END SELECT
-
-    ENDIF
-
-    IF ( rans_tke_e )  THEN
-
-       IF ( turbulent_inflow .OR. turbulent_outflow )  THEN
-          message_string = 'turbulent inflow/outflow is not yet '//            &
-                           'implemented for TKE-e closure'
-          CALL message( 'tcm_check_parameters', 'TG0002', 1, 2, 0, 6, 0 )
-       ENDIF
 
     ENDIF
@@ -379,23 +412,21 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  unit     !< 
-    CHARACTER (LEN=*) ::  var      !<
-
-    INTEGER(iwp) ::  i      !<
-    INTEGER(iwp) ::  ilen   !<
-    INTEGER(iwp) ::  k      !<
+    CHARACTER (LEN=*) ::  unit     !< unit of output variable
+    CHARACTER (LEN=*) ::  var      !< name of output variable
+
+    INTEGER(iwp) ::  i      !< index of var in data_output
+    INTEGER(iwp) ::  ilen   !< length of var string
+    INTEGER(iwp) ::  k      !< flag if var contains one of '_xy', '_xz' or '_yz'
 
     SELECT CASE ( TRIM( var ) )
 
        CASE ( 'diss' )
-          IF ( .NOT.  rans_tke_e )  THEN
-             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                      'res TKE-e closure for RANS mode.'
-             CALL message( 'tcm_check_data_output', 'TG0101', 1, 2, 0, 6, 0 )
-          ENDIF
           unit = 'm2/s3'
 
-       CASE ( 'dummy2', 'dummy3', 'dummy1' )
-          unit = 'mixing length'
+       CASE ( 'diss1', 'diss2',                         &                      !> @todo remove later
+              'diss_prod1', 'diss_adve1', 'diss_diff1', &
+              'diss_prod2', 'diss_adve2', 'diss_diff2', &
+              'diss_prod3', 'diss_adve3', 'diss_diff3', 'dummy3'  )
+          unit = 'debug output'
 
        CASE ( 'kh', 'km' )
@@ -420,11 +451,11 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x   !< 
-    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y   !< 
-    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z   !< 
-    CHARACTER (LEN=*), INTENT(IN)  ::  var      !<
-    
-    LOGICAL, INTENT(OUT) ::  found   !<
-    
+    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x   !< x grid of output variable
+    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y   !< y grid of output variable
+    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z   !< z grid of output variable
+    CHARACTER (LEN=*), INTENT(IN)  ::  var      !< name of output variable
+
+    LOGICAL, INTENT(OUT) ::  found   !< flag if output variable is found
+
     found  = .TRUE.
 
@@ -438,5 +469,8 @@
           grid_z = 'zu'
 
-       CASE ( 'dummy2', 'dummy3', 'dummy1' )                                    !### remove later
+       CASE ( 'diss1', 'diss2',                         &                       !> @todo remove later
+              'diss_prod1', 'diss_adve1', 'diss_diff1', &
+              'diss_prod2', 'diss_adve2', 'diss_diff2', &
+              'diss_prod3', 'diss_adve3', 'diss_diff3', 'dummy3' )
           grid_x = 'x'
           grid_y = 'y'
@@ -480,10 +514,10 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  mode       !< 
-    CHARACTER (LEN=*) ::  variable   !< 
-
-    INTEGER(iwp) ::  i   !< 
-    INTEGER(iwp) ::  j   !< 
-    INTEGER(iwp) ::  k   !< 
+    CHARACTER (LEN=*) ::  mode       !< flag defining mode 'allocate', 'sum' or 'average'
+    CHARACTER (LEN=*) ::  variable   !< name of variable
+
+    INTEGER(iwp) ::  i   !< loop index
+    INTEGER(iwp) ::  j   !< loop index
+    INTEGER(iwp) ::  k   !< loop index
 
     IF ( mode == 'allocate' )  THEN
@@ -616,16 +650,16 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  grid       !< 
-    CHARACTER (LEN=*) ::  mode       !< 
-    CHARACTER (LEN=*) ::  variable   !< 
-
-    INTEGER(iwp) ::  av   !< 
-    INTEGER(iwp) ::  i    !< 
-    INTEGER(iwp) ::  j    !< 
-    INTEGER(iwp) ::  k    !< 
-    INTEGER(iwp) ::  nzb_do   !< 
-    INTEGER(iwp) ::  nzt_do   !< 
-
-    LOGICAL ::  found   !< 
+    CHARACTER (LEN=*) ::  grid       !< name of vertical grid
+    CHARACTER (LEN=*) ::  mode       !< either 'xy', 'xz' or 'yz'
+    CHARACTER (LEN=*) ::  variable   !< name of variable
+
+    INTEGER(iwp) ::  av   !< flag for (non-)average output
+    INTEGER(iwp) ::  i    !< loop index
+    INTEGER(iwp) ::  j    !< loop index
+    INTEGER(iwp) ::  k    !< loop index
+    INTEGER(iwp) ::  nzb_do   !< vertical output index (bottom)
+    INTEGER(iwp) ::  nzt_do   !< vertical output index (top)
+
+    LOGICAL ::  found   !< flag if output variable is found
     LOGICAL ::  two_d   !< flag parameter that indicates 2D variables (horizontal cross sections)
 
@@ -737,14 +771,14 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  variable   !< 
-
-    INTEGER(iwp) ::  av     !< 
-    INTEGER(iwp) ::  i      !< 
-    INTEGER(iwp) ::  j      !< 
-    INTEGER(iwp) ::  k      !< 
+    CHARACTER (LEN=*) ::  variable   !< name of variable
+
+    INTEGER(iwp) ::  av     !< flag for (non-)average output
+    INTEGER(iwp) ::  i      !< loop index
+    INTEGER(iwp) ::  j      !< loop index
+    INTEGER(iwp) ::  k      !< loop index
     INTEGER(iwp) ::  nzb_do !< lower limit of the data output (usually 0)
     INTEGER(iwp) ::  nzt_do !< vertical upper limit of the data output (usually nz_do3d)
 
-    LOGICAL ::  found   !< 
+    LOGICAL ::  found   !< flag if output variable is found
 
     REAL(wp) ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
@@ -829,5 +863,16 @@
           ENDIF
 
-       CASE ( 'dummy1' )                                                        !### remove later
+       CASE ( 'dummy3' )                                                        !> @todo remove later
+          IF ( av == 0 )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_do, nzt_do
+                      local_pf(i,j,k) = dummy3(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'diss1' )                                                         !> @todo remove later
           IF ( av == 0 )  THEN
              DO  i = nxl, nxr
@@ -840,5 +885,5 @@
           ENDIF
 
-       CASE ( 'dummy2' )                                                        !### remove later
+       CASE ( 'diss2' )                                                         !> @todo remove later
           IF ( av == 0 )  THEN
              DO  i = nxl, nxr
@@ -851,10 +896,10 @@
           ENDIF
 
-       CASE ( 'dummy3' )                                                        !### remove later
+       CASE ( 'diss_prod1' )                                                    !> @todo remove later
           IF ( av == 0 )  THEN
              DO  i = nxl, nxr
                 DO  j = nys, nyn
                    DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = dummy3(k,j,i)
+                      local_pf(i,j,k) = diss_prod1(k,j,i)
                    ENDDO
                 ENDDO
@@ -862,4 +907,92 @@
           ENDIF
 
+       CASE ( 'diss_adve1' )                                                    !> @todo remove later
+          IF ( av == 0 )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_do, nzt_do
+                      local_pf(i,j,k) = diss_adve1(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'diss_diff1' )                                                    !> @todo remove later
+          IF ( av == 0 )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_do, nzt_do
+                      local_pf(i,j,k) = diss_diff1(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'diss_prod2' )                                                    !> @todo remove later
+          IF ( av == 0 )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_do, nzt_do
+                      local_pf(i,j,k) = diss_prod2(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'diss_adve2' )                                                    !> @todo remove later
+          IF ( av == 0 )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_do, nzt_do
+                      local_pf(i,j,k) = diss_adve2(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'diss_diff2' )                                                    !> @todo remove later
+          IF ( av == 0 )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_do, nzt_do
+                      local_pf(i,j,k) = diss_diff2(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'diss_prod3' )                                                    !> @todo remove later
+          IF ( av == 0 )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_do, nzt_do
+                      local_pf(i,j,k) = diss_prod3(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'diss_adve3' )                                                    !> @todo remove later
+          IF ( av == 0 )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_do, nzt_do
+                      local_pf(i,j,k) = diss_adve3(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'diss_diff3' )                                                    !> @todo remove later
+          IF ( av == 0 )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_do, nzt_do
+                      local_pf(i,j,k) = diss_diff3(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+          
        CASE DEFAULT
           found = .FALSE.
@@ -891,9 +1024,28 @@
     ALLOCATE( km(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 
-    ALLOCATE( dummy1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )                           !### remove later
+    ALLOCATE( dummy1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )                           !> @todo remove later
     ALLOCATE( dummy2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     ALLOCATE( dummy3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-
-    IF ( rans_mode )  ALLOCATE( l_black(nzb:nzt+1) )
+    ALLOCATE( diss_adve1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+    ALLOCATE( diss_adve2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+    ALLOCATE( diss_adve3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+    ALLOCATE( diss_prod1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+    ALLOCATE( diss_prod2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+    ALLOCATE( diss_prod3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+    ALLOCATE( diss_diff1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+    ALLOCATE( diss_diff2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+    ALLOCATE( diss_diff3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+    dummy1 = 0.0_wp
+    dummy2 = 0.0_wp
+    dummy3 = 0.0_wp
+    diss_adve1 = 0.0_wp
+    diss_adve2 = 0.0_wp
+    diss_adve3 = 0.0_wp
+    diss_prod1 = 0.0_wp
+    diss_prod2 = 0.0_wp
+    diss_prod3 = 0.0_wp
+    diss_diff1 = 0.0_wp
+    diss_diff2 = 0.0_wp
+    diss_diff3 = 0.0_wp
 
 #if defined( __nopointer )
@@ -912,5 +1064,5 @@
 !-- they do not necessarily need to be transferred, which is attributed to 
 !-- the design of the model coupler which allocates memory for each variable. 
-    IF ( rans_tke_e  .OR.  use_sgs_for_particles  .OR.  wang_kernel  .OR.      &
+    IF ( rans_mode  .OR.  use_sgs_for_particles  .OR.  wang_kernel  .OR.       &
          collision_turbulence  .OR.  nested_run )  THEN
 #if defined( __nopointer )
@@ -934,5 +1086,5 @@
     e  => e_1;   e_p  => e_2;   te_m  => e_3
 
-    IF ( rans_tke_e  .OR.  use_sgs_for_particles  .OR.     &
+    IF ( rans_mode  .OR.  use_sgs_for_particles  .OR.     &
          wang_kernel  .OR.  collision_turbulence  .OR.  nested_run )  THEN
        diss => diss_1
@@ -967,10 +1119,11 @@
     INTEGER(iwp) :: j            !< loop index
     INTEGER(iwp) :: k            !< loop index
-    INTEGER(iwp) :: nz_s_shift   !<
-    INTEGER(iwp) :: nz_s_shift_l !< 
+    INTEGER(iwp) :: nz_s_shift   !< lower shift index for scalars
+    INTEGER(iwp) :: nz_s_shift_l !< local lower shift index in case of turbulent inflow
 
 !
 !-- Initialize mixing length
     CALL tcm_init_mixing_length
+    dummy3 = l_wall                 !> @todo remove later
 
 !
@@ -995,7 +1148,7 @@
 
           IF ( rans_tke_e )  THEN
-             IF ( dissipation_1d == 'prognostic' )  THEN    !### Why must this be checked?
-                DO  i = nxlg, nxrg                          !### Should 'diss' not always
-                   DO  j = nysg, nyng                       !### be prognostic in case rans_tke_e?
+             IF ( dissipation_1d == 'prognostic' )  THEN    !> @query Why must this be checked?
+                DO  i = nxlg, nxrg                          !>   Should 'diss' not always
+                   DO  j = nysg, nyng                       !>   be prognostic in case rans_tke_e?
                       diss(:,j,i) = diss1d
                    ENDDO
@@ -1005,5 +1158,5 @@
                    DO  j = nysg, nyng
                       DO  k = nzb+1, nzt
-                         diss(k,j,i) = e(k,j,i) * SQRT( e(k,j,i) ) / l1d(k)
+                         diss(k,j,i) = c_0**4 * e(k,j,i)**2 / km1d(k)
                       ENDDO
                    ENDDO
@@ -1016,15 +1169,19 @@
 
           IF ( constant_diffusion )  THEN
-             km   = km_constant
-             kh   = km / prandtl_number
-             e    = 0.0_wp
+             km = km_constant
+             kh = km / prandtl_number
+             e  = 0.0_wp
           ELSEIF ( e_init > 0.0_wp )  THEN
-             DO  k = nzb+1, nzt
-                km(k,:,:) = c_m * l_grid(k) * SQRT( e_init )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb+1, nzt
+                      km(k,j,i) = c_0 * l_wall(k,j,i) * SQRT( e_init )
+                   ENDDO
+                ENDDO
              ENDDO
              km(nzb,:,:)   = km(nzb+1,:,:)
              km(nzt+1,:,:) = km(nzt,:,:)
-             kh   = km / prandtl_number
-             e    = e_init
+             kh = km / prandtl_number
+             e  = e_init
           ELSE
              IF ( .NOT. ocean )  THEN
@@ -1040,4 +1197,16 @@
           ENDIF
 
+          IF ( rans_tke_e )  THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb+1, nzt
+                      diss(k,j,i) = c_0**4 * e(k,j,i)**2 / km(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+             diss(nzb,:,:) = diss(nzb+1,:,:)
+             diss(nzt+1,:,:) = diss(nzt,:,:)
+          ENDIF
+
        ENDIF
 !
@@ -1073,4 +1242,13 @@
              ENDDO
           ENDDO
+          IF ( rans_tke_e )  THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   nz_s_shift = get_topography_top_index_ji( j, i, 's' )
+
+                   diss(nz_s_shift:nzt+1,j,i) = diss(0:nzt+1-nz_s_shift,j,i)
+                ENDDO
+             ENDDO
+          ENDIF
        ENDIF
 
@@ -1084,5 +1262,6 @@
 !--       boundary and adjust mean inflow profiles
           IF ( complex_terrain )  THEN
-             IF ( nxlg <= 0 .AND. nxrg >= 0 .AND. nysg <= 0 .AND. nyng >= 0 )  THEN
+             IF ( nxlg <= 0 .AND. nxrg >= 0 .AND.  &
+                  nysg <= 0 .AND. nyng >= 0        )  THEN
                 nz_s_shift_l = get_topography_top_index_ji( 0, 0, 's' )
              ELSE
@@ -1095,5 +1274,6 @@
              nz_s_shift = nz_s_shift_l
 #endif
-             mean_inflow_profiles(nz_s_shift:nzt+1,5) = hom_sum(0:nzt+1-nz_s_shift,8,0)  ! e
+             mean_inflow_profiles(nz_s_shift:nzt+1,5) =  &
+                hom_sum(0:nzt+1-nz_s_shift,8,0)  ! e
           ENDIF
 !
@@ -1114,5 +1294,5 @@
 !
 !--       Inside buildings set TKE back to zero.
-!--       Other scalars (km, kh, diss, ...) are ignored at present,
+!--       Other scalars (km, kh,...) are ignored at present,
 !--       maybe revise later.
           DO  i = nxlg, nxrg
@@ -1121,10 +1301,18 @@
                    e(k,j,i)     = MERGE( e(k,j,i), 0.0_wp,                     &
                                          BTEST( wall_flags_0(k,j,i), 0 ) )
-                   te_m(k,j,i)  = MERGE( te_m(k,j,i), 0.0_wp,                  &
-                                         BTEST( wall_flags_0(k,j,i), 0 ) )
                 ENDDO
              ENDDO
           ENDDO
 
+          IF ( rans_tke_e )  THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt
+                      diss(k,j,i)    = MERGE( diss(k,j,i), 0.0_wp,             &
+                                              BTEST( wall_flags_0(k,j,i), 0 ) )
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
        ENDIF
 !
@@ -1134,7 +1322,12 @@
 !
 !--    Allthough tendency arrays are set in prognostic_equations, they have
-!--    have to be predefined here because they are used (but multiplied with 0)
-!--    there before they are set. 
+!--    to be predefined here because there they are used (but multiplied with 0)
+!--    before they are set.
        te_m = 0.0_wp
+
+       IF ( rans_tke_e )  THEN
+          diss_p = diss
+          tdiss_m = 0.0_wp
+       ENDIF
 
     ENDIF
@@ -1215,11 +1408,11 @@
 
        DO  k = 1, nzt
-          IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR.           &
+          IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR.            &
                l_grid(k) > 1.5_wp * dy * wall_adjustment_factor )  THEN
-             WRITE( message_string, * ) 'grid anisotropy exceeds ',            &
-                                        'threshold given by only local',       &
-                                        ' horizontal reduction of near_wall ', &
-                                        'mixing length l_wall',                &
-                                        ' starting from height level k = ', k, &
+             WRITE( message_string, * ) 'grid anisotropy exceeds ',             &
+                                        'threshold given by only local',        &
+                                        ' &horizontal reduction of near_wall ', &
+                                        'mixing length l_wall',                 &
+                                        ' &starting from height level k = ', k, &
                                         '.'
              CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 )
@@ -1320,10 +1513,11 @@
 !
 !-- Initialize the mixing length in case of a RANS simulation
+       ALLOCATE( l_black(nzb:nzt+1) )
 
 !
 !--    Calculate mixing length according to Blackadar (1962)
        IF ( f /= 0.0_wp )  THEN
-          l_max = 2.7E-4 * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) /               &
-                  ABS( f ) + 1E-10_wp
+          l_max = 2.7E-4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) /            &
+                  ABS( f ) + 1.0E-10_wp
        ELSE
           l_max = 30.0_wp
@@ -1338,11 +1532,11 @@
 !
 !--    Gather topography information of whole domain
-       !## TODO: reduce amount of data sent by MPI call
-       !##  By now, a whole global 3D-array is sent and received with
-       !##  MPI_ALLREDUCE although most of the array is 0. This can be
-       !##  drastically reduced if only the local subarray is sent and stored
-       !##  in a global array. For that, an MPI data type or subarray must be
-       !##  defined.
-       !##  2018-03-19, gronemeier
+       !> @todo reduce amount of data sent by MPI call
+       !>   By now, a whole global 3D-array is sent and received with
+       !>   MPI_ALLREDUCE although most of the array is 0. This can be
+       !>   drastically reduced if only the local subarray is sent and stored
+       !>   in a global array. For that, an MPI data type or subarray must be
+       !>   defined.
+       !>   2018-03-19, gronemeier
        ALLOCATE( wall_flags_0_global(nzb:nzt+1,0:ny,0:nx) )
 
@@ -1373,4 +1567,7 @@
           ENDDO
        ENDDO
+
+       l_wall(nzb,:,:) = l_black(nzb)
+       l_wall(nzt+1,:,:) = l_black(nzt+1)
 !
 !--    Limit mixing length to either nearest wall or Blackadar mixing length.
@@ -1415,10 +1612,10 @@
           IF ( rad_k_b /= 0 .OR. rad_k_t /= 0 )  THEN
 
-             !## NOTE: shape of vicinity is larger in z direction
-             !##  Shape of vicinity is two grid points larger than actual search
-             !##  radius in vertical direction. The first and last grid point is
-             !##   always set to 1 to asure correct detection of topography. See
-             !##  function "shortest_distance" for details.
-             !##  2018-03-16, gronemeier
+             !> @note shape of vicinity is larger in z direction
+             !>   Shape of vicinity is two grid points larger than actual search
+             !>   radius in vertical direction. The first and last grid point is
+             !>   always set to 1 to asure correct detection of topography. See
+             !>   function "shortest_distance" for details.
+             !>   2018-03-16, gronemeier
              ALLOCATE( vicinity(-rad_k-1:rad_k+1,-rad_j:rad_j,-rad_i:rad_i) )
              ALLOCATE( vic_yz(0:rad_k+1,0:rad_j) )
@@ -1601,5 +1798,5 @@
                    ELSE  !Check if (i,j,k) belongs to atmosphere
 
-                      l_wall(k,j,i) = -999.0
+                      l_wall(k,j,i) = l_black(k)
 
                    ENDIF
@@ -1630,4 +1827,6 @@
 !> (pos_i/jj/kk), where (jj/kk) is the position of the maximum of 'array'
 !> closest to the origin (0/0) of 'array'.
+!> @todo this part of PALM does not reproduce the same results for optimized
+!>   and debug options for the compiler. This should be fixed
 !------------------------------------------------------------------------------!
     REAL FUNCTION shortest_distance( array, orientation, pos_i )
@@ -1659,28 +1858,28 @@
        IF ( orientation ) THEN  !if array is oriented upwards
           DO  jj = 0, rad_j
-             shortest_distance = MIN( shortest_distance,                        &
-                                      SQRT( MAX(pos_i*dx-0.5*dx,0.0)**2         &
-                                          + MAX(jj*dy-0.5*dy,0.0)**2            &
-                                          + MAX(zw(loc_k(jj)+k-1)-zu(k),0.0)**2 &
-                                          )                                     &
-                                    )
+             shortest_distance =                                               &
+                MIN( shortest_distance,                                        &
+                     SQRT( MAX(REAL(pos_i, KIND=wp)*dx-0.5_wp*dx, 0.0_wp)**2   &
+                         + MAX(REAL(jj, KIND=wp)*dy-0.5_wp*dy, 0.0_wp)**2      &
+                         + MAX(zw(loc_k(jj)+k-1)-zu(k), 0.0_wp)**2             &
+                         )                                                     &
+                   )
           ENDDO
        ELSE  !if array is oriented downwards
-          !## NOTE: MAX within zw required to circumvent error at domain border
-          !##  At the domain border, if non-cyclic boundary is present, the
-          !##  index for zw could be -1, which will be errorneous (zw(-1) does
-          !##  not exist). The MAX function limits the index to be at least 0.
+          !> @note MAX within zw required to circumvent error at domain border
+          !>   At the domain border, if non-cyclic boundary is present, the
+          !>   index for zw could be -1, which will be errorneous (zw(-1) does
+          !>   not exist). The MAX function limits the index to be at least 0.
           DO  jj = 0, rad_j
-             shortest_distance = MIN( shortest_distance,                       &
-                                      SQRT( MAX(pos_i*dx-0.5*dx,0.0)**2        &
-                                          + MAX(jj*dy-0.5*dy,0.0)**2           &
-                                          + MAX(zu(k)-zw(MAX(k-loc_k(jj),      &
-                                                             0_iwp)),          &
-                                                0.0)**2                        &
-                                          )                                    &
-                                    )
+             shortest_distance =                                               &
+                MIN( shortest_distance,                                        &
+                     SQRT( MAX(REAL(pos_i, KIND=wp)*dx-0.5_wp*dx, 0.0_wp)**2   &
+                         + MAX(REAL(jj, KIND=wp)*dy-0.5_wp*dy, 0.0_wp)**2      &
+                         + MAX(zu(k)-zw(MAX(k-loc_k(jj),0_iwp)), 0.0_wp)**2    &
+                         )                                                     &
+                   )
           ENDDO
        ENDIF
-
+       
     END FUNCTION
 
@@ -2004,4 +2203,6 @@
 !
 !--    Use special boundary condition in case of TKE-e closure
+       !> @todo do the same for usm and lsm surfaces
+       !>   2018-06-05, gronemeier
        IF ( rans_tke_e )  THEN
           DO  i = nxl, nxr
@@ -2011,5 +2212,5 @@
                 DO  m = surf_s, surf_e
                    k = surf_def_h(0)%k(m)
-                   e_p(k,j,i) = surf_def_h(0)%us(m)**2 / c_m**2
+                   e_p(k,j,i) = surf_def_h(0)%us(m)**2 / c_0**2
                 ENDDO
              ENDDO
@@ -2092,5 +2293,5 @@
              DO  k = nzb+1, nzt
 !                tend(k,j,i) = tend(k,j,i) + c_1 * diss(k,j,i) / ( e(k,j,i) + 1.0E-20_wp ) * produc(k)
-                tend(k,j,i) = tend(k,j,i) + c_1 * c_mu * f / c_h               &  !### needs revision
+                tend(k,j,i) = tend(k,j,i) + c_1 * c_0**4 * f / c_4               &  !> @todo needs revision
                       / surf_def_h(0)%us(surf_def_h(0)%start_index(j,i))       &
                       * SQRT(e(k,j,i)) * produc(k,j,i)
@@ -2103,7 +2304,7 @@
 !
 !--    Additional sink term for flows through plant canopies
-!        IF ( plant_canopy )  CALL pcm_tendency( ? )                            !### what to do with this?
-
-!        CALL user_actions( 'diss-tendency' )                                   !### not yet implemented
+!        IF ( plant_canopy )  CALL pcm_tendency( ? )                            !> @query what to do with this?
+
+!        CALL user_actions( 'diss-tendency' )                                   !> @todo not yet implemented
 
 !
@@ -2181,8 +2382,12 @@
     USE arrays_3d,                                                             &
         ONLY:  ddzu, diss_l_diss, diss_l_e, diss_s_diss, diss_s_e,             &
-               flux_l_diss, flux_l_e, flux_s_diss, flux_s_e
+               flux_l_diss, flux_l_e, flux_s_diss, flux_s_e,&
+               u_p,v_p,w_p
 
     USE control_parameters,                                                    &
         ONLY:  f, tsc
+
+    USE grid_variables,                                                        &
+        ONLY:  dx, dy
 
     USE surface_mod,                                                           &
@@ -2190,17 +2395,30 @@
                 surf_usm_v 
 
+    use indices, only: nx, ny
+
     IMPLICIT NONE
 
     INTEGER(iwp) ::  i       !< loop index x direction
-    INTEGER(iwp) ::  i_omp   !< 
+    INTEGER(iwp) ::  i_omp   !< first loop index of i-loop in prognostic_equations
     INTEGER(iwp) ::  j       !< loop index y direction
     INTEGER(iwp) ::  k       !< loop index z direction
+    INTEGER(iwp) ::  l       !< loop index
     INTEGER(iwp) ::  m       !< loop index
     INTEGER(iwp) ::  surf_e  !< end index of surface elements at given i-j position
     INTEGER(iwp) ::  surf_s  !< start index of surface elements at given i-j position
-    INTEGER(iwp) ::  tn      !<
-
-    REAL(wp), DIMENSION(nzb:nzt+1) :: advec  !< advection term of TKE tendency
-    REAL(wp), DIMENSION(nzb:nzt+1) :: produc !< production term of TKE tendency
+    INTEGER(iwp) ::  tn      !< task number of openmp task
+
+    INTEGER(iwp) :: pis = 32 !< debug variable, print from i=pis                !> @todo remove later
+    INTEGER(iwp) :: pie = 32 !< debug variable, print until i=pie               !> @todo remove later
+    INTEGER(iwp) :: pjs = 26 !< debug variable, print from j=pjs                !> @todo remove later
+    INTEGER(iwp) :: pje = 26 !< debug variable, print until j=pje               !> @todo remove later
+    INTEGER(iwp) :: pkb = 1  !< debug variable, print from k=pkb                !> @todo remove later
+    INTEGER(iwp) :: pkt = 7  !< debug variable, print until k=pkt               !> @todo remove later
+
+    REAL(wp), DIMENSION(nzb:nzt+1) :: dum_adv   !< debug variable               !> @todo remove later
+    REAL(wp), DIMENSION(nzb:nzt+1) :: dum_pro   !< debug variable               !> @todo remove later
+    REAL(wp), DIMENSION(nzb:nzt+1) :: dum_dif   !< debug variable               !> @todo remove later
+
+5555 FORMAT(A,7(1X,E12.5))   !> @todo remove later
 
 !
@@ -2224,9 +2442,9 @@
        ENDIF
 
-       advec(:) = tend(:,j,i)
-
-       CALL production_e( i, j )
-
-       produc(:) = tend(:,j,i) - advec(:)
+       dum_adv = tend(:,j,i)                                                    !> @todo remove later
+
+       CALL production_e( i, j, .FALSE. )
+
+       dum_pro = tend(:,j,i) - dum_adv                                          !> @todo remove later
 
        IF ( .NOT. humidity )  THEN
@@ -2239,4 +2457,6 @@
           CALL diffusion_e( i, j, vpt, pt_reference )
        ENDIF
+
+       dum_dif = tend(:,j,i) - dum_adv - dum_pro                                !> @todo remove later
 
 !
@@ -2268,5 +2488,14 @@
           DO  m = surf_s, surf_e
              k = surf_def_h(0)%k(m)
-             e_p(k,j,i) = surf_def_h(0)%us(m)**2 / c_m**2
+             e_p(k,j,i) = ( surf_def_h(0)%us(m) / c_0 )**2
+          ENDDO
+
+          DO  l = 0, 3
+             surf_s = surf_def_v(l)%start_index(j,i)
+             surf_e = surf_def_v(l)%end_index(j,i)
+             DO  m = surf_s, surf_e
+                k = surf_def_v(l)%k(m)
+                e_p(k,j,i) = ( surf_def_v(l)%us(m) / c_0 )**2
+             ENDDO
           ENDDO
        ENDIF
@@ -2287,4 +2516,21 @@
           ENDIF
        ENDIF
+
+!        if ( i >= pis .and. i <= pie .and. j >= pjs .and. j <= pje ) then        !> @todo remove later
+!           WRITE(9, *) '------'
+!           WRITE(9, '(A,F8.3,1X,F8.3,1X,I2)') 't, dt, int_ts:', simulated_time, dt_3d, intermediate_timestep_count
+!           WRITE(9, *) 'i:', i
+!           WRITE(9, *) 'j:', j
+!           WRITE(9, *) 'k:', pkb, ' - ', pkt
+!           WRITE(9, *) '---'
+!           WRITE(9, *) 'e:'
+!           WRITE(9, 5555) 'adv :', dum_adv(pkb:pkt)
+!           WRITE(9, 5555) 'pro :', dum_pro(pkb:pkt)
+!           WRITE(9, 5555) 'dif :', dum_dif(pkb:pkt)
+!           WRITE(9, 5555) 'tend:', tend(pkb:pkt,j,i)
+!           WRITE(9, 5555) 'e_p :', e_p(pkb:pkt,j,i)
+!           WRITE(9, 5555) 'e   :', e(pkb:pkt,j,i)
+!           FLUSH(9)
+!        endif
 
     ENDIF   ! TKE equation
@@ -2309,20 +2555,34 @@
        ENDIF
 
+       IF ( intermediate_timestep_count == 1 )  diss_adve1(:,j,i) = tend(:,j,i) !> @todo remove later
+       IF ( intermediate_timestep_count == 2 )  diss_adve2(:,j,i) = tend(:,j,i)
+       IF ( intermediate_timestep_count == 3 )  diss_adve3(:,j,i) = tend(:,j,i)
+
 !
 !--    Production of TKE dissipation rate
-       DO  k = nzb+1, nzt
-!              tend(k,j,i) = tend(k,j,i) + c_1 * diss(k,j,i) / ( e(k,j,i) + 1.0E-20_wp ) * produc(k)
-          tend(k,j,i) = tend(k,j,i) + c_1 * c_mu * f / c_h                     &  !### needs revision
-                      / surf_def_h(0)%us(surf_def_h(0)%start_index(j,i))       &
-                      * SQRT(e(k,j,i)) * produc(k)
-       ENDDO
-
+       CALL production_e( i, j, .TRUE. )
+
+       IF ( intermediate_timestep_count == 1 )  diss_prod1(:,j,i) = tend(:,j,i) - diss_adve1(:,j,i) !> @todo remove later
+       IF ( intermediate_timestep_count == 2 )  diss_prod2(:,j,i) = tend(:,j,i) - diss_adve2(:,j,i)
+       IF ( intermediate_timestep_count == 3 )  diss_prod3(:,j,i) = tend(:,j,i) - diss_adve3(:,j,i)
+
+       dum_pro = tend(:,j,i) - dum_adv                                          !> @todo remove later
+
+!
+!--    Diffusion term of TKE dissipation rate
        CALL diffusion_diss( i, j )
 
+       IF ( intermediate_timestep_count == 1 )  diss_diff1(:,j,i) = tend(:,j,i) - diss_adve1(:,j,i) - diss_prod1(:,j,i) !> @todo remove later
+       IF ( intermediate_timestep_count == 2 )  diss_diff2(:,j,i) = tend(:,j,i) - diss_adve2(:,j,i) - diss_prod2(:,j,i)
+       IF ( intermediate_timestep_count == 3 )  diss_diff3(:,j,i) = tend(:,j,i) - diss_adve3(:,j,i) - diss_prod3(:,j,i)
+       IF ( intermediate_timestep_count == 3 )  dummy3(:,j,i) = km(:,j,i)
+
+       dum_dif = tend(:,j,i) - dum_adv - dum_pro                                !> @todo remove later
+
 !
 !--    Additional sink term for flows through plant canopies
-!        IF ( plant_canopy )  CALL pcm_tendency( i, j, ? )                      !### not yet implemented
-
-!        CALL user_actions( i, j, 'diss-tendency' )                             !### not yet implemented
+!        IF ( plant_canopy )  CALL pcm_tendency( i, j, ? )                      !> @todo not yet implemented
+
+!        CALL user_actions( i, j, 'diss-tendency' )                             !> @todo not yet implemented
 
 !
@@ -2338,5 +2598,4 @@
                                                 BTEST( wall_flags_0(k,j,i), 0 )&
                                                )
-          IF ( diss_p(k,j,i) <= 0.0_wp )  diss_p(k,j,i) = 0.1_wp * diss(k,j,i)
        ENDDO
 
@@ -2350,4 +2609,21 @@
        ENDDO
 
+       DO  l = 0, 1
+          surf_s = surf_def_v(l)%start_index(j,i)
+          surf_e = surf_def_v(l)%end_index(j,i)
+          DO  m = surf_s, surf_e
+             k = surf_def_v(l)%k(m)
+             diss_p(k,j,i) = surf_def_v(l)%us(m)**3 / ( kappa * 0.5_wp * dy )
+          ENDDO
+       ENDDO
+
+       DO  l = 2, 3
+          surf_s = surf_def_v(l)%start_index(j,i)
+          surf_e = surf_def_v(l)%end_index(j,i)
+          DO  m = surf_s, surf_e
+             k = surf_def_v(l)%k(m)
+             diss_p(k,j,i) = surf_def_v(l)%us(m)**3 / ( kappa * 0.5_wp * dx )
+          ENDDO
+       ENDDO
 !
 !--    Calculate tendencies for the next Runge-Kutta step
@@ -2365,8 +2641,30 @@
           ENDIF
        ENDIF
-
-!        IF ( intermediate_timestep_count == 1 )  dummy1(:,j,i) = e_p(:,j,i)
-!        IF ( intermediate_timestep_count == 2 )  dummy2(:,j,i) = e_p(:,j,i)
-!        IF ( intermediate_timestep_count == 3 )  dummy3(:,j,i) = e_p(:,j,i)
+!
+!--    Limit change of diss to be between -90% and +100%. Also, set an absolute
+!--    minimum value
+       DO  k = nzb, nzt+1
+          diss_p(k,j,i) = MIN( MAX( diss_p(k,j,i),         &
+                                    0.1_wp * diss(k,j,i),  &
+                                    0.0001_wp ),           &
+                               2.0_wp * diss(k,j,i) )
+       ENDDO
+
+       IF ( intermediate_timestep_count == 1 )  dummy1(:,j,i) = diss_p(:,j,i)   !> @todo remove later
+       IF ( intermediate_timestep_count == 2 )  dummy2(:,j,i) = diss_p(:,j,i)
+
+!        if ( i >= pis .and. i <= pie .and. j >= pjs .and. j <= pje ) then        !> @todo remove later
+!           WRITE(9, *) '---'
+!           WRITE(9, *) 'diss:'
+!           WRITE(9, 5555) 'adv   :', dum_adv(pkb:pkt)
+!           WRITE(9, 5555) 'pro   :', dum_pro(pkb:pkt)
+!           WRITE(9, 5555) 'dif   :', dum_dif(pkb:pkt)
+!           WRITE(9, 5555) 'tend  :', tend(pkb:pkt,j,i)
+!           WRITE(9, 5555) 'diss_p:', diss_p(pkb:pkt,j,i)
+!           WRITE(9, 5555) 'diss  :', diss(pkb:pkt,j,i)
+!           WRITE(9, *) '---'
+!           WRITE(9, 5555) 'km    :', km(pkb:pkt,j,i)
+!           flush(9)
+!        endif
 
     ENDIF   ! dissipation equation
@@ -2382,4 +2680,5 @@
 !> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is
 !>          not considered well!
+!> @todo Adjust production term in case of rans_tke_e simulation
 !------------------------------------------------------------------------------!
  SUBROUTINE production_e
@@ -3099,6 +3398,7 @@
 !> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is
 !>          not considered well!
+!> @todo non-neutral case is not yet considered for RANS mode
 !------------------------------------------------------------------------------!
- SUBROUTINE production_e_ij( i, j )
+ SUBROUTINE production_e_ij( i, j, diss_production )
 
     USE arrays_3d,                                                             &
@@ -3121,4 +3421,6 @@
 
     IMPLICIT NONE
+
+    LOGICAL :: diss_production
 
     INTEGER(iwp) ::  i       !< running index x-direction
@@ -3153,5 +3455,5 @@
     REAL(wp), DIMENSION(nzb+1:nzt)  ::  dwdy        !< Gradient of w-component in y-direction
     REAL(wp), DIMENSION(nzb+1:nzt)  ::  dwdz        !< Gradient of w-component in z-direction
-
+    REAL(wp), DIMENSION(nzb+1:nzt)  ::  tend_temp   !< temporal tendency
 
     IF ( constant_flux_layer )  THEN
@@ -3348,23 +3650,14 @@
        ENDDO
 
+!        IF ( .NOT. rans_tke_e )  THEN
+
        DO  k = nzb+1, nzt
 
-          def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) +            &
-                           dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2   +            &
-                           dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2   +            &
-                2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + dwdy(k)*dvdz(k) )
-
-!
-!--       Production term according to Kato and Launder (1993)
-!           def = SQRT( ( dudx(k) + dudy(k) + dudz(k) +                       &
-!                         dvdx(k) + dvdy(k) + dvdz(k) +                       &
-!                         dwdx(k) + dwdy(k) + dwdz(k)                         &
-!                       )**4 -                                                &
-!                       ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 +              &
-!                         2.0_wp * ( dudy(k) * dvdx(k) +                      &
-!                                    dudz(k) * dwdx(k) +                      &
-!                                    dvdz(k) * dwdy(k) )                      &
-!                       )**2                                                  &
-!                     )
+          def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) +         &
+                           dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2   +         &
+                           dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2   +         &
+                2.0_wp * ( dvdx(k)*dudy(k) +                                &
+                           dwdx(k)*dudz(k) +                                &
+                           dwdy(k)*dvdz(k) )
 
           IF ( def < 0.0_wp )  def = 0.0_wp
@@ -3372,59 +3665,150 @@
           flag  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 
 
-          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
+          tend_temp(k) = km(k,j,i) * def * flag
 
        ENDDO
 
+!        ELSE
+! 
+!           DO  k = nzb+1, nzt
+! !
+! !--          Production term according to Kato and Launder (1993)
+!              def = SQRT( ( dudy(k)**2 + dvdz(k)**2 + dwdx(k)**2 +              &
+!                            dudz(k)**2 + dvdx(k)**2 + dwdy(k)**2 +              &
+!                            2.0_wp * ( dudy(k) * dvdx(k) +                      &
+!                                       dvdz(k) * dwdy(k) +                      &
+!                                       dwdx(k) * dudz(k) )       )              &
+!                        * ( dudy(k)**2 + dvdz(k)**2 + dwdx(k)**2 +              &
+!                            dudz(k)**2 + dvdx(k)**2 + dwdy(k)**2 -              &
+!                            2.0_wp * ( dudy(k) * dvdx(k) +                      &
+!                                       dvdz(k) * dwdy(k) +                      &
+!                                       dwdx(k) * dudz(k) )       )              &
+!                        )
+! 
+!              IF ( def < 0.0_wp )  def = 0.0_wp
+! 
+!              flag  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 
+! 
+!              tend_temp(k) = km(k,j,i) * def * flag
+! 
+!           ENDDO
+! 
+!        ENDIF
+
+    ELSE  ! not constant_flux_layer
+
+!        IF ( .NOT. rans_tke_e )  THEN
+!
+!--       Calculate TKE production by shear. Here, no additional
+!--       wall-bounded code is considered. 
+!--       Why?
+          DO  k = nzb+1, nzt
+
+             dudx(k)  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+             dudy(k)  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -                &
+                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+             dudz(k)  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -                &
+                                    u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+             dvdx(k)  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -                &
+                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+             dvdy(k)  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+             dvdz(k)  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -                &
+                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+             dwdx(k)  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -                &
+                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+             dwdy(k)  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -                &
+                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+             dwdz(k)  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+             def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) +         &
+                              dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2   +         &
+                              dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2   +         &
+                   2.0_wp * ( dvdx(k)*dudy(k) +                                &
+                              dwdx(k)*dudz(k) +                                &
+                              dwdy(k)*dvdz(k) )
+
+             IF ( def < 0.0_wp )  def = 0.0_wp
+
+             flag        = MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(k,j,i), 29 ) ) 
+             tend_temp(k) = km(k,j,i) * def * flag
+
+          ENDDO
+
+!        ELSE
+! 
+!           DO  k = nzb+1, nzt
+! 
+!              dudx(k)  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+!              dudy(k)  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -                &
+!                                     u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+!              dudz(k)  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -                &
+!                                     u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+! 
+!              dvdx(k)  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -                &
+!                                     v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+!              dvdy(k)  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+!              dvdz(k)  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -                &
+!                                     v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+! 
+!              dwdx(k)  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -                &
+!                                     w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+!              dwdy(k)  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -                &
+!                                     w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+!              dwdz(k)  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+! !
+! !--          Production term according to Kato and Launder (1993)
+!              def = SQRT( ( dudy(k)**2 + dvdz(k)**2 + dwdx(k)**2 +              &
+!                            dudz(k)**2 + dvdx(k)**2 + dwdy(k)**2 +              &
+!                            2.0_wp * ( dudy(k) * dvdx(k) +                      &
+!                                       dvdz(k) * dwdy(k) +                      &
+!                                       dwdx(k) * dudz(k) )       )              &
+!                        * ( dudy(k)**2 + dvdz(k)**2 + dwdx(k)**2 +              &
+!                            dudz(k)**2 + dvdx(k)**2 + dwdy(k)**2 -              &
+!                            2.0_wp * ( dudy(k) * dvdx(k) +                      &
+!                                       dvdz(k) * dwdy(k) +                      &
+!                                       dwdx(k) * dudz(k) )       )              &
+!                        )
+! 
+!              IF ( def < 0.0_wp )  def = 0.0_wp
+! 
+!              flag        = MERGE( 1.0_wp, 0.0_wp,                              &
+!                                   BTEST( wall_flags_0(k,j,i), 29 ) ) 
+!              tend_temp(k) = km(k,j,i) * def * flag
+! 
+!           ENDDO
+! 
+!        ENDIF
+
+    ENDIF
+
+    IF ( .NOT. diss_production )  THEN
+!
+!--    Production term in case of TKE production
+       DO  k = nzb+1, nzt
+          tend(k,j,i) = tend(k,j,i) + tend_temp(k)
+       ENDDO
     ELSE
 !
-!--    Calculate TKE production by shear. Here, no additional
-!--    wall-bounded code is considered. 
-!--    Why?
+!--    Production term in case of dissipation-rate production (rans_tke_e)
        DO  k = nzb+1, nzt
 
-          dudx(k)  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-          dudy(k)  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -                   &
-                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-          dudz(k)  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -                   &
-                                 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-          dvdx(k)  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -                   &
-                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-          dvdy(k)  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-          dvdz(k)  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -                   &
-                                 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-          dwdx(k)  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -                   &
-                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-          dwdy(k)  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -                   &
-                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-          dwdz(k)  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-          def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) +            &
-                           dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2   +            &
-                           dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2   +            &
-                2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + dwdy(k)*dvdz(k) )
-
-!
-!--       Production term according to Kato and Launder (1993)
-!           def = SQRT( ( dudx(k) + dudy(k) + dudz(k) +                        &
-!                         dvdx(k) + dvdy(k) + dvdz(k) +                        &
-!                         dwdx(k) + dwdy(k) + dwdz(k)                          &
-!                       )**4 -                                                 &
-!                       ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 +               &
-!                         2.0_wp * ( dudy(k) * dvdx(k) +                       &
-!                                    dudz(k) * dwdx(k) +                       &
-!                                    dvdz(k) * dwdy(k) )                       &
-!                       )**2                                                   &
-!                     )
-
-          IF ( def < 0.0_wp )  def = 0.0_wp
-
-          flag        = MERGE( 1.0_wp, 0.0_wp,                                 &
-                               BTEST( wall_flags_0(k,j,i), 29 ) ) 
-          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
-
+          ! Standard TKE-e closure
+          tend(k,j,i) = tend(k,j,i) + tend_temp(k) * diss(k,j,i)               &
+                                    /( e(k,j,i) + 1.0E-20_wp )                 &
+                                    * c_1
+!           ! Production according to Koblitz (2013)
+!           tend(k,j,i) = tend(k,j,i) + tend_temp(k) * diss(k,j,i)               &
+!                                     /( e(k,j,i) + 1.0E-20_wp )                 &
+!                                     * ( c_1 + ( c_2 - c_1 )                    &
+!                                             * l_wall(k,j,i) / l_max )
+!           ! Production according to Detering and Etling (1985)
+!           !> @todo us is not correct if there are vertical walls
+!           tend(k,j,i) = tend(k,j,i) + tend_temp(k) * SQRT(e(k,j,i))            &
+!                                     * c_1 * c_0**3 / c_4 * f                   &
+!              / surf_def_h(0)%us(surf_def_h(0)%start_index(j,i))
        ENDDO
-
     ENDIF
 
@@ -3822,10 +4206,10 @@
     REAL(wp)     ::  l              !< mixing length
     REAL(wp)     ::  ll             !< adjusted l
-    REAL(wp)     ::  var_reference  !< 
+    REAL(wp)     ::  var_reference  !< reference temperature
 
 #if defined( __nopointer )
-    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !< 
+    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !< temperature
 #else
-    REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< 
+    REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< temperature
 #endif
     REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dissipation  !< TKE dissipation
@@ -3837,4 +4221,7 @@
        DO  j = nys, nyn
           DO  k = nzb+1, nzt
+!
+!--          Predetermine flag to mask topography
+             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
 
 !
@@ -3851,5 +4238,7 @@
                 CALL mixing_length_rans( i, j, k, l, ll, var, var_reference )
 
-                dissipation(k,j) = c_m**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll
+                dissipation(k,j) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll
+
+                diss(k,j,i) = dissipation(k,j) * flag
 
              ELSEIF ( rans_tke_e )  THEN
@@ -3859,23 +4248,20 @@
              ENDIF
 
-!
-!--          Predetermine flag to mask topography
-             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-
-             tend(k,j,i) = tend(k,j,i)                                         &
-                                     + (                                       &
+             tend(k,j,i) = tend(k,j,i) + (                                     &
+                                           (                                   &
                        ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )     &
                      - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )     &
-                                       ) * ddx2  * flag                        &
-                                     + (                                       &
+                                           ) * ddx2  * flag                    &
+                                         + (                                   &
                        ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )     &
                      - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )     &
-                                       ) * ddy2  * flag                        &
-                                     + (                                       &
+                                           ) * ddy2  * flag                    &
+                                         + (                                   &
             ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1)    &
                                                           * rho_air_zw(k)      &
           - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)      &
                                                           * rho_air_zw(k-1)    &
-                                       ) * ddzw(k) * drho_air(k) * flag        &
+                                           ) * ddzw(k) * drho_air(k)           &
+                                         ) * flag * dsig_e                     &
                           - dissipation(k,j) * flag
 
@@ -3959,10 +4345,10 @@
     REAL(wp)     ::  l              !< mixing length
     REAL(wp)     ::  ll             !< adjusted l
-    REAL(wp)     ::  var_reference  !< 
+    REAL(wp)     ::  var_reference  !< reference temperature
 
 #if defined( __nopointer )
-    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !<
+    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !< temperature
 #else
-    REAL(wp), DIMENSION(:,:,:), POINTER ::  var     !<
+    REAL(wp), DIMENSION(:,:,:), POINTER ::  var     !< temperature
 #endif
     REAL(wp), DIMENSION(nzb+1:nzt) ::  dissipation  !< dissipation of TKE
@@ -3991,5 +4377,7 @@
           CALL mixing_length_rans( i, j, k, l, ll, var, var_reference  )
 
-          dissipation(k) = c_m**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll
+          dissipation(k) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll
+
+          diss(k,j,i) = dissipation(k) * flag
 
        ELSEIF ( rans_tke_e )  THEN
@@ -4001,13 +4389,13 @@
 !
 !--    Calculate the tendency term
-       tend(k,j,i) = tend(k,j,i)                                               &
-                                    + (                                        &
+       tend(k,j,i) = tend(k,j,i) + (                                           &
+                                      (                                        &
                       ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )      &
                     - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )      &
-                                      ) * ddx2  * flag / sig_e                 &
+                                      ) * ddx2                                 &
                                     + (                                        &
                       ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )      &
                     - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )      &
-                                      ) * ddy2  * flag / sig_e                 &
+                                      ) * ddy2                                 &
                                     + (                                        &
            ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1)     &
@@ -4015,6 +4403,7 @@
          - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)       &
                                                          * rho_air_zw(k-1)     &
-                                      ) * ddzw(k) * drho_air(k) * flag / sig_e &
-                                    - dissipation(k) * flag
+                                      ) * ddzw(k) * drho_air(k)                &
+                                   ) * flag * dsig_e                           &
+                                 - dissipation(k) * flag
 
     ENDDO
@@ -4082,13 +4471,13 @@
              flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
 
-             tend(k,j,i) = tend(k,j,i)                                         &
-                               + (                                             &
+             tend(k,j,i) = tend(k,j,i) +                                       &
+                         (       (                                             &
                  ( km(k,j,i)+km(k,j,i+1) ) * ( diss(k,j,i+1)-diss(k,j,i) )     &
                - ( km(k,j,i)+km(k,j,i-1) ) * ( diss(k,j,i)-diss(k,j,i-1) )     &
-                                 ) * ddx2  * flag                              &
+                                 ) * ddx2                                      &
                                + (                                             &
                  ( km(k,j,i)+km(k,j+1,i) ) * ( diss(k,j+1,i)-diss(k,j,i) )     &
                - ( km(k,j,i)+km(k,j-1,i) ) * ( diss(k,j,i)-diss(k,j-1,i) )     &
-                                 ) * ddy2  * flag                              &
+                                 ) * ddy2                                      &
                                + (                                             &
       ( km(k,j,i)+km(k+1,j,i) ) * ( diss(k+1,j,i)-diss(k,j,i) ) * ddzu(k+1)    &
@@ -4096,7 +4485,8 @@
     - ( km(k,j,i)+km(k-1,j,i) ) * ( diss(k,j,i)-diss(k-1,j,i) ) * ddzu(k)      &
                                                     * rho_air_zw(k-1)          &
-                                 ) * ddzw(k) * drho_air(k) * flag              &
-                               - c_2 * diss(k,j,i)**2                          &
-                                     / ( e(k,j,i) + 1.0E-20_wp ) * flag
+                                 ) * ddzw(k) * drho_air(k)                     &
+                         ) * flag * dsig_diss                                  &
+                         - c_2 * diss(k,j,i)**2                                &
+                               / ( e(k,j,i) + 1.0E-20_wp ) * flag
 
           ENDDO
@@ -4123,11 +4513,9 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  i              !< running index x direction
-    INTEGER(iwp) ::  j              !< running index y direction
-    INTEGER(iwp) ::  k              !< running index z direction
-
-    REAL(wp)     ::  flag           !< flag to mask topography
-
-    REAL(wp), DIMENSION(nzb+1:nzt) :: tend_temp
+    INTEGER(iwp) ::  i         !< running index x direction
+    INTEGER(iwp) ::  j         !< running index y direction
+    INTEGER(iwp) ::  k         !< running index z direction
+
+    REAL(wp)     ::  flag      !< flag to mask topography
 
 !
@@ -4141,12 +4529,13 @@
 !
 !--    Calculate the tendency term
-       tend_temp(k) =           (                                              &
+       tend(k,j,i) =  tend(k,j,i) +                                            &
+                   (            (                                              &
                 ( km(k,j,i)+km(k,j,i+1) ) * ( diss(k,j,i+1)-diss(k,j,i) )      &
               - ( km(k,j,i)+km(k,j,i-1) ) * ( diss(k,j,i)-diss(k,j,i-1) )      &
-                                ) * ddx2  * flag / sig_diss                    &
+                                ) * ddx2                                       &
                               + (                                              &
                 ( km(k,j,i)+km(k,j+1,i) ) * ( diss(k,j+1,i)-diss(k,j,i) )      &
               - ( km(k,j,i)+km(k,j-1,i) ) * ( diss(k,j,i)-diss(k,j-1,i) )      &
-                                ) * ddy2  * flag / sig_diss                    &
+                                ) * ddy2                                       &
                               + (                                              &
      ( km(k,j,i)+km(k+1,j,i) ) * ( diss(k+1,j,i)-diss(k,j,i) ) * ddzu(k+1)     &
@@ -4154,9 +4543,7 @@
    - ( km(k,j,i)+km(k-1,j,i) ) * ( diss(k,j,i)-diss(k-1,j,i) ) * ddzu(k)       &
                                                    * rho_air_zw(k-1)           &
-                                ) * ddzw(k) * drho_air(k) * flag / sig_diss    &
-                              - c_2 * diss(k,j,i)**2                           &
-                                    / ( e(k,j,i) + 1.0E-20_wp ) * flag
-
-       tend(k,j,i) = tend(k,j,i) + tend_temp(k)
+                                ) * ddzw(k) * drho_air(k)                      &
+                   ) * flag * dsig_diss                                        &
+                   - c_2 * diss(k,j,i)**2 / ( e(k,j,i) + 1.0E-20_wp ) * flag
 
     ENDDO
@@ -4290,4 +4677,5 @@
 !> Computation of the turbulent diffusion coefficients for momentum and heat 
 !> according to Prandtl-Kolmogorov.
+!> @todo consider non-default surfaces
 !------------------------------------------------------------------------------!
  SUBROUTINE tcm_diffusivities( var, var_reference )
@@ -4297,30 +4685,33 @@
         ONLY:  e_min, outflow_l, outflow_n, outflow_r, outflow_s
 
+    USE grid_variables,                                                        &
+        ONLY:  dx, dy
+
     USE statistics,                                                            &
         ONLY :  rmask, sums_l_l
 
     USE surface_mod,                                                           &
-        ONLY :  bc_h, surf_def_h
+        ONLY :  bc_h, surf_def_h, surf_def_v
 
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  i                   !<
-    INTEGER(iwp) ::  j                   !<
-    INTEGER(iwp) ::  k                   !<
-    INTEGER(iwp) ::  m                   !<
-    INTEGER(iwp) ::  n                   !<
-    INTEGER(iwp) ::  omp_get_thread_num  !<
-    INTEGER(iwp) ::  sr                  !<
-    INTEGER(iwp) ::  tn                  !<
-
-    REAL(wp)     ::  flag                !<
-    REAL(wp)     ::  l                   !<
-    REAL(wp)     ::  ll                  !<
-    REAL(wp)     ::  var_reference       !<
+    INTEGER(iwp) ::  i                   !< loop index
+    INTEGER(iwp) ::  j                   !< loop index
+    INTEGER(iwp) ::  k                   !< loop index
+    INTEGER(iwp) ::  m                   !< loop index
+    INTEGER(iwp) ::  n                   !< loop index
+    INTEGER(iwp) ::  omp_get_thread_num  !< opemmp function to get thread number
+    INTEGER(iwp) ::  sr                  !< statistic region
+    INTEGER(iwp) ::  tn                  !< thread number
+
+    REAL(wp)     ::  flag                !< topography flag
+    REAL(wp)     ::  l                   !< mixing length
+    REAL(wp)     ::  ll                  !< adjusted mixing length
+    REAL(wp)     ::  var_reference       !< reference temperature
 
 #if defined( __nopointer )
-    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !< 
+    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !< temperature
 #else
-    REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< 
+    REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< temperature
 #endif
 
@@ -4365,5 +4756,5 @@
 !
 !--             Compute diffusion coefficients for momentum and heat
-                km(k,j,i) = c_m * l * SQRT( e(k,j,i) ) * flag
+                km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag
                 kh(k,j,i) = ( 1.0_wp + 2.0_wp * l / ll ) * km(k,j,i) * flag
 !
@@ -4391,5 +4782,5 @@
 !
 !--             Compute diffusion coefficients for momentum and heat
-                km(k,j,i) = c_m * l * SQRT( e(k,j,i) ) * flag
+                km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag
                 kh(k,j,i) = km(k,j,i) / prandtl_number * flag
 !
@@ -4414,6 +4805,13 @@
 !
 !--             Compute diffusion coefficients for momentum and heat
-                km(k,j,i) = c_mu * e(k,j,i)**2 / ( diss(k,j,i) + 1.E-10 ) * flag
+                km(k,j,i) = c_0**4 * e(k,j,i)**2 / ( diss(k,j,i) + 1.0E-30_wp ) * flag
                 kh(k,j,i) = km(k,j,i) / prandtl_number * flag
+!
+!--             Summation for averaged profile of mixing length (cf. flow_statistics)
+                DO  sr = 0, statistic_regions
+                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) +                     &
+                      c_0**3 * e(k,j,i) * SQRT(e(k,j,i)) /                     &
+                      ( diss(k,j,i) + 1.0E-30_wp ) * rmask(j,i,sr) * flag
+                ENDDO
 
              ENDDO
@@ -4434,5 +4832,4 @@
 !-- so far vertical surfaces require usage of a Prandtl-layer where the boundary
 !-- values of the diffusivities are not needed.
-
     IF ( .NOT. rans_tke_e )  THEN
 !
@@ -4468,4 +4865,26 @@
           ENDDO
        ENDDO
+!
+!--    North- and southward facing surfaces
+       DO  n = 0, 1
+          DO  m = 1, surf_def_v(n)%ns
+             i = surf_def_v(n)%i(m)
+             j = surf_def_v(n)%j(m)
+             k = surf_def_v(n)%k(m)
+             km(k,j,i) = kappa * surf_def_v(n)%us(m) * 0.5_wp * dy
+             kh(k,j,i) = 1.35_wp * km(k,j,i)
+          ENDDO
+       ENDDO
+!
+!--    West- and eastward facing surfaces
+       DO  n = 2, 3
+          DO  m = 1, surf_def_v(n)%ns
+             i = surf_def_v(n)%i(m)
+             j = surf_def_v(n)%j(m)
+             k = surf_def_v(n)%k(m)
+             km(k,j,i) = kappa * surf_def_v(n)%us(m) * 0.5_wp * dx
+             kh(k,j,i) = 1.35_wp * km(k,j,i)
+          ENDDO
+       ENDDO
 
        CALL exchange_horiz( km, nbgp )
@@ -4473,4 +4892,5 @@
 
     ENDIF
+
 !
 !-- Model top
@@ -4518,5 +4938,5 @@
     INTEGER(iwp) ::  j      !< loop index y direction
     INTEGER(iwp) ::  k      !< loop index z direction
-    INTEGER, INTENT(IN) ::  mod_count  !<
+    INTEGER, INTENT(IN) ::  mod_count  !< flag defining where pointers point to
 
 #if defined( __nopointer )
Index: palm/trunk/SOURCE/vertical_nesting_mod.f90
===================================================================
--- palm/trunk/SOURCE/vertical_nesting_mod.f90	(revision 3069)
+++ palm/trunk/SOURCE/vertical_nesting_mod.f90	(revision 3083)
@@ -4058,4 +4058,9 @@
 !-- Identical timestep for coarse and fine grids
           dt_3d = MIN( dtc, dtf )
+          !> @fixme setting old_dt might be obsolete at this point
+          !>   Due to changes in timestep routine, setting of old_dt might be
+          !>   not necessary any more at this point. However, could not be
+          !>   tested so far.
+          !>   2018-05-18, gronemeier
           old_dt = dt_3d
 #endif
