Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 4474)
+++ palm/trunk/SOURCE/modules.f90	(revision 4476)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Renamed variables for subgrids-scale model steering
+! 
+! 4473 2020-03-25 21:04:07Z gronemeier
 ! moved wall_adjustment_factor to turbulence_closure_mod
 !
@@ -548,5 +551,5 @@
     CHARACTER (LEN=20)   ::  reference_state = 'initial_profile'          !< namelist parameter
     CHARACTER (LEN=20)   ::  timestep_scheme = 'runge-kutta-3'            !< namelist parameter
-    CHARACTER (LEN=20)   ::  turbulence_closure = 'Moeng_Wyngaard'        !< namelist parameter
+    CHARACTER (LEN=20)   ::  turbulence_closure = '1.5-order'             !< namelist parameter
     CHARACTER (LEN=40)   ::  topography = 'flat'                          !< namelist parameter
     CHARACTER (LEN=64)   ::  host = '????'                                !< configuration identifier as given by palmrun option -c, ENVPAR namelist parameter provided by palmrun
@@ -750,6 +753,7 @@
     LOGICAL ::  large_scale_subsidence = .FALSE.                 !< namelist parameter
     LOGICAL ::  land_surface = .FALSE.                           !< use land surface model?
+    LOGICAL ::  les_dai = .FALSE.                                !< use Dai et al. turbulence closure (modified 1.5-order closure) for LES mode. Shall replace the default 1.5-order closure
     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 ::  les_default = .FALSE.                            !< use 1.5-order default turbulence closure for LES mode
     LOGICAL ::  lsf_exception = .FALSE.                          !< use of lsf with buildings (temporary)?
     LOGICAL ::  lsf_surf = .TRUE.                                !< use surface forcing (large scale forcing)?
Index: palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 4474)
+++ palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 4476)
@@ -25,4 +25,18 @@
 ! -----------------
 ! $Id$
+! - added new LES closure after Dai et al. (2020), which provides much better grid
+!   convergence in stable boundary layer runs. The implementation is experimental
+!   at the moment and should be used with special care. The new SGS closure can be
+!   switched on via turbulence_closure = '1.5-order-dai'
+! - variable ml_wall_adjusted renamed to delta as it represents a grid size and
+!   not a mixing length (see Equations 14 and 18 in Maronga et al. 2015, GMD)
+! - nameing of turbulence closures revised:
+!        'Moeng_Wyngaard' to '1.5-order'
+!        'TKE-l' to 'tke-l'
+!        'TKE-e' to 'tke-e'
+! - LOGICAL steering variable renamed:
+!        les_mw to les_default
+! 
+! 4473 2020-03-25 21:04:07Z gronemeier
 ! - rename l-grid to gridsize-geometric-mean
 !          l-wall to ml-wall-adjusted
@@ -140,5 +154,6 @@
                initializing_actions, intermediate_timestep_count,             &
                intermediate_timestep_count_max, km_constant,                  &
-               les_dynamic, les_mw, ocean_mode, plant_canopy, prandtl_number, &
+               les_dai, les_dynamic, les_default,                             &
+               ocean_mode, plant_canopy, prandtl_number,                      &
                pt_reference, rans_mode, rans_tke_e, rans_tke_l,               &
                timestep_scheme, turbulence_closure,                           &
@@ -210,5 +225,6 @@
 
     REAL(wp), ALLOCATABLE, DIMENSION(:)     ::  ml_blackadar             !< mixing length according to Blackadar
-    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ml_wall_adjusted         !< mixing length adjusted according to near-by walls
+    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  delta                    !< grid size, possibly limited by wall adjustment factor
+    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  distance_to_wall         !< distance to the surface/wall
 
 !
@@ -584,15 +600,18 @@
     SELECT CASE ( TRIM( turbulence_closure ) )
 
+       CASE ( '1.5-order' )
+          les_default = .TRUE.
+
+       CASE ( '1.5-order-dai' )
+          les_dai = .TRUE.
+
        CASE ( 'dynamic' )
           les_dynamic = .TRUE.
 
