Index: /palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- /palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3397)
+++ /palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3398)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Refactored production_e and production_e_ij (removed code redundancy)
+! 
+! 3386 2018-10-19 16:28:22Z gronemeier
 ! Renamed tcm_prognostic to tcm_prognostic_equations
 ! 
@@ -227,5 +230,5 @@
     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_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
@@ -234,5 +237,5 @@
 
     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
+       (/ 0.55_wp, 1.44_wp, 1.92_wp, 1.44_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)
@@ -423,5 +426,5 @@
        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_3 = rans_const_c(3)   !> @todo clarify how to switch between different models
        c_4 = rans_const_c(4)
 
@@ -2094,5 +2097,5 @@
        IF ( rans_tke_e )  advec = tend
 
-       CALL production_e
+       CALL production_e( .FALSE. )
 
 !
@@ -2547,5 +2550,5 @@
 !> @todo Adjust production term in case of rans_tke_e simulation
 !------------------------------------------------------------------------------!
- SUBROUTINE production_e
+ SUBROUTINE production_e( diss_production )
 
     USE arrays_3d,                                                             &
@@ -2568,4 +2571,6 @@
 
     IMPLICIT NONE
+
+    LOGICAL :: diss_production
 
     INTEGER(iwp) ::  i       !< running index x-direction
@@ -2600,4 +2605,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) ::  tmp_flux  !< temporary flux-array in z-direction
 
 
@@ -2824,5 +2830,16 @@
              flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),flag_nr) )
 
-             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
+             IF ( .NOT. diss_production )  THEN
+
+!--             Compute temdency for TKE-production from shear
+                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
+
+             ELSE
+
+!--             RANS mode: Compute temdency for dissipation-rate-production from shear
+                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag *           &
+                              diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * c_1
+
+             ENDIF
 
           ENDDO
@@ -2844,17 +2861,7 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
+
                    DO  k = nzb+1, nzt
-                      tend(k,j,i) = tend(k,j,i) +                              &
-                                    kh(k,j,i) * g /                            &
-                           MERGE( rho_reference, prho(k,j,i),                  &
-                                  use_single_reference_value ) *               &
-                                    ( prho(k+1,j,i) - prho(k-1,j,i) ) *        &
-                                    dd2zu(k) *                                 &
-                                MERGE( 1.0_wp, 0.0_wp,                         &
-                                       BTEST( wall_flags_0(k,j,i), 30 )        &
-                                     )                            *            &
-                                MERGE( 1.0_wp, 0.0_wp,                         &
-                                       BTEST( wall_flags_0(k,j,i), 9 )         &
-                                     )
+                      tmp_flux(k) = kh(k,j,i) * ( prho(k+1,j,i) - prho(k-1,j,i) ) * dd2zu(k)
                    ENDDO
 !
@@ -2867,12 +2874,7 @@
                          DO  m = surf_s, surf_e
                             k = surf_def_h(l)%k(m)
-                            tend(k,j,i) = tend(k,j,i) + g /                    &
-                                      MERGE( rho_reference, prho(k,j,i),       &
-                                             use_single_reference_value ) *    &
-                                      drho_air_zw(k-1) *                       &
-                                      surf_def_h(l)%shf(m)
+                            tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
                          ENDDO
                       ENDDO
-
                    ENDIF
 
@@ -2882,12 +2884,32 @@
                       DO  m = surf_s, surf_e
                          k = surf_def_h(2)%k(m)
-                         tend(k,j,i) = tend(k,j,i) + g /                       &
-                                      MERGE( rho_reference, prho(k,j,i),       &
-                                             use_single_reference_value ) *    &
-                                      drho_air_zw(k) *                         &
-                                      surf_def_h(2)%shf(m)
+                         tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
                       ENDDO
                    ENDIF
 
+                   IF ( .NOT. diss_production )  THEN
+
+!--                   Compute temdency for TKE-production from shear
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
+                         tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
+                                       MERGE( rho_reference, prho(k,j,i),       &
+                                              use_single_reference_value ) )
+                      ENDDO
+
+                   ELSE
+
+!--                   RANS mode: Compute temdency for dissipation-rate-production from shear
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
+                         tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
+                                       MERGE( rho_reference, prho(k,j,i),       &
+                                              use_single_reference_value ) ) *  &
+                                       diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *  &
+                                       c_3
+                      ENDDO
+
+                   ENDIF
+
                 ENDDO
              ENDDO
@@ -2899,19 +2921,5 @@
 
                    DO  k = nzb+1, nzt
-!
-!--                   Flag 9 is used to mask top fluxes, flag 30 to mask
-!--                   surface fluxes
-                      tend(k,j,i) = tend(k,j,i) -                              &
-                                    kh(k,j,i) * g /                            &
-                                MERGE( pt_reference, pt(k,j,i),                &
-                                        use_single_reference_value ) *         &
-                                    ( pt(k+1,j,i) - pt(k-1,j,i) ) *            &
-                                    dd2zu(k)                      *            &
-                                MERGE( 1.0_wp, 0.0_wp,                         &
-                                       BTEST( wall_flags_0(k,j,i), 30 )        &
-                                     )                            *            &
-                                MERGE( 1.0_wp, 0.0_wp,                         &
-                                       BTEST( wall_flags_0(k,j,i), 9 )         &
-                                     )
+                      tmp_flux(k) = -1.0_wp * kh(k,j,i) * ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
                    ENDDO
 
@@ -2924,9 +2932,5 @@
                          DO  m = surf_s, surf_e
                             k = surf_def_h(l)%k(m)
-                            tend(k,j,i) = tend(k,j,i) + g /                    &
-                                 MERGE( pt_reference, pt(k,j,i),               &
-                                        use_single_reference_value )           &
-                                                   * drho_air_zw(k-1)          &
-                                                   * surf_def_h(l)%shf(m)
+                            tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
                          ENDDO
                       ENDDO
@@ -2937,9 +2941,5 @@
                       DO  m = surf_s, surf_e
                          k = surf_lsm_h%k(m)
