Index: palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3129)
+++ palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3130)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! - move boundary condition of km and kh to tcm_diffusivities
+! - calculate km at boundaries according to MOST
+! - move phi_m to surface_layer_fluxes_mod
+! 
+! 3129 2018-07-16 07:45:13Z gronemeier
 ! - move limitation of diss to boundary_conds
 ! - move boundary conditions for e and diss to boundary_conds
@@ -1910,6 +1915,4 @@
     INTEGER(iwp) ::  m      !< running index surface elements
 
-    REAL(wp)     :: km_sfc  !< km according to MOST
-
     IF ( constant_flux_layer )  THEN
 !
@@ -1922,6 +1925,4 @@
 !--    surface winds).
 !--    not available in case of non-cyclic boundary conditions.
-!--    WARNING: the exact analytical solution would require the determination
-!--             of the eddy diffusivity by km = u* * kappa * zp / phi_m.
 !--    Default surfaces, upward-facing
        !$OMP PARALLEL DO PRIVATE(i,j,k,m)
@@ -1932,19 +1933,16 @@
           k = surf_def_h(0)%k(m)
 !
-!--       Note, calculatione of u_0 and v_0 is not fully accurate, as u/v 
+!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v 
 !--       and km are not on the same grid. Actually, a further 
 !--       interpolation of km onto the u/v-grid is necessary. However, the
 !--       effect of this error is negligible.
-          km_sfc = surf_def_h(0)%us(m) * kappa * surf_def_h(0)%z_mo(m) /       &
-                   phi_m( surf_def_h(0)%z_mo(m) / surf_def_h(0)%ol(m) )
-
           surf_def_h(0)%u_0(m) = u(k+1,j,i) + surf_def_h(0)%usws(m) *          &
                                      drho_air_zw(k-1)               *          &
                                      ( zu(k+1) - zu(k-1)    )       /          &
-                                     ( km_sfc  + 1.0E-20_wp )
+                                     ( km(k,j,i)  + 1.0E-20_wp )
           surf_def_h(0)%v_0(m) = v(k+1,j,i) + surf_def_h(0)%vsws(m) *          &
                                      drho_air_zw(k-1)               *          &
                                      ( zu(k+1) - zu(k-1)    )       /          &
-                                     ( km_sfc  + 1.0E-20_wp )   
+                                     ( km(k,j,i)  + 1.0E-20_wp )   
 
           IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) )  >                     &
@@ -1965,5 +1963,9 @@
           j = surf_def_h(1)%j(m)
           k = surf_def_h(1)%k(m)
-
+!
+!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v 
+!--       and km are not on the same grid. Actually, a further 
+!--       interpolation of km onto the u/v-grid is necessary. However, the
+!--       effect of this error is negligible.
           surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) *          &
                                      drho_air_zw(k-1) *                        &
@@ -1989,24 +1991,20 @@
        DO  m = 1, surf_lsm_h%ns
 
-          i = surf_lsm_h%i(m)            
+          i = surf_lsm_h%i(m)
           j = surf_lsm_h%j(m)
           k = surf_lsm_h%k(m)
 !
-!--       Note, calculatione of u_0 and v_0 is not fully accurate, as u/v 
+!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v 
 !--       and km are not on the same grid. Actually, a further 
 !--       interpolation of km onto the u/v-grid is necessary. However, the
-!--       effect of this error is negligible. 
-!--       effect of this error is negligible. 
-          km_sfc = surf_lsm_h%us(m) * kappa * surf_lsm_h%z_mo(m) /             &
-                   phi_m( surf_lsm_h%z_mo(m) / surf_lsm_h%ol(m) )
-
+!--       effect of this error is negligible.
           surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m)    *             &
                                         drho_air_zw(k-1)         *             &
-                                        ( zu(k+1) - zu(k-1)    ) /             &
-                                        ( km_sfc  + 1.0E-20_wp )  
+                                        ( zu(k+1)   - zu(k-1)    ) /           &
+                                        ( km(k,j,i) + 1.0E-20_wp )  
           surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m)    *             &
                                         drho_air_zw(k-1)         *             &
                                         ( zu(k+1) - zu(k-1)    ) /             &