-       CASE ( 'Moeng_Wyngaard' )
-          les_mw = .TRUE.
-
-       CASE ( 'TKE-l' )
+       CASE ( 'tke-l' )
           rans_tke_l = .TRUE.
           rans_mode = .TRUE.
 
-       CASE ( 'TKE-e' )
+       CASE ( 'tke-e' )
           rans_tke_e = .TRUE.
           rans_mode = .TRUE.
@@ -1157,5 +1176,5 @@
                    DO  j = nysg, nyng
                       DO  k = nzb+1, nzt
-                         km(k,j,i) = c_0 * SQRT( e_init ) * MIN( ml_wall_adjusted(k,j,i), &
+                         km(k,j,i) = c_0 * SQRT( e_init ) * MIN( delta(k,j,i), &
                                                                  ml_blackadar(k) )
                       ENDDO
@@ -1202,5 +1221,5 @@
                 DO  j = nysg, nyng
                    DO  k = nzb+1, nzt
-                      km(k,j,i) = c_0 * SQRT( e_init ) * ml_wall_adjusted(k,j,i)
+                      km(k,j,i) = c_0 * SQRT( e_init ) * delta(k,j,i)
                    ENDDO
                 ENDDO
@@ -1423,5 +1442,5 @@
     REAL(wp) :: distance_edge_yz_down  !< distance of grid box center to its boundary along yz diagonal (down)
     REAL(wp) :: distance_edge_yz_up    !< distance of grid box center to its boundary along yz diagonal (up)
-    REAL(wp) :: distance_edge_xz_down  !< distance of grid box center to its boundary along xz diagonal (down)
+    REAL(wp) :: distance_edge_xz_down  !< distance of grid box center to its boundary along xz diagonal
     REAL(wp) :: distance_edge_xz_up    !< distance of grid box center to its boundary along xz diagonal (up)
     REAL(wp) :: distance_edge_xy       !< distance of grid box center to its boundary along xy diagonal
@@ -1432,6 +1451,11 @@
     REAL(wp), DIMENSION(nzb:nzt+1) ::  gridsize_geometric_mean  !< geometric mean of grid sizes dx, dy, dz
 
-    ! ALLOCATE( gridsize_geometric_mean(nzb:nzt+1) )
-    ALLOCATE( ml_wall_adjusted(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+
+    ALLOCATE( delta(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+
+
+    IF ( les_dai )  THEN
+       ALLOCATE ( distance_to_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+    ENDIF
 !
 !-- Initialize the mixing length in case of an LES-simulation
@@ -1454,7 +1478,13 @@
 
 !
-!--    Initialize near-wall mixing length ml_wall_adjusted
+!--    Initialize the grid size variable (delta) and the distance to the wall, if needed
        DO  k = nzb, nzt+1
-          ml_wall_adjusted(k,:,:) = gridsize_geometric_mean(k)
+          delta(k,:,:) = gridsize_geometric_mean(k)
+
+!
+!--       If Dai et al. (2020) closure is used, the distance to the wall must be stored
+          IF ( les_dai )  THEN
+             distance_to_wall(k,:,:) = zu(k)
+          ENDIF
        ENDDO
 
@@ -1569,7 +1599,9 @@
                          distance_corners_up = 9999999.9_wp
                       ENDIF
-!
-!--                   Adjust mixing length by wall-adjustment factor
-                      ml_wall_adjusted(k,j,i) = wall_adjustment_factor * MIN(  &
+
+!
+!--                   Calculate the minimum distance from the wall and store it
+!--                   temporarily in the array delta
+                      delta(k,j,i) = MIN(                                      &
                          distance_up, distance_down, distance_ns, distance_lr, &
                          distance_edge_yz_down, distance_edge_yz_up,           &
@@ -1578,8 +1610,24 @@
                          distance_corners_down, distance_corners_up )
 
+!
+!--                   If Dai et al. (2020) closure is used, the distance to the wall
+!--                   must be permanently stored
+                      IF ( les_dai  .AND.  delta(k,j,i) /= 9999999.9_wp )  THEN
+                         distance_to_wall(k,j,i) = delta(k,j,i)
+                      ENDIF
+
+!
+!--                   Close to the surface, the maximum mixing length is limited to 1.8 z
+                      delta(k,j,i) = wall_adjustment_factor * delta(k,j,i)
+
                   ENDIF  !if grid point belongs to atmosphere
 
-                  ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i), &
-                                                 gridsize_geometric_mean(k) )
+
+
+!
+!--               The grid size (delta) is defined as the the minimum of the distance to
+!--               the nearest wall * 1.8 and the geometric mean grid size.
+                  delta(k,j,i) = MIN( delta(k,j,i), gridsize_geometric_mean(k) )
+
 
                 ENDDO  !k loop
@@ -1622,6 +1670,6 @@
        ENDDO
 
-       ml_wall_adjusted(nzb,:,:) = ml_blackadar(nzb)
-       ml_wall_adjusted(nzt+1,:,:) = ml_blackadar(nzt+1)
+       delta(nzb,:,:) = ml_blackadar(nzb)
+       delta(nzt+1,:,:) = ml_blackadar(nzt+1)
 !
 !--    Limit mixing length to either nearest wall or Blackadar mixing length.
@@ -1634,6 +1682,6 @@
           radius = ml_blackadar(k)  ! radius within walls are searched
 !
-!--       Set ml_wall_adjusted to its default maximum value (ml_blackadar)
-          ml_wall_adjusted(k,:,:) = radius
+!--       Set delta to its default maximum value (ml_blackadar)
+          delta(k,:,:) = radius
 
 !
@@ -1730,5 +1778,5 @@
 !--                            Search for walls within octant (+++)
                                vic_yz = vicinity(0:rad_k+1,0:rad_j,ii)
-                               ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),             &
+                               delta(k,j,i) = MIN( delta(k,j,i),             &
                                        shortest_distance( vic_yz, .TRUE., ii ) )
 !
@@ -1738,5 +1786,5 @@
 !--                            shortest_distance").
                                vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii)
-                               ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),             &
+                               delta(k,j,i) = MIN( delta(k,j,i),             &
                                        shortest_distance( vic_yz, .TRUE., ii ) )
 
@@ -1745,14 +1793,14 @@
 !--                         Search for walls within octant (+--)
                             vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii)
-                            ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),                &
+                            delta(k,j,i) = MIN( delta(k,j,i),                &
                                       shortest_distance( vic_yz, .FALSE., ii ) )
 !
 !--                         Search for walls within octant (++-)
                             vic_yz = vicinity(0:-rad_k-1:-1,0:rad_j,ii)