-                         tend(k,j,i) = tend(k,j,i) + g /                       &
-                                 MERGE( pt_reference, pt(k,j,i),               &
-                                        use_single_reference_value )           &
-                                                   * drho_air_zw(k-1)          &
-                                                   * surf_lsm_h%shf(m)
+                         tmp_flux(k) = drho_air_zw(k-1) * surf_lsm_h%shf(m)
                       ENDDO
 !
@@ -2949,9 +2949,5 @@
                       DO  m = surf_s, surf_e
                          k = surf_usm_h%k(m)
-                         tend(k,j,i) = tend(k,j,i) + g /                       &
-                                 MERGE( pt_reference, pt(k,j,i),               &
-                                        use_single_reference_value )           &
-                                                   * drho_air_zw(k-1)          &
-                                                   * surf_usm_h%shf(m)
+                         tmp_flux(k) = drho_air_zw(k-1) * surf_usm_h%shf(m)
                       ENDDO
                    ENDIF
@@ -2962,12 +2958,32 @@
                       DO  m = surf_s, surf_e
                          k = surf_def_h(2)%k(m)
-                         tend(k,j,i) = tend(k,j,i) + g /                       &
-                                 MERGE( pt_reference, pt(k,j,i),               &
-                                        use_single_reference_value )           &
-                                                   * drho_air_zw(k)            &
-                                                   * surf_def_h(2)%shf(m)
+                         tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
                       ENDDO
                    ENDIF
 
+                   IF ( .NOT. diss_production )  THEN
+
+!--                   Compute temdency for TKE-production from shear
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
+                         tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
+                                       MERGE( pt_reference, pt(k,j,i),          &
+                                              use_single_reference_value ) )
+                      ENDDO
+
+                   ELSE
+
+!--                   RANS mode: Compute temdency for dissipation-rate-production from shear
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
+                         tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
+                                       MERGE( pt_reference, pt(k,j,i),          &
+                                              use_single_reference_value ) ) *  &
+                                       diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *  &
+                                       c_3
+                      ENDDO
+
+                   ENDIF
+
                 ENDDO
              ENDDO
@@ -2981,23 +2997,12 @@
 
                 DO  k = nzb+1, nzt
-!
-!--                Flag 9 is used to mask top fluxes, flag 30 to mask
-!--                surface fluxes
+
                    IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
                       k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                       k2 = 0.61_wp * pt(k,j,i)
-                      tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *                  &
-                                      g /                                      &
-                                 MERGE( vpt_reference, vpt(k,j,i),             &
-                                        use_single_reference_value ) *         &
-                                      ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +     &
-                                        k2 * ( q(k+1,j,i) - q(k-1,j,i) )       &
-                                      ) * dd2zu(k) *                           &
-                                   MERGE( 1.0_wp, 0.0_wp,                      &
-                                          BTEST( wall_flags_0(k,j,i), 30 )     &
-                                        )          *                           &
-                                   MERGE( 1.0_wp, 0.0_wp,                      &
-                                          BTEST( wall_flags_0(k,j,i), 9 )      &
-                                        )
+                      tmp_flux(k) = -1.0_wp * kh(k,j,i) *                      &
+                                      ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
+                                        k2 * ( q(k+1,j,i)  - q(k-1,j,i) )      &
+                                      ) * dd2zu(k)
                    ELSE IF ( bulk_cloud_model )  THEN
                       IF ( ql(k,j,i) == 0.0_wp )  THEN
@@ -3014,34 +3019,16 @@
                          k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
                       ENDIF
-                      tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *                  &
-                                      g /                                      &
-                                 MERGE( vpt_reference, vpt(k,j,i),             &
-                                        use_single_reference_value ) *         &
-                                      ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +     &
-                                        k2 * ( q(k+1,j,i) - q(k-1,j,i) )       &
-                                      ) * dd2zu(k) *                           &
-                                   MERGE( 1.0_wp, 0.0_wp,                      &
-                                          BTEST( wall_flags_0(k,j,i), 30 )     &
-                                        )          *                           &
-                                   MERGE( 1.0_wp, 0.0_wp,                      &
-                                          BTEST( wall_flags_0(k,j,i), 9 )      &
-                                        )
+                      tmp_flux(k) = -1.0_wp * kh(k,j,i) *                      &
+                                      ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
+                                        k2 * ( q(k+1,j,i)  - q(k-1,j,i) )      &
+                                      ) * dd2zu(k)
                    ELSE IF ( cloud_droplets )  THEN
                       k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
                       k2 = 0.61_wp * pt(k,j,i)
-                      tend(k,j,i) = tend(k,j,i) -                              &
-                                    kh(k,j,i) * g /                            &
-                                 MERGE( vpt_reference, vpt(k,j,i),             &
-                                        use_single_reference_value ) *         &
-                                    ( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) +      &
-                                      k2 * ( q(k+1,j,i) -  q(k-1,j,i) ) -      &
-                                      pt(k,j,i) * ( ql(k+1,j,i) -              &
-                                      ql(k-1,j,i) ) ) * dd2zu(k) *             &
-                                   MERGE( 1.0_wp, 0.0_wp,                      &
-                                          BTEST( wall_flags_0(k,j,i), 30 )     &
-                                        )                        *             &
-                                   MERGE( 1.0_wp, 0.0_wp,                      &
-                                          BTEST( wall_flags_0(k,j,i), 9 )      &
-                                        )
+                      tmp_flux(k) = -1.0_wp * kh(k,j,i) * &
+                                      ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
+                                        k2 * ( q(k+1,j,i)  - q(k-1,j,i) ) -    &
+                                        pt(k,j,i) * ( ql(k+1,j,i) -            &
+                                        ql(k-1,j,i) ) ) * dd2zu(k)
                    ENDIF
 
@@ -3079,10 +3066,7 @@
                          ENDIF
 