-                                        ( km_sfc  + 1.0E-20_wp )
+                                        ( km(k,j,i)  + 1.0E-20_wp )
 
           IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) )  >                        &
@@ -2024,23 +2022,20 @@
        DO  m = 1, surf_usm_h%ns
 
-          i = surf_usm_h%i(m)            
+          i = surf_usm_h%i(m)
           j = surf_usm_h%j(m)
           k = surf_usm_h%k(m)
 !
-!--       Note, calculatione of u_0 and v_0 is not fully accurate, as u/v 
+!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v 
 !--       and km are not on the same grid. Actually, a further 
 !--       interpolation of km onto the u/v-grid is necessary. However, the
-!--       effect of this error is negligible. 
-          km_sfc = surf_usm_h%us(m) * kappa * surf_usm_h%z_mo(m) /             &
-                   phi_m( surf_usm_h%z_mo(m) / surf_usm_h%ol(m) )
-
+!--       effect of this error is negligible.
           surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m)    *             &
                                         drho_air_zw(k-1)         *             &
                                         ( zu(k+1) - zu(k-1)    ) /             &
-                                        ( km_sfc  + 1.0E-20_wp )
+                                        ( km(k,j,i)  + 1.0E-20_wp )
           surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m)    *             &
                                         drho_air_zw(k-1)         *             &
                                         ( zu(k+1) - zu(k-1)    ) /             &
-                                        ( km_sfc  + 1.0E-20_wp )
+                                        ( km(k,j,i)  + 1.0E-20_wp )
 
           IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) )  >                        &
@@ -2057,31 +2052,4 @@
 
  END SUBROUTINE production_e_init
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Calculates stability function for momentum
-!>
-!> @author Hauke Wurps
-!------------------------------------------------------------------------------!
- FUNCTION phi_m( zeta ) 
-
-    IMPLICIT NONE 
-
-    REAL(wp)            :: phi_m         !< Value of the function
-    REAL(wp)            :: zeta          !< Stability parameter z/L
-
-    REAL(wp), PARAMETER :: a = 16.0_wp   !< constant
-    REAL(wp), PARAMETER :: b = -0.25_wp  !< constant
-    REAL(wp), PARAMETER :: c = 5.0_wp    !< constant
-
-    IF ( zeta < 0.0_wp )  THEN
-       phi_m = ( 1.0_wp - a * zeta )**b
-    ELSE
-       phi_m = 1.0_wp + c * zeta
-    ENDIF
-
- END FUNCTION phi_m
 
 
@@ -4626,4 +4594,20 @@
  SUBROUTINE tcm_diffusivities( var, var_reference )
 
+    USE control_parameters,                                                    &
+        ONLY:  e_min, outflow_l, outflow_n, outflow_r, outflow_s
+
+    USE surface_layer_fluxes_mod,                                              &
+        ONLY:  phi_m
+
+    USE surface_mod,                                                           &
+        ONLY :  bc_h, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v,          &
+                surf_usm_h, surf_usm_v
+
+    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
+
     REAL(wp) ::  var_reference       !< reference temperature
 
@@ -4633,77 +4617,9 @@
     REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< temperature
 #endif
-!
-!-- Call default diffusivities routine. This is always used to calculate kh.
-    CALL tcm_diffusivities_default( var, var_reference )
-!
-!-- Call dynamic subgrid model to calculate km.
-    IF ( les_dynamic )  THEN
-       CALL tcm_diffusivities_dynamic
-    ENDIF
-
- END SUBROUTINE tcm_diffusivities
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Computation of the turbulent diffusion coefficients for momentum and heat 
-!> according to Prandtl-Kolmogorov.
-!------------------------------------------------------------------------------!
- SUBROUTINE tcm_diffusivities_default( var, var_reference )
- 
-
-    USE control_parameters,                                                    &
-        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, surf_def_v, surf_lsm_h, surf_lsm_v,          &
-                surf_usm_h, surf_usm_v
-
-    IMPLICIT NONE
-
-    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  !< temperature
-#else
-    REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< temperature
-#endif
-
-!
-!-- Default thread number in case of one thread
-    tn = 0
-
-!
-!-- Initialization for calculation of the mixing length profile
-    sums_l_l = 0.0_wp
-
-!
-!-- Compute the turbulent diffusion coefficient for momentum
-    !$OMP PARALLEL PRIVATE (i,j,k,l,ll,sr,tn,flag)
-!$  tn = omp_get_thread_num()
 
 !
 !-- Introduce an optional minimum tke
     IF ( e_min > 0.0_wp )  THEN
