Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 3119)
+++ /palm/trunk/SOURCE/modules.f90	(revision 3120)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! +les_dynamic
+! 
+! 3083 2018-06-19 14:03:12Z gronemeier
 ! set dt_3d = 0.01
 ! 
@@ -1306,4 +1309,5 @@
     LOGICAL ::  large_scale_subsidence = .FALSE.                 !< namelist parameter
     LOGICAL ::  land_surface = .FALSE.                           !< use land surface model?
+    LOGICAL ::  les_dynamic = .FALSE.                            !< use dynamic subgrid model as turbulence closure for LES mode
     LOGICAL ::  les_mw = .FALSE.                                 !< use Moeng-Wyngaard turbulence closure for LES mode
     LOGICAL ::  lsf_exception = .FALSE.                          !< use of lsf with buildings (temporary)?
Index: /palm/trunk/SOURCE/timestep.f90
===================================================================
--- /palm/trunk/SOURCE/timestep.f90	(revision 3119)
+++ /palm/trunk/SOURCE/timestep.f90	(revision 3120)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Put ABS( km ) in computation of time step according to the diffusion criterion
+! 
+! 3084 2018-06-19 15:30:55Z gronemeier
 ! limit increase of dt_3d only in case of RANS mode
 ! 
@@ -302,5 +305,6 @@
 !
 !--    Compute time step according to the diffusion criterion.
-!--    First calculate minimum grid spacing which only depends on index k
+!--    First calculate minimum grid spacing which only depends on index k.
+!--    When using the dynamic subgrid model, negative km are possible.
        dt_diff_l = 999999.0_wp
 
@@ -315,5 +319,6 @@
              DO  k = nzb+1, nzt
                 dt_diff_l = MIN( dt_diff_l, dxyz2_min(k) /                     &
-                                    ( MAX( kh(k,j,i), km(k,j,i) ) + 1E-20_wp ) )
+                                    ( MAX( kh(k,j,i), ABS( km(k,j,i) ) )       &
+                                      + 1E-20_wp ) )
              ENDDO
           ENDDO
Index: /palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- /palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3119)
+++ /palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3120)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! - changed tcm_diffusivities to tcm_diffusivities_default
+! - created subroutine tcm_diffusivities that calls tcm_diffusivities_default
+!   and tcm_diffusivities_dynamic
+! 
+! 3086 2018-06-25 09:08:04Z gronemeier
 ! bugfix: set rans_const_sigma(1) = 1.3
 ! 
@@ -95,5 +100,5 @@
 ! --------
 ! @author Tobias Gronemeier
-!
+! @author Hauke Wurps
 !
 ! Description:
@@ -126,9 +131,10 @@
         ONLY:  constant_diffusion, dt_3d, e_init, humidity, inflow_l,          &
                initializing_actions, intermediate_timestep_count,              &
-               intermediate_timestep_count_max, kappa, km_constant, les_mw,    &
-               ocean, plant_canopy, prandtl_number, prho_reference,            &
-               pt_reference, rans_mode, rans_tke_e, rans_tke_l, simulated_time,&
-               timestep_scheme, turbulence_closure, turbulent_inflow,          &
-               use_upstream_for_tke, vpt_reference, ws_scheme_sca
+               intermediate_timestep_count_max, kappa, km_constant,            &
+               les_dynamic, les_mw, ocean, plant_canopy, prandtl_number,       &
+               prho_reference, pt_reference, rans_mode, rans_tke_e, rans_tke_l,&
+               simulated_time,timestep_scheme, turbulence_closure,             &
+               turbulent_inflow, use_upstream_for_tke, vpt_reference,          &
+               ws_scheme_sca
 
     USE advec_ws,                                                              &
@@ -148,6 +154,5 @@
 
     USE indices,                                                               &
-        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg,               &
-               nzb, nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt,   &
+        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,     &
                wall_flags_0
 
@@ -296,8 +301,27 @@
 
 !
-!-- Calculate diffusivities
+!-- Call tcm_diffusivities_default and tcm_diffusivities_dynamic
     INTERFACE tcm_diffusivities
        MODULE PROCEDURE tcm_diffusivities
     END INTERFACE tcm_diffusivities