-                         tend(k,j,i) = tend(k,j,i) + g /                       &
-                                 MERGE( vpt_reference, vpt(k,j,i),             &
-                                        use_single_reference_value ) *         &
-                                            ( k1 * surf_def_h(l)%shf(m) +      &
-                                              k2 * surf_def_h(l)%qsws(m)       &
-                                            ) * drho_air_zw(k-1)
+                         tmp_flux(k) = ( k1 * surf_def_h(l)%shf(m) +           &
+                                         k2 * surf_def_h(l)%qsws(m)            &
+                                       ) * drho_air_zw(k-1)
                       ENDDO
                    ENDDO
@@ -3116,10 +3100,7 @@
                       ENDIF
 
-                      tend(k,j,i) = tend(k,j,i) + g /                          &
-                                 MERGE( vpt_reference, vpt(k,j,i),             &
-                                        use_single_reference_value ) *         &
-                                            ( k1 * surf_lsm_h%shf(m) +         &
-                                              k2 * surf_lsm_h%qsws(m)          &
-                                            ) * drho_air_zw(k-1)
+                      tmp_flux(k) = ( k1 * surf_lsm_h%shf(m) +                 &
+                                      k2 * surf_lsm_h%qsws(m)                  &
+                                    ) * drho_air_zw(k-1)
                    ENDDO
 !
@@ -3152,10 +3133,7 @@
                       ENDIF
 
-                      tend(k,j,i) = tend(k,j,i) + g /                          &
-                                 MERGE( vpt_reference, vpt(k,j,i),             &
-                                        use_single_reference_value ) *         &
-                                            ( k1 * surf_usm_h%shf(m) +         &
-                                              k2 * surf_usm_h%qsws(m)          &
-                                            ) * drho_air_zw(k-1)
+                      tmp_flux(k) = ( k1 * surf_usm_h%shf(m) +                 &
+                                      k2 * surf_usm_h%qsws(m)                  &
+                                    ) * drho_air_zw(k-1)
                    ENDDO
 
@@ -3191,14 +3169,35 @@
                       ENDIF
 
-                      tend(k,j,i) = tend(k,j,i) + g /                          &
-                                 MERGE( vpt_reference, vpt(k,j,i),             &
-                                        use_single_reference_value ) *         &
-                                            ( k1 * surf_def_h(2)%shf(m) +      &
-                                              k2 * surf_def_h(2)%qsws(m)       &
-                                            ) * drho_air_zw(k)
+                      tmp_flux(k) = ( k1 * surf_def_h(2)%shf(m) +              &
+                                      k2 * surf_def_h(2)%qsws(m)               &
+                                    ) * drho_air_zw(k)
 
                    ENDDO
 
                 ENDIF ! from IF ( use_top_fluxes )  THEN
+
+                IF ( .NOT. diss_production )  THEN
+
+!--                Compute temdency for TKE-production from shear
+                   DO  k = nzb+1, nzt
+                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
+                      tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
+                                    MERGE( vpt_reference, vpt(k,j,i),          &
+                                           use_single_reference_value ) )
+                   ENDDO
+
+                ELSE
+
+!--                RANS mode: Compute temdency for dissipation-rate-production from shear
+                   DO  k = nzb+1, nzt
+                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
+                      tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g /   &
+                                    MERGE( vpt_reference, vpt(k,j,i),          &
+                                           use_single_reference_value ) ) *    &
+                                    diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *    &
+                                    c_3
+                   ENDDO
+
+                ENDIF
 
              ENDDO
@@ -3227,5 +3226,5 @@
 
     USE control_parameters,                                                    &
-        ONLY:  cloud_droplets, constant_flux_layer, neutral,                &
+        ONLY:  cloud_droplets, constant_flux_layer, neutral,                   &
                rho_reference, use_single_reference_value, use_surface_fluxes,  &
                use_top_fluxes
@@ -3239,5 +3238,5 @@
     USE surface_mod,                                                           &
         ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
-                surf_usm_v 
+                surf_usm_v
 
     IMPLICIT NONE
@@ -3252,4 +3251,5 @@
     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) ::  flag_nr !< number of masking flag
 
     REAL(wp)     ::  def         !<
@@ -3260,5 +3260,5 @@
     REAL(wp)     ::  theta       !<
     REAL(wp)     ::  temp        !<
-    REAL(wp)     ::  sign_dir    !< sign of wall-tke flux, depending on wall orientation 
+    REAL(wp)     ::  sign_dir    !< sign of wall-tke flux, depending on wall orientation
     REAL(wp)     ::  usvs        !< momentum flux u"v"
     REAL(wp)     ::  vsus        !< momentum flux v"u"
@@ -3266,48 +3266,57 @@
     REAL(wp)     ::  wsvs        !< momentum flux w"v"
 
-
-    REAL(wp), DIMENSION(nzb+1:nzt)  ::  dudx        !< Gradient of u-component in x-direction
-    REAL(wp), DIMENSION(nzb+1:nzt)  ::  dudy        !< Gradient of u-component in y-direction
-    REAL(wp), DIMENSION(nzb+1:nzt)  ::  dudz        !< Gradient of u-component in z-direction
-    REAL(wp), DIMENSION(nzb+1:nzt)  ::  dvdx        !< Gradient of v-component in x-direction
-    REAL(wp), DIMENSION(nzb+1:nzt)  ::  dvdy        !< Gradient of v-component in y-direction
-    REAL(wp), DIMENSION(nzb+1:nzt)  ::  dvdz        !< Gradient of v-component in z-direction
-    REAL(wp), DIMENSION(nzb+1:nzt)  ::  dwdx        !< Gradient of w-component in x-direction
-    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
+    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudx  !< Gradient of u-component in x-direction
+    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudy  !< Gradient of u-component in y-direction
+    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudz  !< Gradient of u-component in z-direction
+    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdx  !< Gradient of v-component in x-direction
+    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdy  !< Gradient of v-component in y-direction
+    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdz  !< Gradient of v-component in z-direction
+    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdx  !< Gradient of w-component in x-direction
+    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) ::  tmp_flux  !< temporary flux-array in z-direction
+
+
+
+!
+!-- Calculate TKE production by shear. Calculate gradients at all grid
+!-- points first, gradients at surface-bounded grid points will be
+!-- overwritten further below.
+    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)
+
+    ENDDO
+
+    flag_nr = 29
 
     IF ( constant_flux_layer )  THEN