-       !$OMP DO
+       !$OMP PARALLEL DO PRIVATE(i,j,k)
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
@@ -4716,183 +4632,102 @@
     ENDIF
 
-    IF ( les_dynamic .OR. les_mw )  THEN
-       !$OMP DO
-       DO  i = nxlg, nxrg
-          DO  j = nysg, nyng
-             DO  k = nzb+1, nzt
-
-                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-
-!
-!--             Determine the mixing length for LES closure
-                CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
-!
-!--             Compute diffusion coefficients for momentum and heat
-                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
-!
-!--             Summation for averaged profile (cf. flow_statistics)
-                DO  sr = 0, statistic_regions
-                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr)   &
-                                                             * flag
-                ENDDO
-
-             ENDDO
-          ENDDO
-       ENDDO
-
-    ELSEIF ( rans_tke_l )  THEN
-
-       !$OMP DO
-       DO  i = nxlg, nxrg
-          DO  j = nysg, nyng
-             DO  k = nzb+1, nzt
-
-                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-!
-!--             Mixing length for RANS mode with TKE-l closure
-                CALL mixing_length_rans( i, j, k, l, ll, var, var_reference )
-!
-!--             Compute diffusion coefficients for momentum and heat
-                km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag
-                kh(k,j,i) = km(k,j,i) / prandtl_number * flag
-!
-!--             Summation for averaged profile (cf. flow_statistics)
-                DO  sr = 0, statistic_regions
-                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr)   &
-                                                             * flag
-                ENDDO
-
-             ENDDO
-          ENDDO
-       ENDDO
-
-    ELSEIF ( rans_tke_e )  THEN
-
-       !$OMP DO
-       DO  i = nxlg, nxrg
-          DO  j = nysg, nyng
-             DO  k = nzb+1, nzt
-
-                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-!
-!--             Compute diffusion coefficients for momentum and heat
-                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
-          ENDDO
-       ENDDO
-
+!
+!-- Call default diffusivities routine. This is always used to calculate kh.
+    CALL tcm_diffusivities_default( var, var_reference )
+!
+!-- Call dynamic subgrid model to calculate km.
+    IF ( les_dynamic )  THEN
+       CALL tcm_diffusivities_dynamic
     ENDIF
 