-                            ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),                &
+                            delta(k,j,i) = MIN( delta(k,j,i),                &
                                       shortest_distance( vic_yz, .FALSE., ii ) )
 !
 !--                         Reduce search along x by already found distance
-                            dist_dx = CEILING( ml_wall_adjusted(k,j,i) / dx )
+                            dist_dx = CEILING( delta(k,j,i) / dx )
 
                          ENDDO
@@ -1767,5 +1815,5 @@
 !--                            Search for walls within octant (-++)
                                vic_yz = vicinity(0:rad_k+1,0:rad_j,ii)
-                               ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),             &
+                               delta(k,j,i) = MIN( delta(k,j,i),             &
                                       shortest_distance( vic_yz, .TRUE., -ii ) )
 !
@@ -1775,5 +1823,5 @@
 !--                            shortest_distance").
                                vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii)
-                               ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),             &
+                               delta(k,j,i) = MIN( delta(k,j,i),             &
                                       shortest_distance( vic_yz, .TRUE., -ii ) )
 
@@ -1782,14 +1830,14 @@
 !--                         Search for walls within octant (---)
                             vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii)
-                            ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),                &
+                            delta(k,j,i) = MIN( delta(k,j,i),                &
                                      shortest_distance( vic_yz, .FALSE., -ii ) )
 !
 !--                         Search for walls within octant (-+-)
                             vic_yz = vicinity(0:-rad_k-1:-1,0:rad_j,ii)