+
+!
+!-- Calculate diffusivities
+    INTERFACE tcm_diffusivities_default
+       MODULE PROCEDURE tcm_diffusivities_default
+    END INTERFACE tcm_diffusivities_default
+
+!
+!-- Calculate diffusivities according to dynamic sgs model
+    INTERFACE tcm_diffusivities_dynamic
+       MODULE PROCEDURE tcm_diffusivities_dynamic
+    END INTERFACE tcm_diffusivities_dynamic
+
+!
+!-- Box-filter method for dynamic sgs model
+    INTERFACE tcm_box_filter_2d
+       MODULE PROCEDURE tcm_box_filter_2d_single
+       MODULE PROCEDURE tcm_box_filter_2d_array
+    END INTERFACE tcm_box_filter_2d
 
 !
@@ -390,4 +414,7 @@
           CASE ( 'Moeng_Wyngaard' )
              les_mw = .TRUE.
+
+          CASE ( 'dynamic' )
+             les_dynamic = .TRUE.
 
           CASE DEFAULT
@@ -1945,8 +1972,8 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  i   !< grid index x-direction
-    INTEGER(iwp) ::  j   !< grid index y-direction
-    INTEGER(iwp) ::  k   !< grid index z-direction
-    INTEGER(iwp) ::  m   !< running index surface elements
+    INTEGER(iwp) ::  i      !< grid index x-direction
+    INTEGER(iwp) ::  j      !< grid index y-direction
+    INTEGER(iwp) ::  k      !< grid index z-direction
+    INTEGER(iwp) ::  m      !< running index surface elements
 
     IF ( constant_flux_layer )  THEN
@@ -1973,13 +2000,13 @@
 !--       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.
           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 )  
+                                     drho_air_zw(k-1)               *          &
+                                     ( zu(k+1) - zu(k-1)    )       /          &
+                                     ( 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(k,j,i)  + 1.0E-20_wp )  
+                                     drho_air_zw(k-1)               *          &
+                                     ( zu(k+1) - zu(k-1)    )       /          &
+                                     ( km(k,j,i)  + 1.0E-20_wp )   
 
           IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) )  >                     &
@@ -2032,12 +2059,13 @@
 !--       interpolation of km onto the u/v-grid is necessary. However, the
 !--       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(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(k,j,i) + 1.0E-20_wp )
+!--       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(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(k,j,i)  + 1.0E-20_wp )
 
           IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) )  >                        &
@@ -2063,12 +2091,12 @@
 !--       interpolation of km onto the u/v-grid is necessary. However, the
 !--       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(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(k,j,i) + 1.0E-20_wp )
+          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 )
+          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 )
 
           IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) )  >                        &
@@ -4230,5 +4258,5 @@
 !
 !--          Calculate dissipation
-             IF ( les_mw )  THEN
+             IF ( les_dynamic .OR. les_mw )  THEN
 
                 CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
@@ -4367,5 +4395,5 @@
 !--    Calculate dissipation...
 !--    ...in case of LES
-       IF ( les_mw )  THEN
+       IF ( les_dynamic .OR. les_mw )  THEN
 
           CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
@@ -4678,9 +4706,35 @@
 ! Description:
 ! ------------
+!> Computation of the turbulent diffusion coefficients for momentum and heat.
+!------------------------------------------------------------------------------!
+ SUBROUTINE tcm_diffusivities( var, var_reference )
+
+    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
+!
+!-- 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.
 !> @todo consider non-default surfaces
 !------------------------------------------------------------------------------!
- SUBROUTINE tcm_diffusivities( var, var_reference )
+ SUBROUTINE tcm_diffusivities_default( var, var_reference )
  
 
@@ -4746,5 +4800,5 @@
     ENDIF
 
-    IF ( les_mw )  THEN
+    IF ( les_dynamic .OR. les_mw )  THEN
        !$OMP DO
        DO  i = nxlg, nxrg
@@ -4926,5 +4980,369 @@
     ENDIF
 
