Index: palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3130)
+++ palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3145)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Use MOST for km only in RANS mode
+! 
+! 3130 2018-07-16 11:08:55Z gronemeier
 ! - move boundary condition of km and kh to tcm_diffusivities
 ! - calculate km at boundaries according to MOST
@@ -1905,4 +1908,7 @@
         ONLY:  constant_flux_layer
 
+    USE surface_layer_fluxes_mod,                                              &
+        ONLY:  phi_m
+
     USE surface_mod,                                                           &
         ONLY :  surf_def_h, surf_lsm_h, surf_usm_h 
@@ -1914,4 +1920,6 @@
     INTEGER(iwp) ::  k      !< grid index z-direction
     INTEGER(iwp) ::  m      !< running index surface elements
+    
+    REAL(wp) ::  km_sfc     !< diffusion coefficient, used to compute virtual velocities
 
     IF ( constant_flux_layer )  THEN
@@ -1937,12 +1945,15 @@
 !--       interpolation of km onto the u/v-grid is necessary. However, the
 !--       effect of this error is negligible.
+          km_sfc = kappa * surf_def_h(0)%us(m) * 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(k,j,i)  + 1.0E-20_wp )
+                                     ( km_sfc  + 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(k,j,i)  + 1.0E-20_wp )   
+                                     ( km_sfc  + 1.0E-20_wp )   
 
           IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) )  >                     &
@@ -1999,12 +2010,15 @@
 !--       interpolation of km onto the u/v-grid is necessary. However, the
 !--       effect of this error is negligible.
+          km_sfc = 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) )
+
           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(k,j,i) + 1.0E-20_wp )  
+                                        ( zu(k+1) - zu(k-1)    ) /             &
+                                        ( km_sfc  + 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(k,j,i)  + 1.0E-20_wp )
+                                        ( km_sfc  + 1.0E-20_wp )
 
           IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) )  >                        &
@@ -2030,12 +2044,15 @@
 !--       interpolation of km onto the u/v-grid is necessary. However, the
 !--       effect of this error is negligible.
+          km_sfc = 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) )
+
           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(k,j,i)  + 1.0E-20_wp )
+                                        ( km_sfc  + 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(k,j,i)  + 1.0E-20_wp )
+                                        ( km_sfc  + 1.0E-20_wp )
 
           IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) )  >                        &
@@ -4591,4 +4608,5 @@
 ! ------------
 !> Computation of the turbulent diffusion coefficients for momentum and heat.
+!> @bug unstable stratification is not properly considered for kh in rans mode.
 !------------------------------------------------------------------------------!
  SUBROUTINE tcm_diffusivities( var, var_reference )
@@ -4642,7 +4660,84 @@
 
 !
-!-- Set boundary values (Neumann conditions at downward facing surfaces,
-!-- according to MOST at upward facing surfaces and vertical surfaces).
-!
+!-- In RANS mode, use MOST to calculate km and kh within the surface layer.
+    IF ( rans_tke_e )  THEN
+!
+!--    Upward facing surfaces
+!--    Default surfaces
+       n = 0
+       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
+       DO  m = 1, surf_def_h(0)%ns
+          i = surf_def_h(0)%i(m)
+          j = surf_def_h(0)%j(m)
+          k = surf_def_h(0)%k(m)
+          km(k,j,i) = kappa * surf_def_h(0)%us(m) * surf_def_h(0)%z_mo(m) /    &
+                      phi_m( surf_def_h(0)%z_mo(m) / surf_def_h(0)%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 not consider stratification at these surfaces.
+       DO  n = 0, 3
+!      
+!--       Default surfaces
+          !$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)
+             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_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
+          !$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)
+             kh(k,j,i) = 1.35_wp * km(k,j,i)
+          ENDDO
+       ENDDO
+       
+       CALL exchange_horiz( km, nbgp )
+       CALL exchange_horiz( kh, nbgp )
+
+    ENDIF
+!
+!-- Set boundary values (Neumann conditions)
 !-- Downward facing surfaces
     !$OMP PARALLEL DO PRIVATE(i,j,k)
@@ -4654,81 +4749,14 @@
        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)
+!
+!-- Downward facing surfaces
+    !$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
-!
-!-- 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
-       !$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
-!
-!--    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
-!
-!--    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
-    ENDDO
-
-    CALL exchange_horiz( km, nbgp )
-    CALL exchange_horiz( kh, nbgp )
-
 !
 !-- Model top