-!
-!--    Calculate TKE production by shear. Calculate gradients at all grid 
-!--    points first, gradients at surface-bounded grid points will be 
-!--    overwritten further below. 
-       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)
-
-       ENDDO
+
+       flag_nr = 0
+
+!--    Position beneath wall
+!--    (2) - Will allways be executed.
+!--    'bottom and wall: use u_0,v_0 and wall functions'
 !
 !--    Compute gradients at north- and south-facing surfaces.
-!--    Note, no vertical natural surfaces so far.
+!--    First, for default surfaces, then for urban surfaces.
+!--    Note, so far no natural vertical surfaces implemented
        DO  l = 0, 1
-!
-!--       Default surfaces
           surf_s = surf_def_v(l)%start_index(j,i)
           surf_e = surf_def_v(l)%end_index(j,i)
@@ -3317,13 +3326,13 @@
              wsvs        = surf_def_v(l)%mom_flux_tke(1,m)
 
-             km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp               &
+             km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
                              * 0.5_wp * dy
 !
 !--          -1.0 for right-facing wall, 1.0 for left-facing wall
-             sign_dir = MERGE( 1.0_wp, -1.0_wp,                                &
-                               BTEST( wall_flags_0(k,j-1,i), 0 ) ) 
+             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
+                               BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
              dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
              dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
-          ENDDO   
+          ENDDO
 !
 !--       Natural surfaces
@@ -3335,13 +3344,13 @@
              wsvs        = surf_lsm_v(l)%mom_flux_tke(1,m)
 
-             km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp               &
+             km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
                              * 0.5_wp * dy
 !
 !--          -1.0 for right-facing wall, 1.0 for left-facing wall
-             sign_dir = MERGE( 1.0_wp, -1.0_wp,                                &
-                               BTEST( wall_flags_0(k,j-1,i), 0 ) ) 
+             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
+                               BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
              dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
              dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
-          ENDDO 
+          ENDDO
 !
 !--       Urban surfaces
@@ -3353,75 +3362,72 @@
              wsvs        = surf_usm_v(l)%mom_flux_tke(1,m)
 
-             km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp               &
+             km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
                              * 0.5_wp * dy
 !
 !--          -1.0 for right-facing wall, 1.0 for left-facing wall
-             sign_dir = MERGE( 1.0_wp, -1.0_wp,                                &
-                               BTEST( wall_flags_0(k,j-1,i), 0 ) ) 
+             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
+                               BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
              dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
              dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
-          ENDDO 
+          ENDDO
        ENDDO
 !
-!--    Compute gradients at east- and west-facing walls
+!--          Compute gradients at east- and west-facing walls
        DO  l = 2, 3
-!
-!--       Default surfaces
           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)
-             vsus        = surf_def_v(l)%mom_flux_tke(0,m)
-             wsus        = surf_def_v(l)%mom_flux_tke(1,m)
-
-             km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp               &
+             k     = surf_def_v(l)%k(m)
+             vsus  = surf_def_v(l)%mom_flux_tke(0,m)
+             wsus  = surf_def_v(l)%mom_flux_tke(1,m)
+
+             km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
                                 * 0.5_wp * dx
 !
 !--          -1.0 for right-facing wall, 1.0 for left-facing wall
-             sign_dir = MERGE( 1.0_wp, -1.0_wp,                                &
-                               BTEST( wall_flags_0(k,j,i-1), 0 ) )
+             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
+                               BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
              dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
              dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
-          ENDDO 
-!
-!--       Natural surfaces
+          ENDDO
+!
+!--             Natural surfaces
           surf_s = surf_lsm_v(l)%start_index(j,i)
           surf_e = surf_lsm_v(l)%end_index(j,i)
           DO  m = surf_s, surf_e
-             k           = surf_lsm_v(l)%k(m)
-             vsus        = surf_lsm_v(l)%mom_flux_tke(0,m)
-             wsus        = surf_lsm_v(l)%mom_flux_tke(1,m)
-
-             km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp               &
+             k     = surf_lsm_v(l)%k(m)
+             vsus  = surf_lsm_v(l)%mom_flux_tke(0,m)
+             wsus  = surf_lsm_v(l)%mom_flux_tke(1,m)
+
+             km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
                                 * 0.5_wp * dx
 !
 !--          -1.0 for right-facing wall, 1.0 for left-facing wall
-             sign_dir = MERGE( 1.0_wp, -1.0_wp,                                &
-                               BTEST( wall_flags_0(k,j,i-1), 0 ) ) 
+             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
+                               BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
              dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
              dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
-          ENDDO 
-!
-!--       Urban surfaces
+          ENDDO
+!
+!--             Urban surfaces
           surf_s = surf_usm_v(l)%start_index(j,i)
           surf_e = surf_usm_v(l)%end_index(j,i)
           DO  m = surf_s, surf_e
-             k           = surf_usm_v(l)%k(m)
-             vsus        = surf_usm_v(l)%mom_flux_tke(0,m)
-             wsus        = surf_usm_v(l)%mom_flux_tke(1,m)
-
-             km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp               &
+             k     = surf_usm_v(l)%k(m)
+             vsus  = surf_usm_v(l)%mom_flux_tke(0,m)
+             wsus  = surf_usm_v(l)%mom_flux_tke(1,m)
+
+             km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
                                 * 0.5_wp * dx
 !
 !--          -1.0 for right-facing wall, 1.0 for left-facing wall
-             sign_dir = MERGE( 1.0_wp, -1.0_wp,                                &
-                               BTEST( wall_flags_0(k,j,i-1), 0 ) ) 
+             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
+                               BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
              dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
              dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
-          ENDDO 
+          ENDDO
        ENDDO
 !