- END SUBROUTINE tcm_diffusivities
+ END SUBROUTINE tcm_diffusivities_default
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculates the eddy viscosity dynamically using the linear dynamic model 
+!> according to 
+!> Heinz, Stefan. "Realizability of dynamic subgrid-scale stress models via 
+!> stochastic analysis." 
+!> Monte Carlo Methods and Applications 14.4 (2008): 311-329.
+!>
+!> Furthermore dynamic bounds are used to limit the absolute value of c* as 
+!> described in
+!> Mokhtarpoor, Reza, and Stefan Heinz. "Dynamic large eddy simulation: 
+!> Stability via realizability." Physics of Fluids 29.10 (2017): 105104.
+!>
+!> @author Hauke Wurps
+!------------------------------------------------------------------------------!
+ SUBROUTINE tcm_diffusivities_dynamic
+        
+    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
+
+    INTEGER(iwp) ::  i           !< running index x-direction
+    INTEGER(iwp) ::  j           !< running index y-direction
+    INTEGER(iwp) ::  k           !< running index z-direction
+    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
+    REAL(wp)     ::  dudy        !< Gradient of u-component in y-direction
+    REAL(wp)     ::  dudz        !< Gradient of u-component in z-direction
+    REAL(wp)     ::  dvdx        !< Gradient of v-component in x-direction
+    REAL(wp)     ::  dvdy        !< Gradient of v-component in y-direction
+    REAL(wp)     ::  dvdz        !< Gradient of v-component in z-direction
+    REAL(wp)     ::  dwdx        !< Gradient of w-component in x-direction
+    REAL(wp)     ::  dwdy        !< Gradient of w-component in y-direction
+    REAL(wp)     ::  dwdz        !< Gradient of w-component in z-direction
+
+    REAL(wp)     ::  uc(-1:1,-1:1)  !< u on grid center
+    REAL(wp)     ::  vc(-1:1,-1:1)  !< v on grid center
+    REAL(wp)     ::  wc(-1:1,-1:1)  !< w on grid center
+    REAL(wp)     ::  u2(-1:1,-1:1)  !< u2 on grid center
+    REAL(wp)     ::  v2(-1:1,-1:1)  !< v2 on grid center
+    REAL(wp)     ::  w2(-1:1,-1:1)  !< w2 on grid center
+    REAL(wp)     ::  uv(-1:1,-1:1)  !< u*v on grid center
+    REAL(wp)     ::  uw(-1:1,-1:1)  !< u*w on grid center
+    REAL(wp)     ::  vw(-1:1,-1:1)  !< v*w on grid center
+
+    REAL(wp)     ::  ut(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !< test filtered u
+    REAL(wp)     ::  vt(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !< test filtered v
+    REAL(wp)     ::  wt(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !< test filtered w
+
+    REAL(wp)     ::  uct         !< test filtered u on grid center
+    REAL(wp)     ::  vct         !< test filtered v on grid center
+    REAL(wp)     ::  wct         !< test filtered w on grid center
+    REAL(wp)     ::  u2t         !< test filtered u**2 on grid center
+    REAL(wp)     ::  v2t         !< test filtered v**2 on grid center
+    REAL(wp)     ::  w2t         !< test filtered w**2 on grid center
+    REAL(wp)     ::  uvt         !< test filtered u*v on grid center
+    REAL(wp)     ::  uwt         !< test filtered u*w on grid center
+    REAL(wp)     ::  vwt         !< test filtered v*w on grid center
+
+    REAL(wp)     ::  sd11        !< deviatoric shear tensor
+    REAL(wp)     ::  sd22        !< deviatoric shear tensor
+    REAL(wp)     ::  sd33        !< deviatoric shear tensor
+    REAL(wp)     ::  sd12        !< deviatoric shear tensor
+    REAL(wp)     ::  sd13        !< deviatoric shear tensor
+    REAL(wp)     ::  sd23        !< deviatoric shear tensor
+
+    REAL(wp)     ::  sd2         !< sum: sd_ij*sd_ij
+
+    REAL(wp)     ::  sdt11       !< filtered deviatoric shear tensor
+    REAL(wp)     ::  sdt22       !< filtered deviatoric shear tensor
+    REAL(wp)     ::  sdt33       !< filtered deviatoric shear tensor
+    REAL(wp)     ::  sdt12       !< filtered deviatoric shear tensor
+    REAL(wp)     ::  sdt13       !< filtered deviatoric shear tensor
+    REAL(wp)     ::  sdt23       !< filtered deviatoric shear tensor
+
+    REAL(wp)     ::  sdt2        !< sum: sdt_ij*sdt_ij
+
+    REAL(wp)     ::  ld11        !< deviatoric stress tensor
+    REAL(wp)     ::  ld22        !< deviatoric stress tensor
+    REAL(wp)     ::  ld33        !< deviatoric stress tensor
+    REAL(wp)     ::  ld12        !< deviatoric stress tensor
+    REAL(wp)     ::  ld13        !< deviatoric stress tensor
+    REAL(wp)     ::  ld23        !< deviatoric stress tensor
+
+    REAL(wp)     ::  lnn         !< sum ld_nn
+    REAL(wp)     ::  ldsd        !< sum: ld_ij*sd_ij
+
+    REAL(wp)     ::  delta       !< grid size
+    REAL(wp)     ::  cst         !< c*
+    REAL(wp)     ::  cstnust_t   !< product c*nu*
+    REAL(wp)     ::  cst_max     !< bounds of c*
+
+    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:
+    CALL tcm_box_filter_2d( u, ut )
+    CALL tcm_box_filter_2d( v, vt )
+    CALL tcm_box_filter_2d( w, wt )
+
+    DO  i = nxl, nxr
+       DO  j = nys, nyn
+          DO  k = nzb+1, nzt
+!
+!--          Compute the deviatoric shear tensor s_ij on grid centers:
+!--          s_ij =  0.5 * ( du_i/dx_j + du_j/dx_i )
+             dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+             dudy  = 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  = 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  = 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  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+             dvdz  = 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  = 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  = 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  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+             sd11 = dudx
+             sd22 = dvdy
+             sd33 = dwdz
+             sd12 = 0.5_wp * ( dudy + dvdx )
+             sd13 = 0.5_wp * ( dudz + dwdx )
+             sd23 = 0.5_wp * ( dvdz + dwdy )
+!
+!--          sum: sd_ij*sd_ij
+             sd2 = sd11**2 + sd22**2 + sd33**2                     &
+                   + 2.0_wp * ( sd12**2 + sd13**2 + sd23**2 )
+!
+!--          The filtered velocities are needed to calculate the filtered shear
+!--          tensor: sdt_ij = 0.5 * ( dut_i/dx_j + dut_j/dx_i )
+             dudx  =           ( ut(k,j,i+1) - ut(k,j,i)     ) * ddx
+             dudy  = 0.25_wp * ( ut(k,j+1,i) + ut(k,j+1,i+1) - &
+                                 ut(k,j-1,i) - ut(k,j-1,i+1) ) * ddy
+             dudz  = 0.5_wp  * ( ut(k+1,j,i) + ut(k+1,j,i+1) - &
+                                 ut(k-1,j,i) - ut(k-1,j,i+1) ) * dd2zu(k)
+
+             dvdx  = 0.25_wp * ( vt(k,j,i+1) + vt(k,j+1,i+1) - &
+                                 vt(k,j,i-1) - vt(k,j+1,i-1) ) * ddx
+             dvdy  =           ( vt(k,j+1,i) - vt(k,j,i)     ) * ddy
+             dvdz  = 0.5_wp  * ( vt(k+1,j,i) + vt(k+1,j+1,i) - &
+                                 vt(k-1,j,i) - vt(k-1,j+1,i) ) * dd2zu(k)
+
+             dwdx  = 0.25_wp * ( wt(k,j,i+1) + wt(k-1,j,i+1) - &
+                                 wt(k,j,i-1) - wt(k-1,j,i-1) ) * ddx
+             dwdy  = 0.25_wp * ( wt(k,j+1,i) + wt(k-1,j+1,i) - &
+                                 wt(k,j-1,i) - wt(k-1,j-1,i) ) * ddy
+             dwdz  =           ( wt(k,j,i)   - wt(k-1,j,i)   ) * ddzw(k)
+
+             sdt11 = dudx
+             sdt22 = dvdy
+             sdt33 = dwdz
+             sdt12 = 0.5_wp * ( dudy + dvdx )
+             sdt13 = 0.5_wp * ( dudz + dwdx )
+             sdt23 = 0.5_wp * ( dvdz + dwdy )
+!
+!--          sum: sd_ij*sd_ij
+             sdt2 = sdt11**2 + sdt22**2 + sdt33**2         &
+                   + 2.0_wp * ( sdt12**2 + sdt13**2 + sdt23**2 )
+!
+!--          Need filtered velocities and filtered squared velocities on grid
+!--          centers. Substraction of geostrophic velocity helps to avoid
+!--          numerical errors in the expression <u**2> - <u>*<u>, which can be
+!--          very small (<...> means filtered).
+             DO  l = -1, 1
+                DO  m = -1, 1
+                   uc(l,m) = 0.5_wp * ( u(k,j+l,i+m) + u(k,j+l,i+m+1) ) - ug(k)
+                   vc(l,m) = 0.5_wp * ( v(k,j+l,i+m) + v(k,j+l+1,i+m) ) - vg(k)
+                   wc(l,m) = 0.5_wp * ( w(k-1,j+l,i+m) + w(k,j+l,i+m) )
+                ENDDO
+             ENDDO
+
+             CALL tcm_box_filter_2d( uc, uct )
+             CALL tcm_box_filter_2d( vc, vct )
+             CALL tcm_box_filter_2d( wc, wct )
+             CALL tcm_box_filter_2d( uc**2, u2t )
+             CALL tcm_box_filter_2d( vc**2, v2t )
+             CALL tcm_box_filter_2d( wc**2, w2t )
+             CALL tcm_box_filter_2d( uc*vc, uvt )
+             CALL tcm_box_filter_2d( uc*wc, uwt )
+             CALL tcm_box_filter_2d( vc*wc, vwt )
+
+             ld11 = u2t - uct*uct
+             ld22 = v2t - vct*vct
+             ld33 = w2t - wct*wct
+             ld12 = uvt - uct*vct
+             ld13 = uwt - uct*wct
+             ld23 = vwt - vct*wct
+
+             lnn = ld11 + ld22 + ld33
+!
+!--          Substract trace to get deviatoric resolved stress
+             ld11 = ld11 - lnn / 3.0_wp
+             ld22 = ld22 - lnn / 3.0_wp
+             ld33 = ld33 - lnn / 3.0_wp
+
+             ldsd = ld11 * sdt11 + ld22 * sdt22 + ld33 * sdt33 + &
+                    2.0_wp * ( ld12 * sdt12 + ld13 * sdt13 + ld23 * sdt23 )
+!
+!--          c* nu*^T is SGS viscosity on test filter level:
+             cstnust_t = -ldsd / sdt2
+!
+!--          The model was only tested for an isotropic grid. The following
+!--          expression was a recommendation of Stefan Heinz.
+             delta = MAX( dx, dy, dzw(k) )
+             cst = cstnust_t / ( 4.0_wp * delta * SQRT( lnn / 2.0_wp ) )
+!
+!--          Calculate border according to Mokhtarpoor and Heinz (2017)
+             cst_max = fac_cmax * SQRT( e(k,j,i) ) / ( delta * SQRT( 2.0_wp * sd2 ) )
+
+             IF ( ABS( cst ) > cst_max )  THEN
+                cst = cst_max * cst / ABS( cst )
+             ENDIF
+
+             km(k,j,i) = cst * delta * SQRT( e(k,j,i) )
+
+          ENDDO
+       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
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This subroutine acts as a box filter with filter width 2 * dx.
+!> Output is only one point.
+!------------------------------------------------------------------------------!
+ SUBROUTINE tcm_box_filter_2d_single( var, var_fil )
+ 
+    IMPLICIT NONE
+
+    REAL(wp)     ::  var(-1:1,-1:1)      !< variable to be filtered
+    REAL(wp)     ::  var_fil             !< filtered variable
+!
+!-- It is assumed that a box with a side length of 2 * dx and centered at the 
+!-- variable*s position contains one half of the four closest neigbours and one
+!-- forth of the diagonally closest neighbours.
+    var_fil = 0.25_wp * ( var(0,0) +                                           &
+                      0.5_wp  * ( var(0,1)  + var(1,0)   +                     &
+                                  var(0,-1) + var(-1,0)  ) +                   &
+                      0.25_wp * ( var(1,1)  + var(1,-1)  +                     &
+                                  var(-1,1) + var(-1,-1) ) )
+
+ END SUBROUTINE tcm_box_filter_2d_single
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This subroutine acts as a box filter with filter width 2 * dx.
+!> The filtered variable var_fil is on the same grid as var.
+!------------------------------------------------------------------------------!
+ SUBROUTINE tcm_box_filter_2d_array( var, var_fil )
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  i    !< running index x-direction
+    INTEGER(iwp) ::  j    !< running index y-direction
+    INTEGER(iwp) ::  k    !< running index z-direction
+
+    REAL(wp)     ::  var(nzb:nzt+1,nysg:nyng,nxlg:nxrg)      !< variable to be filtered
+    REAL(wp)     ::  var_fil(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !< filtered variable
+!
+!-- It is assumed that a box with a side length of 2 * dx and centered at the 
+!-- variable's position contains one half of the four closest neigbours and one
+!-- forth of the diagonally closest neighbours.
+    DO  i = nxlg+1, nxrg-1
+       DO  j = nysg+1, nyng-1
+          DO  k = nzb, nzt+1
+             var_fil(k,j,i) = 0.25_wp * ( var(k,j,i) +                         &
+                              0.5_wp  * ( var(k,j,i+1)   + var(k,j+1,i)   +    &
+                                          var(k,j,i-1)   + var(k,j-1,i)     ) +&
+                              0.25_wp * ( var(k,j+1,i+1) + var(k,j+1,i-1) +    &
+                                          var(k,j-1,i+1) + var(k,j-1,i-1)   )  )
+          END DO
+       END DO
+    END DO
+
+ END SUBROUTINE tcm_box_filter_2d_array
 
 