-    sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn)   ! quasi boundary-condition for
-                                                ! data output
-!$OMP END PARALLEL
-
-!
-!-- Set vertical boundary values (Neumann conditions both at upward- and 
-!-- downward facing walls. To set wall-boundary values, the surface data type 
-!-- is applied.
-!-- Horizontal boundary conditions at vertical walls are not set because
-!-- 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
-!
-!--    Upward-facing
-       !$OMP PARALLEL DO PRIVATE( i, j, k )
-       DO  m = 1, bc_h(0)%ns
-          i = bc_h(0)%i(m)            
-          j = bc_h(0)%j(m)
-          k = bc_h(0)%k(m)
-          km(k-1,j,i) = km(k,j,i)
-          kh(k-1,j,i) = kh(k,j,i)
-       ENDDO
-!
-!--    Downward facing surfaces
-       !$OMP PARALLEL DO PRIVATE( i, j, k )
-       DO  m = 1, bc_h(1)%ns
-          i = bc_h(1)%i(m)            
-          j = bc_h(1)%j(m)
-          k = bc_h(1)%k(m)
-          km(k+1,j,i) = km(k,j,i)
-          kh(k+1,j,i) = kh(k,j,i)
-       ENDDO
-    ELSE
-!
-!--    Up- and downward facing surfaces
+!
+!-- Set boundary values (Neumann conditions at downward facing surfaces,
+!-- according to MOST at upward facing surfaces and vertical surfaces).
+!
+!-- Downward facing surfaces
+    !$OMP PARALLEL DO PRIVATE(i,j,k)
+    DO  m = 1, bc_h(1)%ns
+       i = bc_h(1)%i(m)
+       j = bc_h(1)%j(m)
+       k = bc_h(1)%k(m)
+       km(k+1,j,i) = km(k,j,i)
+       kh(k+1,j,i) = kh(k,j,i)
+    ENDDO
+
+!-- Upward facing surfaces
+!-- Default surfaces
+    n = 0
+    !$OMP PARALLEL DO PRIVATE(i,j,k,m)
+    DO  m = 1, surf_def_h(n)%ns
+       i = surf_def_h(n)%i(m)
+       j = surf_def_h(n)%j(m)
+       k = surf_def_h(n)%k(m)
+       km(k,j,i) = kappa * surf_def_h(n)%us(m) * surf_def_h(n)%z_mo(m) /       &
+                   phi_m( surf_def_h(n)%z_mo(m) / surf_def_h(n)%ol(m) )
+       kh(k,j,i) = 1.35_wp * km(k,j,i)
+    ENDDO
+!
+!-- Natural surfaces
+    !$OMP PARALLEL DO PRIVATE(i,j,k,m)
+    DO  m = 1, surf_lsm_h%ns
+       i = surf_lsm_h%i(m)
+       j = surf_lsm_h%j(m)
+       k = surf_lsm_h%k(m)
+       km(k,j,i) = kappa * surf_lsm_h%us(m) * surf_lsm_h%z_mo(m) /             &
+                   phi_m( surf_lsm_h%z_mo(m) / surf_lsm_h%ol(m) )
+       kh(k,j,i) = 1.35_wp * km(k,j,i)
+    ENDDO
+!
+!-- Urban surfaces
+    !$OMP PARALLEL DO PRIVATE(i,j,k,m)
+    DO  m = 1, surf_usm_h%ns
+       i = surf_usm_h%i(m)
+       j = surf_usm_h%j(m)
+       k = surf_usm_h%k(m)
+       km(k,j,i) = kappa * surf_usm_h%us(m) * surf_usm_h%z_mo(m) /             &
+                   phi_m( surf_usm_h%z_mo(m) / surf_usm_h%ol(m) )
+       kh(k,j,i) = 1.35_wp * km(k,j,i)
+    ENDDO
+
+!
+!-- North-, south-, west and eastward facing surfaces
+    DO  n = 0, 3
+!
 !--    Default surfaces
-       DO  n = 0, 1
-          DO  m = 1, surf_def_h(n)%ns
-             i = surf_def_h(n)%i(m)
-             j = surf_def_h(n)%j(m)
-             k = surf_def_h(n)%k(m)
-             km(k,j,i) = kappa * surf_def_h(n)%us(m) * surf_def_h(n)%z_mo(m)
-             kh(k,j,i) = 1.35_wp * km(k,j,i)
-          ENDDO
-       ENDDO
-!
-!--    Upward facing surfaces
-!--    Natural surfaces
-       DO  m = 1, surf_lsm_h%ns
-          i = surf_lsm_h%i(m)
-          j = surf_lsm_h%j(m)
-          k = surf_lsm_h%k(m)
-          km(k,j,i) = kappa * surf_lsm_h%us(m) * surf_lsm_h%z_mo(m)
+       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
+       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) * surf_def_v(n)%z_mo(m) /    &
+                      phi_m( surf_def_v(n)%z_mo(m) / surf_def_v(n)%ol(m) )
           kh(k,j,i) = 1.35_wp * km(k,j,i)
        ENDDO
 !
-!--    Urban surfaces
-       DO  m = 1, surf_usm_h%ns
-          i = surf_usm_h%i(m)
-          j = surf_usm_h%j(m)
-          k = surf_usm_h%k(m)
-          km(k,j,i) = kappa * surf_usm_h%us(m) * surf_usm_h%z_mo(m)
+!--    Natural surfaces
+       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
+       DO  m = 1, surf_lsm_v(n)%ns
+          i = surf_lsm_v(n)%i(m)
+          j = surf_lsm_v(n)%j(m)
+          k = surf_lsm_v(n)%k(m)
+          km(k,j,i) = kappa * surf_lsm_v(n)%us(m) * surf_lsm_v(n)%z_mo(m) /    &
+                      phi_m( surf_lsm_v(n)%z_mo(m) / surf_lsm_v(n)%ol(m) )
           kh(k,j,i) = 1.35_wp * km(k,j,i)
        ENDDO