-!--    Compute gradients at upward-facing walls, first for 
-!--    non-natural default surfaces
+!--          Compute gradients at upward-facing surfaces
        surf_s = surf_def_h(0)%start_index(j,i)
        surf_e = surf_def_h(0)%end_index(j,i)
@@ -3429,11 +3435,11 @@
           k = surf_def_h(0)%k(m)
 !
-!--       Please note, actually, an interpolation of u_0 and v_0 
-!--       onto the grid center would be required. However, this 
+!--       Please note, actually, an interpolation of u_0 and v_0
+!--       onto the grid center would be required. However, this
 !--       would require several data transfers between 2D-grid and
-!--       wall type. The effect of this missing interpolation is 
+!--       wall type. The effect of this missing interpolation is
 !--       negligible. (See also production_e_init).
-          dudz(k)     = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k)   
-          dvdz(k)     = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k)
+          dudz(k) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k)
+          dvdz(k) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k)
 
        ENDDO
@@ -3445,6 +3451,7 @@
           k = surf_lsm_h%k(m)
 
-          dudz(k)     = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k)   
-          dvdz(k)     = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k)
+          dudz(k) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k)
+          dvdz(k) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k)
+
        ENDDO
 !
@@ -3455,9 +3462,10 @@
           k = surf_usm_h%k(m)
 
-          dudz(k)     = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k)   
-          dvdz(k)     = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k)
+          dudz(k) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k)
+          dvdz(k) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k)
+
        ENDDO
 !
-!--    Compute gradients at downward-facing walls, only for 
+!--    Compute gradients at downward-facing walls, only for
 !--    non-natural default surfaces
        surf_s = surf_def_h(1)%start_index(j,i)
@@ -3466,171 +3474,38 @@
           k = surf_def_h(1)%k(m)
 
-          dudz(k)     = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k)   
-          dvdz(k)     = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k) 
+          dudz(k) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k)
+          dvdz(k) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k)
 
        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) )
-
-          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
-
-!        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
-!
-!--    Production term in case of dissipation-rate production (rans_tke_e)
-       DO  k = nzb+1, nzt
-
-          ! 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
+    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) )
+
+       IF ( def < 0.0_wp )  def = 0.0_wp
+
+       flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),flag_nr) )
+
+       IF ( .NOT. diss_production )  THEN
+
+!--       Compute temdency for TKE-production from shear
+          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
+
+       ELSE
+
+!--       RANS mode: Compute temdency for dissipation-rate-production from shear
+          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag *           &
+                        diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * c_1
+
+       ENDIF
+
+    ENDDO
+
 
 !
@@ -3642,20 +3517,62 @@
           IF ( ocean_mode )  THEN
 !
-!--          So far in the ocean no special treatment of density flux in
-!--          the bottom and top surface layer
+!--          So far in the ocean no special treatment of density flux
+!--          in the bottom and top surface layer
+
              DO  k = nzb+1, nzt
- 
-                tend(k,j,i) = tend(k,j,i) +                                    &
-                              kh(k,j,i) * g /                                  &
-                              MERGE( rho_reference, prho(k,j,i),               &
-                                     use_single_reference_value ) *            &
-                              ( prho(k+1,j,i) - prho(k-1,j,i) ) *              &
-                              dd2zu(k) *                                       &
-                                MERGE( 1.0_wp, 0.0_wp,                         &
-                                       BTEST( wall_flags_0(k,j,i), 30 )        &
-                                     ) *                                       &
-                                MERGE( 1.0_wp, 0.0_wp,                         &
-                                       BTEST( wall_flags_0(k,j,i), 9 )         &
-                                     )           
+                tmp_flux(k) = kh(k,j,i) * ( prho(k+1,j,i) - prho(k-1,j,i) ) * dd2zu(k)
+             ENDDO
+!
+!--          Treatment of near-surface grid points, at up- and down-
+!--          ward facing surfaces
+             IF ( use_surface_fluxes )  THEN
+                DO  l = 0, 1
+                   surf_s = surf_def_h(l)%start_index(j,i)
+                   surf_e = surf_def_h(l)%end_index(j,i)
+                   DO  m = surf_s, surf_e
+                      k = surf_def_h(l)%k(m)
+                      tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
+                   ENDDO
+                ENDDO
+             ENDIF
+
+             IF ( use_top_fluxes )  THEN
+                surf_s = surf_def_h(2)%start_index(j,i)
+                surf_e = surf_def_h(2)%end_index(j,i)
+                DO  m = surf_s, surf_e
+                   k = surf_def_h(2)%k(m)
+                   tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
+                ENDDO
+             ENDIF
+
+             IF ( .NOT. diss_production )  THEN
+
+!--             Compute temdency for TKE-production from shear
+                DO  k = nzb+1, nzt
+                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
+                   tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
+                                 MERGE( rho_reference, prho(k,j,i),       &
+                                        use_single_reference_value ) )
+                ENDDO
+
+             ELSE
+
+!--             RANS mode: Compute temdency for dissipation-rate-production from shear
+                DO  k = nzb+1, nzt
+                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
+                   tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
+                                 MERGE( rho_reference, prho(k,j,i),       &
+                                        use_single_reference_value ) ) *  &
+                                 diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *  &
+                                 c_3
+                ENDDO
+
+             ENDIF
+
+
+          ELSE ! or IF ( .NOT. ocean_mode )  THEN
+
+             DO  k = nzb+1, nzt
+                tmp_flux(k) = -1.0_wp * kh(k,j,i) * ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
              ENDDO
 
@@ -3668,12 +3585,23 @@
                    DO  m = surf_s, surf_e
                       k = surf_def_h(l)%k(m)
-                      tend(k,j,i) = tend(k,j,i) + g /                          &
-                                MERGE( rho_reference, prho(k,j,i),             &
-                                       use_single_reference_value ) *          &
-                                drho_air_zw(k-1) *                             &
-                                surf_def_h(l)%shf(m)
+                      tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
                    ENDDO
                 ENDDO