-                            ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),                &
+                            delta(k,j,i) = MIN( delta(k,j,i),                &
                                      shortest_distance( vic_yz, .FALSE., -ii ) )
 !
 !--                         Reduce search along x by already found distance
-                            dist_dx = CEILING( ml_wall_adjusted(k,j,i) / dx )
+                            dist_dx = CEILING( delta(k,j,i) / dx )
 
                          ENDDO
@@ -1799,5 +1847,5 @@
                    ELSE  !Check if (i,j,k) belongs to atmosphere
 
-                      ml_wall_adjusted(k,j,i) = ml_blackadar(k)
+                      delta(k,j,i) = ml_blackadar(k)
 
                    ENDIF
@@ -1818,8 +1866,14 @@
 
 !
-!-- Set lateral boundary conditions for ml_wall_adjusted
-    CALL exchange_horiz( ml_wall_adjusted, nbgp )
-
-    !$ACC ENTER DATA COPYIN(ml_wall_adjusted(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
+!-- Set lateral boundary conditions for delta
+    CALL exchange_horiz( delta, nbgp )
+    !$ACC ENTER DATA COPYIN(delta(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
+
+    IF ( les_dai )  THEN
+       CALL exchange_horiz( distance_to_wall, nbgp )
+       !$ACC ENTER DATA COPYIN(distance_to_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
+    ENDIF
+
+
 
     CONTAINS
@@ -4032,10 +4086,10 @@
 !
 !-- Calculate the dissipation
-    IF ( les_dynamic .OR. les_mw )  THEN
+    IF( les_default  .OR.  les_dynamic )  THEN
 
        !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
        !$ACC PRIVATE(ml, ml_stratification, dvar_dz) &
        !$ACC PRESENT(diss, e, var, wall_flags_total_0) &
-       !$ACC PRESENT(dd2zu, ml_wall_adjusted)
+       !$ACC PRESENT(dd2zu, delta)
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -4053,16 +4107,19 @@
                    ENDIF
                 ELSE
-                   ml_stratification(k) = ml_wall_adjusted(k,j,i)
+                   ml_stratification(k) = delta(k,j,i)
                 ENDIF
 
              ENDDO
 
+!
+!--          ATTENTION: Don't merge the following loop with the previous one, because this would prohibit proper vectorization by
+!--          the Intel18 compiler. This behaviour might change for future compiler versions.
              !$ACC LOOP PRIVATE(k)
              !DIR$ IVDEP
              DO  k = nzb+1, nzt
 
-                ml = MIN( ml_wall_adjusted(k,j,i), ml_stratification(k) )
-
-                diss(k,j,i) = ( 0.19_wp + 0.74_wp * ml / ml_wall_adjusted(k,j,i) )           &
+                ml = MIN( delta(k,j,i), ml_stratification(k) )
+
+                diss(k,j,i) = ( 0.19_wp + 0.74_wp * ml / delta(k,j,i) )           &
                               * e(k,j,i) * SQRT( e(k,j,i) ) / ml                             &
                               * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
@@ -4073,4 +4130,51 @@
        ENDDO
 
+    ELSEIF ( les_dai )  THEN
+
+       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
+       !$ACC PRIVATE(ml, ml_stratification, dvar_dz) &
+       !$ACC PRESENT(diss, e, var, wall_flags_total_0) &
+       !$ACC PRESENT(dd2zu, delta)
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             !$ACC LOOP PRIVATE(k)
+             DO  k = nzb+1, nzt
+                dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+                IF ( dvar_dz > 0.0_wp ) THEN
+!
+!--                The mixing length is calculated as 1/l = 1/(kappa*z) + 1/Lb, where Lb is
+!--                the stratification term. 1E-5 is added as l is zero at the beginning of
+!--                the simulation.
+                   IF ( use_single_reference_value )  THEN
+                      ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) ) &
+                                             / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
+                   ELSE
+                      ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) ) &
+                                             / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
+                   ENDIF
+                   ml_stratification(k) = 1.0_wp / ( 1.0_wp / ( kappa * distance_to_wall(k,j,i) ) &
+                                                   + 1.0_wp / ml_stratification(k) )
+                ELSE
+                   ml_stratification(k) = delta(k,j,i)
+                ENDIF
+
+             ENDDO
+
+!
+!--          ATTENTION: Don't merge the following loop with the previous one, because this would prohibit proper vectorization by
+!--          the Intel18 compiler. This behaviour might change for future compiler versions.
+             !$ACC LOOP PRIVATE(k)
+             !DIR$ IVDEP
+             DO  k = nzb+1, nzt
+
+                ml          = ml_stratification(k)
+                diss(k,j,i) = ( 0.19_wp + 0.74_wp * ml / delta(k,j,i) )                      &
+                              * e(k,j,i) * SQRT( e(k,j,i) ) / ml                             &
+                              * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+
+             ENDDO
+
+          ENDDO
+       ENDDO
     ELSEIF ( rans_tke_l )  THEN
 
@@ -4078,5 +4182,5 @@
       !$ACC PRIVATE(ml_stratification, duv2_dz2, rif, dvar_dz) &
       !$ACC PRESENT(diss, e, u, v, var, wall_flags_total_0) &
-      !$ACC PRESENT(dd2zu, ml_blackadar, ml_wall_adjusted)
+      !$ACC PRESENT(dd2zu, ml_blackadar, delta)
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -4084,4 +4188,5 @@
 !--          Calculate Richardson-flux number
              IF ( use_single_reference_value )  THEN
+
                 !$ACC LOOP PRIVATE(k)
                 DO  k = nzb+1, nzt
@@ -4096,4 +4201,7 @@
                 ENDDO
              ELSE
+!
+!--             ATTENTION: Don't merge the following loops with the previous one, because this would prohibit proper vectorization
+!--             by the Intel18 compiler. This behaviour might change for future compiler versions.
                 !$ACC LOOP PRIVATE(k)
                 DO  k = nzb+1, nzt
@@ -4124,5 +4232,5 @@
 
                 diss(k,j,i) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) )                         &