-
-!
-!--    North-, south-, west and eastward facing surfaces
-       DO  n = 0, 3
-!
-!--       Default surfaces
-          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) * surf_def_v(n)%z_mo(m)
-             kh(k,j,i) = 1.35_wp * km(k,j,i)
-          ENDDO
-!
-!--       Natural surfaces
-          DO  m = 1, surf_lsm_v(n)%ns
-             i = surf_lsm_v(n)%i(m)
-             j = surf_lsm_v(n)%j(m)
-             k = surf_lsm_v(n)%k(m)
-             km(k,j,i) = kappa * surf_lsm_v(n)%us(m) * surf_lsm_v(n)%z_mo(m)
-             kh(k,j,i) = 1.35_wp * km(k,j,i)
-          ENDDO
-!
-!--       Urban surfaces
-          DO  m = 1, surf_usm_v(n)%ns
-             i = surf_usm_v(n)%i(m)
-             j = surf_usm_v(n)%j(m)
-             k = surf_usm_v(n)%k(m)
-             km(k,j,i) = kappa * surf_usm_v(n)%us(m) * surf_usm_v(n)%z_mo(m)
-             kh(k,j,i) = 1.35_wp * km(k,j,i)
-          ENDDO
+!
+!--    Urban surfaces
+       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
+       DO  m = 1, surf_usm_v(n)%ns
+          i = surf_usm_v(n)%i(m)
+          j = surf_usm_v(n)%j(m)
+          k = surf_usm_v(n)%k(m)
+          km(k,j,i) = kappa * surf_usm_v(n)%us(m) * surf_usm_v(n)%z_mo(m) /    &
+                      phi_m( surf_usm_v(n)%z_mo(m) / surf_usm_v(n)%ol(m) )
+          kh(k,j,i) = 1.35_wp * km(k,j,i)
        ENDDO
-
-       CALL exchange_horiz( km, nbgp )
-       CALL exchange_horiz( kh, nbgp )
-
-    ENDIF
+    ENDDO
+
+    CALL exchange_horiz( km, nbgp )
+    CALL exchange_horiz( kh, nbgp )
 
 !
@@ -4926,4 +4761,133 @@
     ENDIF
 
+ END SUBROUTINE tcm_diffusivities
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Computation of the turbulent diffusion coefficients for momentum and heat 
+!> according to Prandtl-Kolmogorov.
+!------------------------------------------------------------------------------!
+ SUBROUTINE tcm_diffusivities_default( var, var_reference )
+ 
+    USE statistics,                                                            &
+        ONLY :  rmask, sums_l_l
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  i                   !< loop index
+    INTEGER(iwp) ::  j                   !< loop index
+    INTEGER(iwp) ::  k                   !< 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  !< temperature
+#else
+    REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< temperature
+#endif
+
+!
+!-- Default thread number in case of one thread
+    tn = 0
+
+!
+!-- Initialization for calculation of the mixing length profile
+    sums_l_l = 0.0_wp
+
+!
+!-- Compute the turbulent diffusion coefficient for momentum
+    !$OMP PARALLEL PRIVATE (i,j,k,l,ll,sr,tn,flag)
+!$  tn = omp_get_thread_num()
+
+    IF ( les_dynamic .OR. les_mw )  THEN
+       !$OMP DO
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
+             DO  k = nzb+1, nzt
+
+                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+
+!
+!--             Determine the mixing length for LES closure
+                CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
+!
+!--             Compute diffusion coefficients for momentum and heat
+                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
+!
+!--             Summation for averaged profile (cf. flow_statistics)
+                DO  sr = 0, statistic_regions
+                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr)   &
+                                                             * flag
+                ENDDO
+
+             ENDDO
+          ENDDO
+       ENDDO
+
+    ELSEIF ( rans_tke_l )  THEN
+
+       !$OMP DO
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
+             DO  k = nzb+1, nzt
+
+                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+!
+!--             Mixing length for RANS mode with TKE-l closure
+                CALL mixing_length_rans( i, j, k, l, ll, var, var_reference )
+!
+!--             Compute diffusion coefficients for momentum and heat
+                km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag
+                kh(k,j,i) = km(k,j,i) / prandtl_number * flag
+!
+!--             Summation for averaged profile (cf. flow_statistics)
+                DO  sr = 0, statistic_regions
+                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr)   &
+                                                             * flag
+                ENDDO
+
+             ENDDO
+          ENDDO
+       ENDDO
+
+    ELSEIF ( rans_tke_e )  THEN
+
+       !$OMP DO
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
+             DO  k = nzb+1, nzt
+
+                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+!
+!--             Compute diffusion coefficients for momentum and heat
+                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
+          ENDDO
+       ENDDO
+
+    ENDIF
+
+    sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn)   ! quasi boundary-condition for
+                                                ! data output
+!$OMP END PARALLEL
+
  END SUBROUTINE tcm_diffusivities_default
 