-
+!
+!--             Natural surfaces
+                surf_s = surf_lsm_h%start_index(j,i)
+                surf_e = surf_lsm_h%end_index(j,i)
+                DO  m = surf_s, surf_e
+                   k = surf_lsm_h%k(m)
+                   tmp_flux(k) = drho_air_zw(k-1) * surf_lsm_h%shf(m)
+                ENDDO
+!
+!--             Urban surfaces
+                surf_s = surf_usm_h%start_index(j,i)
+                surf_e = surf_usm_h%end_index(j,i)
+                DO  m = surf_s, surf_e
+                   k = surf_usm_h%k(m)
+                   tmp_flux(k) = drho_air_zw(k-1) * surf_usm_h%shf(m)
+                ENDDO
              ENDIF
 
@@ -3683,109 +3611,45 @@
                 DO  m = surf_s, surf_e
                    k = surf_def_h(2)%k(m)
-                   tend(k,j,i) = tend(k,j,i) + g /                             &
-                                MERGE( rho_reference, prho(k,j,i),             &
-                                       use_single_reference_value ) *          &
-                                drho_air_zw(k) *                               &
-                                surf_def_h(2)%shf(m) 
+                   tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
                 ENDDO
              ENDIF
 
-          ELSE
-
-             DO  k = nzb+1, nzt
-!
-!--             Flag 9 is used to mask top fluxes, flag 30 to mask
-!--             surface fluxes
-                tend(k,j,i) = tend(k,j,i) -                                    &
-                              kh(k,j,i) * g /                                  &
-                                MERGE( pt_reference, pt(k,j,i),                &
-                                       use_single_reference_value ) *          &
-                              ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) *       &
-                                MERGE( 1.0_wp, 0.0_wp,                         &
-                                       BTEST( wall_flags_0(k,j,i), 30 )        &
-                                     ) *                                       &
-                                MERGE( 1.0_wp, 0.0_wp,                         &
-                                       BTEST( wall_flags_0(k,j,i), 9 )         &
-                                     )
-
-             ENDDO
-
-             IF ( use_surface_fluxes )  THEN
-!
-!--             Default surfaces, up- and downward-facing
-                DO  l = 0, 1
-                   surf_s = surf_def_h(l)%start_index(j,i)
-                   surf_e = surf_def_h(l)%end_index(j,i)
-                   DO  m = surf_s, surf_e
-                      k = surf_def_h(l)%k(m)
-                      tend(k,j,i) = tend(k,j,i) + g /                          &
-                                MERGE( pt_reference, pt(k,j,i),                &
-                                       use_single_reference_value ) *          &
-                                drho_air_zw(k-1) *                             &
-                                surf_def_h(l)%shf(m)
-                   ENDDO
+             IF ( .NOT. diss_production )  THEN
+
+!--             Compute temdency for TKE-production from shear
+                DO  k = nzb+1, nzt
+                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
+                   tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
+                                 MERGE( pt_reference, pt(k,j,i),          &
+                                        use_single_reference_value ) )
                 ENDDO
-!
-!--             Natural surfaces
-                surf_s = surf_lsm_h%start_index(j,i)
-                surf_e = surf_lsm_h%end_index(j,i)
-                DO  m = surf_s, surf_e
-                   k = surf_lsm_h%k(m)
-                   tend(k,j,i) = tend(k,j,i) + g /                             &
-                                MERGE( pt_reference, pt(k,j,i),                &
-                                       use_single_reference_value ) *          &
-                                drho_air_zw(k-1) *                             &
-                                surf_lsm_h%shf(m) 
+
+             ELSE
+
+!--             RANS mode: Compute temdency for dissipation-rate-production from shear
+                DO  k = nzb+1, nzt
+                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
+                   tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
+                                 MERGE( pt_reference, pt(k,j,i),          &
+                                        use_single_reference_value ) ) *  &
+                                 diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *  &
+                                 c_3
                 ENDDO
-!
-!--             Urban surfaces
-                surf_s = surf_usm_h%start_index(j,i)
-                surf_e = surf_usm_h%end_index(j,i)
-                DO  m = surf_s, surf_e
-                   k = surf_usm_h%k(m)
-                   tend(k,j,i) = tend(k,j,i) + g /                             &
-                                MERGE( pt_reference, pt(k,j,i),                &
-                                       use_single_reference_value ) *          &
-                                drho_air_zw(k-1) *                             &
-                                surf_usm_h%shf(m) 
-                ENDDO
+
              ENDIF
 
-             IF ( use_top_fluxes )  THEN
-                surf_s = surf_def_h(2)%start_index(j,i)
-                surf_e = surf_def_h(2)%end_index(j,i)
-                DO  m = surf_s, surf_e
-                   k = surf_def_h(2)%k(m)
-                   tend(k,j,i) = tend(k,j,i) + g /                             &
-                                MERGE( pt_reference, pt(k,j,i),                &
-                                       use_single_reference_value ) *          &
-                                drho_air_zw(k) *                               &
-                                surf_def_h(2)%shf(m) 
-                ENDDO
-             ENDIF
-
-          ENDIF
-
-       ELSE
+          ENDIF ! from IF ( .NOT. ocean_mode )
+
+       ELSE ! or IF ( humidity )  THEN
 
           DO  k = nzb+1, nzt
-!
-!--          Flag 9 is used to mask top fluxes, flag 30 to mask surface fluxes
-             IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets )  THEN
+
+             IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
                 k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                 k2 = 0.61_wp * pt(k,j,i)
-                tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g /                    &
-                                MERGE( vpt_reference, vpt(k,j,i),              &
-                                       use_single_reference_value ) *          &
-                                      ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +     &
-                                        k2 * ( q(k+1,j,i) - q(k-1,j,i) )       &
-                                      ) * dd2zu(k) *                           &
-                                   MERGE( 1.0_wp, 0.0_wp,                      &
-                                          BTEST( wall_flags_0(k,j,i), 30 )     &
-                                        )          *                           &
-                                   MERGE( 1.0_wp, 0.0_wp,                      &
-                                          BTEST( wall_flags_0(k,j,i), 9 )      &
-                                        )
-
+                tmp_flux(k) = -1.0_wp * kh(k,j,i) *                      &
+                                ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
+                                  k2 * ( q(k+1,j,i)  - q(k-1,j,i) )      &
+                                ) * dd2zu(k)
              ELSE IF ( bulk_cloud_model )  THEN
                 IF ( ql(k,j,i) == 0.0_wp )  THEN