-                            / MIN( ml_stratification(k), ml_wall_adjusted(k,j,i) )         &
+                            / MIN( ml_stratification(k), delta(k,j,i) )                    &
                             * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
 
@@ -4243,5 +4351,5 @@
 !-- Calculate the mixing length and dissipation...
 !-- ...in case of LES
-    IF ( les_dynamic .OR. les_mw )  THEN
+    IF ( les_default  .OR.  les_dynamic )  THEN
 
        DO  k = nzb+1, nzt
@@ -4256,5 +4364,5 @@
              ENDIF
           ELSE
-             ml_stratification(k) = ml_wall_adjusted(k,j,i)
+             ml_stratification(k) = delta(k,j,i)
           ENDIF
        ENDDO
@@ -4262,10 +4370,47 @@
        !DIR$ IVDEP
        DO  k = nzb+1, nzt
-          ml  = MIN( ml_wall_adjusted(k,j,i), ml_stratification(k) )
-
-          diss(k,j,i) = ( 0.19_wp + 0.74_wp * ml / ml_wall_adjusted(k,j,i) )           &
+          ml  = MIN( delta(k,j,i), ml_stratification(k) )
+
+          diss(k,j,i) = ( 0.19_wp + 0.74_wp * ml / delta(k,j,i) )                      &
                         * e(k,j,i) * SQRT( e(k,j,i) ) / ml                             &
                         * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
        ENDDO