@@ -4949,13 +4913,7 @@
     USE arrays_3d,                                                             &
         ONLY:  ddzw, dzw, dd2zu, w, ug, vg
-
-    USE control_parameters,                                                    &
-        ONLY:  e_min, outflow_l, outflow_n, outflow_r, outflow_s
     
     USE grid_variables,                                                        &
         ONLY : ddx, ddy, dx, dy
-
-    USE surface_mod,                                                           &
-        ONLY :  bc_h
 
     IMPLICIT NONE
@@ -4966,5 +4924,4 @@
     INTEGER(iwp) ::  l           !< running index
     INTEGER(iwp) ::  m           !< running index
-    INTEGER(iwp) ::  n           !< running index
 
     REAL(wp)     ::  dudx        !< Gradient of u-component in x-direction
@@ -5037,17 +4994,4 @@
     REAL(wp), PARAMETER :: fac_cmax = 23.0_wp/(24.0_wp*sqrt(3.0_wp))  !< constant
 
-!
-!-- Introduce an optional minimum tke
-    IF ( e_min > 0.0_wp )  THEN
-       !$OMP DO
-       DO  i = nxlg, nxrg
-          DO  j = nysg, nyng
-             DO  k = nzb+1, nzt
-                e(k,j,i) = MAX( e(k,j,i), e_min ) *                            &
-                        MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-             ENDDO
-          ENDDO
-       ENDDO
-    ENDIF
 !
 !-- velocities on grid centers:
@@ -5181,53 +5125,4 @@
        ENDDO
     ENDDO
-!
-!-- Set vertical boundary values (Neumann conditions both at upward- and 
-!-- downward facing walls. To set wall-boundary values, the surface data type 
-!-- is applied.
-!-- Horizontal boundary conditions at vertical walls are not set because
-!-- so far vertical surfaces require usage of a constant-flux layer where the
-!-- boundary values of the diffusivities are not needed.
-
-!
-!-- Upward-facing
-    !$OMP PARALLEL DO PRIVATE( i, j, k )
-    DO  m = 1, bc_h(0)%ns
-        i = bc_h(0)%i(m)            
-        j = bc_h(0)%j(m)
-        k = bc_h(0)%k(m)
-        km(k-1,j,i) = km(k,j,i)
-    ENDDO
-!
-!-- Downward facing surfaces
-    !$OMP PARALLEL DO PRIVATE( i, j, k )
-    DO  m = 1, bc_h(1)%ns
-        i = bc_h(1)%i(m)            
-        j = bc_h(1)%j(m)
-        k = bc_h(1)%k(m)
-        km(k+1,j,i) = km(k,j,i)
-    ENDDO
-!
-!-- Model top
-    !$OMP PARALLEL DO
-    DO  i = nxlg, nxrg
-       DO  j = nysg, nyng
-          km(nzt+1,j,i) = km(nzt,j,i)
-       ENDDO
-    ENDDO
-!
-!-- Set Neumann boundary conditions at the outflow boundaries in case of
-!-- non-cyclic lateral boundaries
-    IF ( outflow_l )  THEN
-       km(:,:,nxl-1) = km(:,:,nxl)
-    ENDIF
-    IF ( outflow_r )  THEN
-       km(:,:,nxr+1) = km(:,:,nxr)
-    ENDIF
-    IF ( outflow_s )  THEN
-       km(:,nys-1,:) = km(:,nys,:)
-    ENDIF
-    IF ( outflow_n )  THEN
-       km(:,nyn+1,:) = km(:,nyn,:)
-    ENDIF
 
  END SUBROUTINE tcm_diffusivities_dynamic