@@ -3795,45 +3659,30 @@
                    theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
                    temp  = theta * exner(k)
-                   k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                        &
-                                 ( q(k,j,i) - ql(k,j,i) ) *                    &
-                        ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /              &
-                        ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *               &
+                   k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                  &
+                                 ( q(k,j,i) - ql(k,j,i) ) *              &
+                        ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /        &
+                        ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *         &
                         ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
                    k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
                 ENDIF
-                tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g /                    &
-                                MERGE( vpt_reference, vpt(k,j,i),              &
-                                       use_single_reference_value ) *          &
-                                      ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +     &
-                                        k2 * ( q(k+1,j,i) - q(k-1,j,i) )       &
-                                      ) * dd2zu(k) *                           &
-                                   MERGE( 1.0_wp, 0.0_wp,                      &
-                                          BTEST( wall_flags_0(k,j,i), 30 )     &
-                                        )          *                           &
-                                   MERGE( 1.0_wp, 0.0_wp,                      &
-                                          BTEST( wall_flags_0(k,j,i), 9 )      &
-                                        )
+                tmp_flux(k) = -1.0_wp * kh(k,j,i) *                      &
+                                ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
+                                  k2 * ( q(k+1,j,i)  - q(k-1,j,i) )      &
+                                ) * dd2zu(k)
              ELSE IF ( cloud_droplets )  THEN
                 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
                 k2 = 0.61_wp * pt(k,j,i)
-                tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g /                    &
-                                MERGE( vpt_reference, vpt(k,j,i),              &
-                                       use_single_reference_value ) *          &
-                                  ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +         &
-                                    k2 * ( q(k+1,j,i) - q(k-1,j,i) ) -         &
-                                    pt(k,j,i) * ( ql(k+1,j,i) -                &
-                                                  ql(k-1,j,i) ) ) * dd2zu(k)   &
-                                 * MERGE( 1.0_wp, 0.0_wp,                      &
-                                          BTEST( wall_flags_0(k,j,i), 30 )     &
-                                        )                                      &
-                                 * MERGE( 1.0_wp, 0.0_wp,                      &
-                                          BTEST( wall_flags_0(k,j,i), 9 )      &
-                                        )    
+                tmp_flux(k) = -1.0_wp * kh(k,j,i) * &
+                                ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
+                                  k2 * ( q(k+1,j,i)  - q(k-1,j,i) ) -    &
+                                  pt(k,j,i) * ( ql(k+1,j,i) -            &
+                                  ql(k-1,j,i) ) ) * dd2zu(k)
              ENDIF
+
           ENDDO
 
           IF ( use_surface_fluxes )  THEN
 !
-!--          Treat horizontal default surfaces, up- and downward-facing
+!--          Treat horizontal default surfaces
              DO  l = 0, 1
                 surf_s = surf_def_h(l)%start_index(j,i)
@@ -3850,12 +3699,12 @@
                          k2 = 0.61_wp * pt(k,j,i)
                       ELSE
-                        theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
-                        temp  = theta * exner(k)
-                        k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                   &
-                                   ( q(k,j,i) - ql(k,j,i) ) *                  &
-                          ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /            &
-                          ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *             &
-                          ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                        k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
+                         theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
+                         temp  = theta * exner(k)
+                         k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *            &
+                                    ( q(k,j,i) - ql(k,j,i) ) *           &
+                           ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /     &
+                           ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *      &
+                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+                         k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
                       ENDIF
                    ELSE IF ( cloud_droplets )  THEN
@@ -3864,10 +3713,7 @@
                    ENDIF
 
-                   tend(k,j,i) = tend(k,j,i) + g /                             &
-                                MERGE( vpt_reference, vpt(k,j,i),              &
-                                       use_single_reference_value ) *          &
-                                      ( k1 * surf_def_h(l)%shf(m) +            &
-                                        k2 * surf_def_h(l)%qsws(m)             &
-                                      ) * drho_air_zw(k-1)
+                   tmp_flux(k) = ( k1 * surf_def_h(l)%shf(m) +           &
+                                   k2 * surf_def_h(l)%qsws(m)            &
+                                 ) * drho_air_zw(k-1)
                 ENDDO
              ENDDO
@@ -3880,80 +3726,4 @@
 
                 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