+
+    ELSEIF ( les_dai )  THEN
+
+       DO  k = nzb+1, nzt
+          dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+          IF ( dvar_dz > 0.0_wp ) THEN
+!
+!--          The mixing length is calculated as 1/l = 1/(kappa*z) + 1/Lb, where Lb is
+!--          the stratification term. 1E-5 is added as l is zero at the beginning of
+!--          the simulation.
+             IF ( use_single_reference_value )  THEN
+                ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) )                           &
+                                       / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
+             ELSE
+                ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) )                           &
+                                       / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
+             ENDIF
+
+             ml_stratification(k) = 1.0_wp / ( 1.0_wp / ( kappa * distance_to_wall(k,j,i) ) &
+                                             + 1.0_wp / ml_stratification(k) )
+
+          ELSE
+             ml_stratification(k) = delta(k,j,i)
+          ENDIF
+       ENDDO
+
+       !DIR$ IVDEP
+       DO  k = nzb+1, nzt
+          ml  = ml_stratification(k)
+
+          diss(k,j,i) = ( 0.19_wp + 0.74_wp * ml / delta(k,j,i) )                      &
+                        * e(k,j,i) * SQRT( e(k,j,i) ) / ml                             &
+                        * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+       ENDDO
+
+
+
 
 !
@@ -4310,5 +4455,5 @@
        DO  k = nzb+1, nzt
           diss(k,j,i) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) )                         &
-                      / MIN( ml_stratification(k), ml_wall_adjusted(k,j,i) )         &
+                      / MIN( ml_stratification(k), delta(k,j,i) )                    &
                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
        ENDDO
@@ -4730,9 +4875,9 @@
 !$  tn = omp_get_thread_num()
 
-    IF ( les_dynamic .OR. les_mw )  THEN
+    IF ( les_default  .OR.  les_dynamic )  THEN
        !$OMP DO
        !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
        !$ACC PRIVATE(dvar_dz, ml, ml_stratification, ml_local_profile) &
-       !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, ml_wall_adjusted) &
+       !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, delta) &
        !$ACC PRESENT(kh, km, sums_l_l, rmask)
        DO  i = nxlg, nxrg
@@ -4745,5 +4890,5 @@
 !--             due to errors when using OpenACC directives. The execution
 !--             crashes reliably if a subroutine is called at this point (the
-!--             reasong for this behaviour is unknown, however).
+!--             reason for this behaviour is unknown, however).
                 dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
                 IF ( dvar_dz > 0.0_wp ) THEN
@@ -4756,5 +4901,5 @@
                    ENDIF
                 ELSE
-                   ml_stratification(k) = ml_wall_adjusted(k,j,i)
+                   ml_stratification(k) = delta(k,j,i)
                 ENDIF
 
@@ -4765,5 +4910,5 @@
              DO  k = nzb+1, nzt
 
-                ml = MIN( ml_wall_adjusted(k,j,i), ml_stratification(k) )         &
+                ml = MIN( delta(k,j,i), ml_stratification(k) )                    &
                    * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
                 ml_local_profile(k) = ml
@@ -4771,5 +4916,5 @@
 !--             Compute diffusion coefficients for momentum and heat
                 km(k,j,i) = c_0 * ml * SQRT( e(k,j,i) )
-                kh(k,j,i) = ( 1.0_wp + 2.0_wp * ml / ml_wall_adjusted(k,j,i) ) * km(k,j,i)
+                kh(k,j,i) = ( 1.0_wp + 2.0_wp * ml / delta(k,j,i) ) * km(k,j,i)
 
              ENDDO
@@ -4786,4 +4931,67 @@
        ENDDO
 