-                    k1 = 1.0_wp + 0.61_wp * q(k,j,i)
-                    k2 = 0.61_wp * pt(k,j,i)
-                ELSE IF ( bulk_cloud_model )  THEN
-                    IF ( ql(k,j,i) == 0.0_wp )  THEN
-                       k1 = 1.0_wp + 0.61_wp * q(k,j,i)
-                       k2 = 0.61_wp * pt(k,j,i)
-                    ELSE
-                       theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
-                       temp  = theta * exner(k)
-                       k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                    &
-                                     ( q(k,j,i) - ql(k,j,i) ) *                &
-                            ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /          &
-                            ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *           &
-                            ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                       k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
-                   ENDIF
-                ELSE IF ( cloud_droplets )  THEN
-                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
-                   k2 = 0.61_wp * pt(k,j,i)
-                ENDIF
-
-                tend(k,j,i) = tend(k,j,i) + g /                                &
-                                MERGE( vpt_reference, vpt(k,j,i),              &
-                                       use_single_reference_value ) *          &
-                                         ( k1 * surf_lsm_h%shf(m) +            &
-                                           k2 * surf_lsm_h%qsws(m)             &
-                                         ) * drho_air_zw(k-1)
-             ENDDO
-!
-!--          Treat horizontal urban surfaces
-             surf_s = surf_usm_h%start_index(j,i)
-             surf_e = surf_usm_h%end_index(j,i)
-             DO  m = surf_s, surf_e
-                k = surf_usm_h%k(m)
-
-                IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
-                    k1 = 1.0_wp + 0.61_wp * q(k,j,i)
-                    k2 = 0.61_wp * pt(k,j,i)
-                ELSE IF ( bulk_cloud_model )  THEN
-                    IF ( ql(k,j,i) == 0.0_wp )  THEN
-                       k1 = 1.0_wp + 0.61_wp * q(k,j,i)
-                       k2 = 0.61_wp * pt(k,j,i)
-                    ELSE
-                       theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
-                       temp  = theta * exner(k)
-                       k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                    &
-                                     ( q(k,j,i) - ql(k,j,i) ) *                &
-                            ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /          &
-                            ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *           &
-                            ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                       k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
-                   ENDIF
-                ELSE IF ( cloud_droplets )  THEN
-                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
-                   k2 = 0.61_wp * pt(k,j,i)
-                ENDIF
-
-                tend(k,j,i) = tend(k,j,i) + g /                                &
-                                MERGE( vpt_reference, vpt(k,j,i),              &
-                                       use_single_reference_value ) *          &
-                                         ( k1 * surf_usm_h%shf(m) +            &
-                                           k2 * surf_usm_h%qsws(m)             &
-                                         ) * drho_air_zw(k-1)
-             ENDDO
-
-          ENDIF
-
-          IF ( use_top_fluxes )  THEN
-             surf_s = surf_def_h(2)%start_index(j,i)
-             surf_e = surf_def_h(2)%end_index(j,i)
-             DO  m = surf_s, surf_e
-                k = surf_def_h(2)%k(m)
-
-
-
-                IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets )  THEN
                    k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                    k2 = 0.61_wp * pt(k,j,i)
@@ -3965,8 +3735,77 @@
                       theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
                       temp  = theta * exner(k)
-                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                     &
-                                 ( q(k,j,i) - ql(k,j,i) ) *                    &
-                        ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /              &
-                        ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *               &
+                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
+                                    ( q(k,j,i) - ql(k,j,i) ) *           &
+                           ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /     &
+                           ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *      &
+                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+                      k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
+                   ENDIF
+                ELSE IF ( cloud_droplets )  THEN
+                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+                   k2 = 0.61_wp * pt(k,j,i)
+                ENDIF
+
+                tmp_flux(k) = ( k1 * surf_lsm_h%shf(m) +                 &
+                                k2 * surf_lsm_h%qsws(m)                  &
+                              ) * drho_air_zw(k-1)
+             ENDDO
+!
+!--          Treat horizontal urban surfaces
+             surf_s = surf_usm_h%start_index(j,i)
+             surf_e = surf_usm_h%end_index(j,i)
+             DO  m = surf_s, surf_e
+                k = surf_usm_h%k(m)
+
+                IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
+                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                   k2 = 0.61_wp * pt(k,j,i)
+                ELSE IF ( bulk_cloud_model )  THEN
+                   IF ( ql(k,j,i) == 0.0_wp )  THEN
+                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                      k2 = 0.61_wp * pt(k,j,i)
+                   ELSE
+                      theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
+                      temp  = theta * exner(k)
+                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
+                                    ( q(k,j,i) - ql(k,j,i) ) *           &
+                           ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /     &
+                           ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *      &
+                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+                      k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
+                   ENDIF
+                ELSE IF ( cloud_droplets )  THEN
+                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+                   k2 = 0.61_wp * pt(k,j,i)
+                ENDIF
+
+                tmp_flux(k) = ( k1 * surf_usm_h%shf(m) +                 &
+                                k2 * surf_usm_h%qsws(m)                  &
+                              ) * drho_air_zw(k-1)
+             ENDDO
+
+          ENDIF ! from IF ( use_surface_fluxes )  THEN
+
+          IF ( use_top_fluxes )  THEN
+
+             surf_s = surf_def_h(2)%start_index(j,i)
+             surf_e = surf_def_h(2)%end_index(j,i)
+             DO  m = surf_s, surf_e
+                k = surf_def_h(2)%k(m)
+
+                IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
+                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                   k2 = 0.61_wp * pt(k,j,i)
+                ELSE IF ( bulk_cloud_model )  THEN
+                   IF ( ql(k,j,i) == 0.0_wp )  THEN
+                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                      k2 = 0.61_wp * pt(k,j,i)
+                   ELSE
+                      theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
+                      temp  = theta * exner(k)
+                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
+                                 ( q(k,j,i) - ql(k,j,i) ) *              &
+                        ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /        &
+                        ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *         &
                         ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
                       k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
@@ -3977,10 +3816,32 @@
                 ENDIF
 
-                tend(k,j,i) = tend(k,j,i) + g /                                &
-                                MERGE( vpt_reference, vpt(k,j,i),              &
-                                       use_single_reference_value ) *          &
-                            ( k1* surf_def_h(2)%shf(m) +                       &
-                              k2 * surf_def_h(2)%qsws(m)                       &
-                            ) * drho_air_zw(k)
+                tmp_flux(k) = ( k1 * surf_def_h(2)%shf(m) +              &
+                                k2 * surf_def_h(2)%qsws(m)               &
+                              ) * drho_air_zw(k)
+
+             ENDDO
+
+          ENDIF ! from IF ( use_top_fluxes )  THEN
+
+          IF ( .NOT. diss_production )  THEN
+
+!--          Compute temdency for TKE-production from shear
+             DO  k = nzb+1, nzt
+                flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
+                tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
+                              MERGE( vpt_reference, vpt(k,j,i),          &
+                                     use_single_reference_value ) )
+             ENDDO
+
+          ELSE
+
+!--          RANS mode: Compute temdency for dissipation-rate-production from shear
+             DO  k = nzb+1, nzt
+                flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
+                tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g /   &
+                              MERGE( vpt_reference, vpt(k,j,i),          &
+                                     use_single_reference_value ) ) *    &
+                              diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *    &
+                              c_3
              ENDDO
 
@@ -3991,5 +3852,5 @@
     ENDIF
 
-  END SUBROUTINE production_e_ij
+ END SUBROUTINE production_e_ij
 
 