+    ELSEIF ( les_dai )  THEN
+       !$OMP DO
+       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
+       !$ACC PRIVATE(dvar_dz, ml, ml_stratification, ml_local_profile) &
+       !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, delta) &
+       !$ACC PRESENT(kh, km, sums_l_l, rmask)
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
+             !$ACC LOOP PRIVATE(k)
+             DO  k = nzb+1, nzt
+
+                dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+                IF ( dvar_dz > 0.0_wp ) THEN
+!
+!--                The mixing length is calculated as 1/l = 1/(kappa*z) + 1/Lb, where Lb is
+!--                the stratification term. 1E-5 is added as l is zero at the beginning of
+!--                the simulation.
+                   IF ( use_single_reference_value )  THEN
+                      ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) )                           &
+                                             / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
+                   ELSE
+                      ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) )                           &
+                                             / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
+                   ENDIF
+
+                   ml_stratification(k) = 1.0_wp / ( 1.0_wp / ( kappa * distance_to_wall(k,j,i) ) &
+                                        + 1.0_wp / ml_stratification(k) )
+                ELSE
+                   ml_stratification(k) = delta(k,j,i)
+                ENDIF
+
+             ENDDO
+
+             !$ACC LOOP PRIVATE(k)
+             !DIR$ IVDEP
+             DO  k = nzb+1, nzt
+
+                ml_local_profile(k) = ml_stratification(k) * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+                ml  = ml_local_profile(k)
+!
+!--             Compute diffusion coefficients for momentum and heat
+                km(k,j,i) = c_0 * ml * SQRT( e(k,j,i) )
+
+                dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+                IF ( dvar_dz > 0.0_wp ) THEN
+                   kh(k,j,i) = km(k,j,i)
+                ELSE
+                   kh(k,j,i) = 3.0_wp * km(k,j,i)
+                ENDIF
+
+             ENDDO
+!
+!--          Summation for averaged profile (cf. flow_statistics)
+             !$ACC LOOP PRIVATE(sr, k)
+             DO  sr = 0, statistic_regions
+                DO  k = nzb+1, nzt
+                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + ml_local_profile(k) * rmask(j,i,sr)
+                ENDDO
+             ENDDO
+
+          ENDDO
+       ENDDO
+
     ELSEIF ( rans_tke_l )  THEN
 
@@ -4791,5 +4999,5 @@
        !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
        !$ACC PRIVATE(dvar_dz, duv2_dz2, ml_stratification, ml_local_profile, rif) &
-       !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, u, v, ml_wall_adjusted, ml_blackadar) &
+       !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, u, v, delta, ml_blackadar) &
        !$ACC PRESENT(kh, km, sums_l_l, rmask)
        DO  i = nxlg, nxrg
@@ -4839,5 +5047,5 @@
              !DIR$ IVDEP
              DO  k = nzb+1, nzt
-                ml_local_profile(k) = MIN( ml_stratification(k), ml_wall_adjusted(k,j,i) ) &
+                ml_local_profile(k) = MIN( ml_stratification(k), delta(k,j,i) ) &
                                     * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
                 km(k,j,i) = c_0 * ml_local_profile(k) * SQRT( e(k,j,i) )
@@ -4993,5 +5201,5 @@
     REAL(wp)     ::  ldsd        !< sum: ld_ij*sd_ij
 
-    REAL(wp)     ::  delta       !< grid size
+    REAL(wp)     ::  delta_max   !< maximum of the grid spacings
     REAL(wp)     ::  cst         !< c*
     REAL(wp)     ::  cstnust_t   !< product c*nu*
@@ -5119,5 +5327,5 @@
 !--          The model was only tested for an isotropic grid. The following
 !--          expression was a recommendation of Stefan Heinz.
-             delta = MAX( dx, dy, dzw(k) )
+             delta_max = MAX( dx, dy, dzw(k) )
 
              IF ( lnn <= 0.0_wp ) THEN
@@ -5125,5 +5333,5 @@
              ELSE
                 cst = cstnust_t /                                              &
-                   ( 4.0_wp * delta * SQRT( lnn / 2.0_wp ) + 1.0E-20_wp )
+                   ( 4.0_wp * delta_max * SQRT( lnn / 2.0_wp ) + 1.0E-20_wp )
              ENDIF
 
@@ -5131,5 +5339,5 @@
 !--          Calculate border according to Mokhtarpoor and Heinz (2017)
              cst_max = fac_cmax * SQRT( e(k,j,i) ) /                           &
-                       ( delta * SQRT( 2.0_wp * sd2 ) + 1.0E-20_wp )
+                       ( delta_max * SQRT( 2.0_wp * sd2 ) + 1.0E-20_wp )
 
              IF ( ABS( cst ) > cst_max )  THEN
@@ -5137,5 +5345,5 @@
              ENDIF
 
-             km(k,j,i) = cst * delta * SQRT( e(k,j,i) ) * flag
+             km(k,j,i) = cst * delta_max * SQRT( e(k,j,i) ) * flag
 
           ENDDO
