Index: /palm/trunk/SOURCE/header.f90
===================================================================
--- /palm/trunk/SOURCE/header.f90	(revision 4472)
+++ /palm/trunk/SOURCE/header.f90	(revision 4473)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! revised message if wall_adjustment is used
+! 
+! 4444 2020-03-05 15:59:50Z raasch
 ! bugfix: cpp-directives for serial mode added
 ! 
@@ -1539,5 +1542,5 @@
        IF ( e_init > 0.0_wp )  WRITE ( io, 455 )  e_init
        IF ( e_min > 0.0_wp )  WRITE ( io, 454 )  e_min
-       IF ( wall_adjustment )  WRITE ( io, 453 )  wall_adjustment_factor
+       IF ( wall_adjustment )  WRITE ( io, 453 )
     ENDIF
     IF ( rans_mode )  THEN
@@ -1884,5 +1887,5 @@
 451 FORMAT ('    Diffusion coefficients are constant:'/ &
             '    Km = ',F6.2,' m**2/s   Kh = ',F6.2,' m**2/s   Pr = ',F5.2)
-453 FORMAT ('    Mixing length is limited to',F5.2,' * z')
+453 FORMAT ('    Mixing length is limited close to surfaces')
 454 FORMAT ('    TKE is not allowed to fall below ',E9.2,' (m/s)**2')
 455 FORMAT ('    initial TKE is prescribed as ',E9.2,' (m/s)**2')
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 4472)
+++ /palm/trunk/SOURCE/modules.f90	(revision 4473)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! moved wall_adjustment_factor to turbulence_closure_mod
+!
+! 4472 2020-03-24 12:21:00Z Giersch
 ! Additional switch added to activate calculations in flow_statistics for the
 ! kolmogorov length scale
@@ -961,5 +964,4 @@
     REAL(wp) ::  v_bulk = 0.0_wp                               !< namelist parameter
     REAL(wp) ::  v_gtrans = 0.0_wp                             !< transformed wind component (galilei transformation)
-    REAL(wp) ::  wall_adjustment_factor = 1.8_wp               !< adjustment factor for mixing length l
     REAL(wp) ::  zeta_max = 20.0_wp                            !< namelist parameter
     REAL(wp) ::  zeta_min = -20.0_wp                           !< namelist parameter
Index: /palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- /palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 4472)
+++ /palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 4473)
@@ -25,6 +25,19 @@
 ! -----------------
 ! $Id$
+! - rename l-grid to gridsize-geometric-mean
+!          l-wall to ml-wall-adjusted
+!          l-stable to ml-stratification
+!          l-black to ml-blackadar
+!          l-v to ml-local-profile
+!          l-max to max-length-scale
+!          l to ml
+! - adjust some comments
+! - corrected definition of wall-adjusted mixing length to include
+!   gridsize-geometric-mean
+! - moved definition of wall_adjustment_factor to this module
+!
+! 4457 2020-03-11 14:20:43Z raasch
 ! use statement for exchange horiz added
-! 
+!
 ! 4433 2020-02-28 22:14:43Z gronemeier
 ! remove warning for newly implemented RANS mode
@@ -97,5 +110,4 @@
 ! ------------
 !> This module contains the available turbulence closures for PALM.
-!>
 !>
 !> @todo test initialization for all possibilities
@@ -181,24 +193,22 @@
 
 
-    REAL(wp) ::  c_0                !< constant used for diffusion coefficient and dissipation (dependent on mode RANS/LES)
-    REAL(wp) ::  c_1                !< model constant for RANS mode
-    REAL(wp) ::  c_2                !< model constant for RANS mode
-    REAL(wp) ::  c_3                !< model constant for RANS mode
-    REAL(wp) ::  c_4                !< model constant for RANS mode
-    REAL(wp) ::  l_max              !< maximum length scale for Blackadar mixing length
-    REAL(wp) ::  dsig_e = 1.0_wp    !< factor to calculate Ke from Km (1/sigma_e)
-    REAL(wp) ::  dsig_diss = 1.0_wp !< factor to calculate K_diss from Km (1/sigma_diss)
-
-    REAL(wp), DIMENSION(0:4) :: rans_const_c = &       !< model constants for RANS mode (namelist param)
+    REAL(wp) ::  c_0                 !< constant used for diffusion coefficient and dissipation (dependent on mode RANS/LES)
+    REAL(wp) ::  c_1                 !< model constant for RANS mode
+    REAL(wp) ::  c_2                 !< model constant for RANS mode
+    REAL(wp) ::  c_3                 !< model constant for RANS mode
+    REAL(wp) ::  c_4                 !< model constant for RANS mode
+    REAL(wp) ::  dsig_e = 1.0_wp     !< factor to calculate Ke from Km (1/sigma_e)
+    REAL(wp) ::  dsig_diss = 1.0_wp  !< factor to calculate K_diss from Km (1/sigma_diss)
+    REAL(wp) ::  length_scale_max    !< maximum length scale for Blackadar mixing length
+    REAL(wp) ::  wall_adjustment_factor = 1.8_wp  !< adjustment factor for mixing length
+
+    REAL(wp), DIMENSION(0:4) :: rans_const_c = &        !< model constants for RANS mode (namelist param)
        (/ 0.55_wp, 1.44_wp, 1.92_wp, 1.44_wp, 0.0_wp /) !> default values fit for standard-tke-e closure
 
-    REAL(wp), DIMENSION(2) :: rans_const_sigma = &     !< model constants for RANS mode, sigma values (sigma_e, sigma_diss) (namelist param)
-       (/ 1.0_wp, 1.30_wp /)
-
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l_black    !< mixing length according to Blackadar
-
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l_grid     !< geometric mean of grid sizes dx, dy, dz
-
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  l_wall !< near-wall mixing length
+    REAL(wp), DIMENSION(2) :: rans_const_sigma = &      !< model constants for RANS mode, sigma values (namelist param)
+       (/ 1.0_wp, 1.30_wp /)                            !> (sigma_e, sigma_diss)
+
+    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
 
 !
@@ -1147,5 +1157,6 @@
                    DO  j = nysg, nyng
                       DO  k = nzb+1, nzt
-                         km(k,j,i) = c_0 * l_wall(k,j,i) * SQRT( e_init )
+                         km(k,j,i) = c_0 * SQRT( e_init ) * MIN( ml_wall_adjusted(k,j,i), &
+                                                                 ml_blackadar(k) )
                       ENDDO
                    ENDDO
@@ -1191,5 +1202,5 @@
                 DO  j = nysg, nyng
                    DO  k = nzb+1, nzt
-                      km(k,j,i) = c_0 * l_wall(k,j,i) * SQRT( e_init )
+                      km(k,j,i) = c_0 * SQRT( e_init ) * ml_wall_adjusted(k,j,i)
                    ENDDO
                 ENDDO
@@ -1368,5 +1379,5 @@
 
     USE control_parameters,                                                    &
-        ONLY:  f, message_string, wall_adjustment, wall_adjustment_factor
+        ONLY:  f, message_string, wall_adjustment
 
     USE exchange_horiz_mod,                                                    &
@@ -1406,45 +1417,51 @@
     INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE :: vicinity !< contains topography information of the vicinity of (i/j/k)
 
-    REAL(wp) :: radius          !< search radius in meter
-
-    ALLOCATE( l_grid(1:nzt) )
-    ALLOCATE( l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+    REAL(wp) :: distance_up            !< distance of grid box center to its boundary in upper direction
+    REAL(wp) :: distance_down          !< distance of grid box center to its boundary in lower direction
+    REAL(wp) :: distance_ns            !< distance of grid box center to its boundary in y direction
+    REAL(wp) :: distance_lr            !< distance of grid box center to its boundary in x direction
+    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_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
+    REAL(wp) :: distance_corners_down  !< distance of grid box center to its upper corners
+    REAL(wp) :: distance_corners_up    !< distance of grid box center to its lower corners
+    REAL(wp) :: radius                 !< search radius in meter
+
+    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) )
 !
 !-- Initialize the mixing length in case of an LES-simulation
     IF ( .NOT. rans_mode )  THEN
 !
-!--    Compute the grid-dependent mixing length.
-       DO  k = 1, nzt
-          l_grid(k)  = ( dx * dy * dzw(k) )**0.33333333333333_wp
+!--    Compute the geometric averaged grid size (used for limiting the mixing length)
+       DO  k = nzb+1, nzt
+          gridsize_geometric_mean(k)  = ( dx * dy * dzw(k) )**0.33333333333333_wp
        ENDDO
-!
-!--    Initialize near-wall mixing length l_wall only in the vertical direction
-!--    for the moment, multiplication with wall_adjustment_factor further below
-       l_wall(nzb,:,:)   = l_grid(1)
-       DO  k = nzb+1, nzt
-          l_wall(k,:,:)  = l_grid(k)
+       gridsize_geometric_mean(nzb)   = gridsize_geometric_mean(1)
+       gridsize_geometric_mean(nzt+1) = gridsize_geometric_mean(nzt)
+
+       IF ( ANY( gridsize_geometric_mean > 1.5_wp * dx * wall_adjustment_factor ) .OR. &
+            ANY( gridsize_geometric_mean > 1.5_wp * dy * wall_adjustment_factor ) )  THEN
+            WRITE( message_string, * ) 'grid anisotropy exceeds threshold',    &
+                                      ' &starting from height level k = ', k, &
+                                      '.'
+           CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 )
+       ENDIF
+
+!
+!--    Initialize near-wall mixing length ml_wall_adjusted
+       DO  k = nzb, nzt+1
+          ml_wall_adjusted(k,:,:) = gridsize_geometric_mean(k)
        ENDDO
-       l_wall(nzt+1,:,:) = l_grid(nzt)
 
        IF ( wall_adjustment )  THEN
-
-          DO  k = 1, nzt
-             IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR.            &
-                  l_grid(k) > 1.5_wp * dy * wall_adjustment_factor )  THEN
-                WRITE( message_string, * ) 'grid anisotropy exceeds ',             &
-                                           'threshold given by only local',        &
-                                           ' &horizontal reduction of near_wall ', &
-                                           'mixing length l_wall',                 &
-                                           ' &starting from height level k = ', k, &
-                                           '.'
-                CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 )
-                EXIT
-             ENDIF
-          ENDDO
-!
-!--       In case of topography: limit near-wall mixing length l_wall further:
-!--       Go through all points of the subdomain one by one and look for the closest
-!--       surface.
-!--       Is this correct in the ocean case?
+!
+!--       In case of topography, adjust mixing length if there is any wall at
+!--       the surrounding grid boxes:
+          !> @todo check if this is correct also for the ocean case
           DO  i = nxl, nxr
              DO  j = nys, nyn
@@ -1454,81 +1471,115 @@
                    IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )  THEN
 !
-!--                   Check for neighbouring grid-points.
-!--                   Vertical distance, down
-                      IF ( .NOT. BTEST( wall_flags_total_0(k-1,j,i), 0 ) )             &
-                         l_wall(k,j,i) = MIN( l_grid(k), zu(k) - zw(k-1) )
-!
-!--                   Vertical distance, up
-                      IF ( .NOT. BTEST( wall_flags_total_0(k+1,j,i), 0 ) )             &
-                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), zw(k) - zu(k) )
-!
-!--                   y-distance
-                      IF ( .NOT. BTEST( wall_flags_total_0(k,j-1,i), 0 )  .OR.         &
-                           .NOT. BTEST( wall_flags_total_0(k,j+1,i), 0 ) )             &
-                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dy )
-!
-!--                   x-distance
-                      IF ( .NOT. BTEST( wall_flags_total_0(k,j,i-1), 0 )  .OR.         &
-                           .NOT. BTEST( wall_flags_total_0(k,j,i+1), 0 ) )             &
-                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dx )
-!
-!--                   yz-distance (vertical edges, down)
-                      IF ( .NOT. BTEST( wall_flags_total_0(k-1,j-1,i), 0 )  .OR.       &
-                           .NOT. BTEST( wall_flags_total_0(k-1,j+1,i), 0 )  )          &
-                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),                &
-                                              SQRT( 0.25_wp * dy**2 +                  &
-                                             ( zu(k) - zw(k-1) )**2 ) )
-!
-!--                   yz-distance (vertical edges, up)
-                      IF ( .NOT. BTEST( wall_flags_total_0(k+1,j-1,i), 0 )  .OR.       &
-                           .NOT. BTEST( wall_flags_total_0(k+1,j+1,i), 0 )  )          &
-                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),                &
-                                              SQRT( 0.25_wp * dy**2 +                  &
-                                             ( zw(k) - zu(k) )**2 ) )
-!
-!--                   xz-distance (vertical edges, down)
-                      IF ( .NOT. BTEST( wall_flags_total_0(k-1,j,i-1), 0 )  .OR.       &
-                           .NOT. BTEST( wall_flags_total_0(k-1,j,i+1), 0 )  )          &
-                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),                &
-                                              SQRT( 0.25_wp * dx**2 +                  &
-                                             ( zu(k) - zw(k-1) )**2 ) )
-!
-!--                   xz-distance (vertical edges, up)
-                      IF ( .NOT. BTEST( wall_flags_total_0(k+1,j,i-1), 0 )  .OR.       &
-                           .NOT. BTEST( wall_flags_total_0(k+1,j,i+1), 0 )  )          &
-                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),                &
-                                              SQRT( 0.25_wp * dx**2 +                  &
-                                             ( zw(k) - zu(k) )**2 ) )
-!
-!--                   xy-distance (horizontal edges)
-                      IF ( .NOT. BTEST( wall_flags_total_0(k,j-1,i-1), 0 )  .OR.        &
-                           .NOT. BTEST( wall_flags_total_0(k,j+1,i-1), 0 )  .OR.        &
-                           .NOT. BTEST( wall_flags_total_0(k,j-1,i+1), 0 )  .OR.        &
-                           .NOT. BTEST( wall_flags_total_0(k,j+1,i+1), 0 ) )            &
-                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),                 &
-                                              SQRT( 0.25_wp * ( dx**2 + dy**2 ) ) )
-!
-!--                   xyz distance (vertical and horizontal edges, down)
-                      IF ( .NOT. BTEST( wall_flags_total_0(k-1,j-1,i-1), 0 )  .OR.      &
-                           .NOT. BTEST( wall_flags_total_0(k-1,j+1,i-1), 0 )  .OR.      &
-                           .NOT. BTEST( wall_flags_total_0(k-1,j-1,i+1), 0 )  .OR.      &
-                           .NOT. BTEST( wall_flags_total_0(k-1,j+1,i+1), 0 ) )          &
-                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),                 &
-                                              SQRT( 0.25_wp * ( dx**2 + dy**2 )         &
-                                                    +  ( zu(k) - zw(k-1) )**2  ) )
-!
-!--                   xyz distance (vertical and horizontal edges, up)
-                      IF ( .NOT. BTEST( wall_flags_total_0(k+1,j-1,i-1), 0 )  .OR.      &
-                           .NOT. BTEST( wall_flags_total_0(k+1,j+1,i-1), 0 )  .OR.      &
-                           .NOT. BTEST( wall_flags_total_0(k+1,j-1,i+1), 0 )  .OR.      &
-                           .NOT. BTEST( wall_flags_total_0(k+1,j+1,i+1), 0 ) )          &
-                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),                 &
-                                              SQRT( 0.25_wp * ( dx**2 + dy**2 )         &
-                                                    +  ( zw(k) - zu(k) )**2  ) )
-
-                   ENDIF
-!
-!--                Adjust mixing length by wall-adjustment factor and limit it by l_grid
-                   l_wall(k,j,i) = MIN( l_wall(k,j,i) * wall_adjustment_factor, l_grid(k) )
+!--                   First, check if grid points directly next to current grid point
+!--                   are surface grid points
+!--                   Check along...
+!--                   ...vertical direction, down
+                      IF ( .NOT. BTEST( wall_flags_total_0(k-1,j,i), 0 ) )  THEN
+                         distance_down =  zu(k) - zw(k-1)
+                      ELSE
+                         distance_down = 9999999.9_wp
+                      ENDIF
+!
+!--                   ...vertical direction, up
+                      IF ( .NOT. BTEST( wall_flags_total_0(k+1,j,i), 0 ) )  THEN
+                         distance_up = zw(k) - zu(k)
+                      ELSE
+                         distance_up = 9999999.9_wp
+                      ENDIF
+!
+!--                   ...y-direction (north/south)
+                      IF ( .NOT. BTEST( wall_flags_total_0(k,j-1,i), 0 )  .OR. &
+                           .NOT. BTEST( wall_flags_total_0(k,j+1,i), 0 ) )  THEN
+                         distance_ns = 0.5_wp * dy
+                      ELSE
+                         distance_ns = 9999999.9_wp
+                      ENDIF
+!
+!--                   ...x-direction (left/right)
+                      IF ( .NOT. BTEST( wall_flags_total_0(k,j,i-1), 0 )  .OR. &
+                           .NOT. BTEST( wall_flags_total_0(k,j,i+1), 0 ) )  THEN
+                         distance_lr = 0.5_wp * dx
+                      ELSE
+                         distance_lr = 9999999.9_wp
+                      ENDIF
+!
+!--                   Now, check the edges along...
+!--                   ...yz-direction (vertical edges, down)
+                      IF ( .NOT. BTEST( wall_flags_total_0(k-1,j-1,i), 0 )  .OR.  &
+                           .NOT. BTEST( wall_flags_total_0(k-1,j+1,i), 0 )  )  THEN
+                         distance_edge_yz_down = SQRT( 0.25_wp * dy**2 + ( zu(k) - zw(k-1) )**2 )
+                      ELSE
+                         distance_edge_yz_down = 9999999.9_wp
+                      ENDIF
+!
+!--                   ...yz-direction (vertical edges, up)
+                      IF ( .NOT. BTEST( wall_flags_total_0(k+1,j-1,i), 0 )  .OR.  &
+                           .NOT. BTEST( wall_flags_total_0(k+1,j+1,i), 0 )  )  THEN
+                         distance_edge_yz_up = SQRT( 0.25_wp * dy**2 + ( zw(k) - zu(k) )**2 )
+                      ELSE
+                         distance_edge_yz_up = 9999999.9_wp
+                      ENDIF
+!
+!--                   ...xz-direction (vertical edges, down)
+                      IF ( .NOT. BTEST( wall_flags_total_0(k-1,j,i-1), 0 )  .OR.  &
+                           .NOT. BTEST( wall_flags_total_0(k-1,j,i+1), 0 )  )  THEN
+                         distance_edge_xz_down = SQRT( 0.25_wp * dx**2 + ( zu(k) - zw(k-1) )**2 )
+                      ELSE
+                         distance_edge_xz_down = 9999999.9_wp
+                      ENDIF
+!
+!--                   ...xz-direction (vertical edges, up)
+                      IF ( .NOT. BTEST( wall_flags_total_0(k+1,j,i-1), 0 )  .OR.  &
+                           .NOT. BTEST( wall_flags_total_0(k+1,j,i+1), 0 )  )  THEN
+                         distance_edge_xz_up = SQRT( 0.25_wp * dx**2 + ( zw(k) - zu(k) )**2 )
+                      ELSE
+                         distance_edge_xz_up = 9999999.9_wp
+                      ENDIF
+!
+!--                   ...xy-direction (horizontal edges)
+                      IF ( .NOT. BTEST( wall_flags_total_0(k,j-1,i-1), 0 )  .OR. &
+                           .NOT. BTEST( wall_flags_total_0(k,j+1,i-1), 0 )  .OR. &
+                           .NOT. BTEST( wall_flags_total_0(k,j-1,i+1), 0 )  .OR. &
+                           .NOT. BTEST( wall_flags_total_0(k,j+1,i+1), 0 ) )  THEN
+                         distance_edge_xy = SQRT( 0.25_wp * ( dx**2 + dy**2 ) )
+                      ELSE
+                         distance_edge_xy = 9999999.9_wp
+                      ENDIF
+!
+!--                   Now, check the corners...
+!--                   ...lower four corners
+                      IF ( .NOT. BTEST( wall_flags_total_0(k-1,j-1,i-1), 0 )  .OR. &
+                           .NOT. BTEST( wall_flags_total_0(k-1,j+1,i-1), 0 )  .OR. &
+                           .NOT. BTEST( wall_flags_total_0(k-1,j-1,i+1), 0 )  .OR. &
+                           .NOT. BTEST( wall_flags_total_0(k-1,j+1,i+1), 0 ) )  THEN
+                         distance_corners_down = SQRT( 0.25_wp * ( dx**2 + dy**2 ) &
+                                                       + ( zu(k) - zw(k-1) )**2 )
+                      ELSE
+                         distance_corners_down = 9999999.9_wp
+                      ENDIF
+!
+!--                   ...upper four corners
+                      IF ( .NOT. BTEST( wall_flags_total_0(k+1,j-1,i-1), 0 )  .OR. &
+                           .NOT. BTEST( wall_flags_total_0(k+1,j+1,i-1), 0 )  .OR. &
+                           .NOT. BTEST( wall_flags_total_0(k+1,j-1,i+1), 0 )  .OR. &
+                           .NOT. BTEST( wall_flags_total_0(k+1,j+1,i+1), 0 ) )  THEN
+                         distance_corners_up = SQRT( 0.25_wp * ( dx**2 + dy**2 )   &
+                                                     + ( zw(k) - zu(k) )**2 )
+                      ELSE
+                         distance_corners_up = 9999999.9_wp
+                      ENDIF
+!
+!--                   Adjust mixing length by wall-adjustment factor
+                      ml_wall_adjusted(k,j,i) = wall_adjustment_factor * MIN(  &
+                         distance_up, distance_down, distance_ns, distance_lr, &
+                         distance_edge_yz_down, distance_edge_yz_up,           &
+                         distance_edge_xz_down, distance_edge_xz_up,           &
+                         distance_edge_xy,                                     &
+                         distance_corners_down, distance_corners_up )
+
+                  ENDIF  !if grid point belongs to atmosphere
+
+                  ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i), &
+                                                 gridsize_geometric_mean(k) )
 
                 ENDDO  !k loop
@@ -1538,23 +1589,23 @@
        ENDIF  !if wall_adjustment
 
-    ELSE
+    ELSE  !--> RANS mode
 !
 !--    Initialize the mixing length in case of a RANS simulation
-       ALLOCATE( l_black(nzb:nzt+1) )
+       ALLOCATE( ml_blackadar(nzb:nzt+1) )
 
 !
 !--    Calculate mixing length according to Blackadar (1962)
        IF ( f /= 0.0_wp )  THEN
-          l_max = 2.7E-4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) /            &
-                  ABS( f ) + 1.0E-10_wp
+          length_scale_max = 2.7E-4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 )   &
+                           / ABS( f ) + 1.0E-10_wp
        ELSE
-          l_max = 30.0_wp
+          length_scale_max = 30.0_wp
        ENDIF
 
        DO  k = nzb, nzt
-          l_black(k) = kappa * zu(k) / ( 1.0_wp + kappa * zu(k) / l_max )
+          ml_blackadar(k) = kappa * zu(k) / ( 1.0_wp + kappa * zu(k) / length_scale_max )
        ENDDO
 
-       l_black(nzt+1) = l_black(nzt)
+       ml_blackadar(nzt+1) = ml_blackadar(nzt)
 
 !
@@ -1571,6 +1622,6 @@
        ENDDO
 
-       l_wall(nzb,:,:) = l_black(nzb)
-       l_wall(nzt+1,:,:) = l_black(nzt+1)
+       ml_wall_adjusted(nzb,:,:) = ml_blackadar(nzb)
+       ml_wall_adjusted(nzt+1,:,:) = ml_blackadar(nzt+1)
 !
 !--    Limit mixing length to either nearest wall or Blackadar mixing length.
@@ -1581,8 +1632,8 @@
        DO  k = nzb+1, nzt
 
-          radius = l_black(k)  ! radius within walls are searched
-!
-!--       Set l_wall to its default maximum value (l_back)
-          l_wall(k,:,:) = radius
+          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
 
 !
@@ -1679,5 +1730,5 @@
 !--                            Search for walls within octant (+++)
                                vic_yz = vicinity(0:rad_k+1,0:rad_j,ii)
-                               l_wall(k,j,i) = MIN( l_wall(k,j,i),             &
+                               ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),             &
                                        shortest_distance( vic_yz, .TRUE., ii ) )
 !
@@ -1687,5 +1738,5 @@
 !--                            shortest_distance").
                                vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii)
-                               l_wall(k,j,i) = MIN( l_wall(k,j,i),             &
+                               ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),             &
                                        shortest_distance( vic_yz, .TRUE., ii ) )
 
@@ -1694,14 +1745,14 @@
 !--                         Search for walls within octant (+--)
                             vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii)
-                            l_wall(k,j,i) = MIN( l_wall(k,j,i),                &
+                            ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(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)
-                            l_wall(k,j,i) = MIN( l_wall(k,j,i),                &
+                            ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),                &
                                       shortest_distance( vic_yz, .FALSE., ii ) )
 !
 !--                         Reduce search along x by already found distance
-                            dist_dx = CEILING( l_wall(k,j,i) / dx )
+                            dist_dx = CEILING( ml_wall_adjusted(k,j,i) / dx )
 
                          ENDDO
@@ -1716,5 +1767,5 @@
 !--                            Search for walls within octant (-++)
                                vic_yz = vicinity(0:rad_k+1,0:rad_j,ii)
-                               l_wall(k,j,i) = MIN( l_wall(k,j,i),             &
+                               ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),             &
                                       shortest_distance( vic_yz, .TRUE., -ii ) )
 !
@@ -1724,5 +1775,5 @@
 !--                            shortest_distance").
                                vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii)
-                               l_wall(k,j,i) = MIN( l_wall(k,j,i),             &
+                               ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),             &
                                       shortest_distance( vic_yz, .TRUE., -ii ) )
 
@@ -1731,14 +1782,14 @@
 !--                         Search for walls within octant (---)
                             vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii)
-                            l_wall(k,j,i) = MIN( l_wall(k,j,i),                &
+                            ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(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)
-                            l_wall(k,j,i) = MIN( l_wall(k,j,i),                &
+                            ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),                &
                                      shortest_distance( vic_yz, .FALSE., -ii ) )
 !
 !--                         Reduce search along x by already found distance
-                            dist_dx = CEILING( l_wall(k,j,i) / dx )
+                            dist_dx = CEILING( ml_wall_adjusted(k,j,i) / dx )
 
                          ENDDO
@@ -1748,5 +1799,5 @@
                    ELSE  !Check if (i,j,k) belongs to atmosphere
 
-                      l_wall(k,j,i) = l_black(k)
+                      ml_wall_adjusted(k,j,i) = ml_blackadar(k)
 
                    ENDIF
@@ -1762,14 +1813,13 @@
        ENDDO  !k loop
 
-       !$ACC ENTER DATA COPYIN(l_black(nzb:nzt+1))
+       !$ACC ENTER DATA COPYIN(ml_blackadar(nzb:nzt+1))
 
     ENDIF  !LES or RANS mode
 
 !
-!-- Set lateral boundary conditions for l_wall
-    CALL exchange_horiz( l_wall, nbgp )
-
-    !$ACC ENTER DATA COPYIN(l_grid(nzb:nzt+1)) &
-    !$ACC COPYIN(l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
+!-- 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))
 
     CONTAINS
@@ -2266,5 +2316,6 @@
        IF ( plant_canopy )  CALL pcm_tendency( 6 )
 
-!       CALL user_actions( 'e-tendency' ) ToDo: find general solution for circular dependency between modules
+       !> @todo find general solution for circular dependency between modules
+       ! CALL user_actions( 'e-tendency' )
 
 !
@@ -2377,7 +2428,7 @@
 !
 !--    Additional sink term for flows through plant canopies
-!        IF ( plant_canopy )  CALL pcm_tendency( ? )         !> @todo not yet implemented
-
-!       CALL user_actions( 'e-tendency' ) ToDo: find general solution for circular dependency between modules
+       ! IF ( plant_canopy )  CALL pcm_tendency( ? )     !> @todo not yet implemented
+
+       ! CALL user_actions( 'e-tendency' ) !> @todo find general solution for circular dependency between modules
 
 !
@@ -2496,5 +2547,5 @@
        IF ( plant_canopy )  CALL pcm_tendency( i, j, 6 )
 
-!       CALL user_actions( i, j, 'e-tendency' ) ToDo: find general solution for circular dependency between modules
+       ! CALL user_actions( i, j, 'e-tendency' ) !> @todo: find general solution for circular dependency between modules
 
 !
@@ -2563,5 +2614,5 @@
 !        IF ( plant_canopy )  CALL pcm_tendency( i, j, ? )     !> @todo not yet implemented
 
-!       CALL user_actions( i, j, 'diss-tendency' ) ToDo: find general solution for circular dependency between modules
+!       CALL user_actions( i, j, 'diss-tendency' ) !> @todo: find general solution for circular dependency between modules
 
 !
@@ -2726,5 +2777,5 @@
 
                    km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
-                                   * 0.5_wp * dy
+                                      * 0.5_wp * dy
 !
 !--                -1.0 for right-facing wall, 1.0 for left-facing wall
@@ -2745,5 +2796,5 @@
 
                    km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
-                                   * 0.5_wp * dy
+                                      * 0.5_wp * dy
 !
 !--                -1.0 for right-facing wall, 1.0 for left-facing wall
@@ -2764,5 +2815,5 @@
 
                    km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
-                                   * 0.5_wp * dy
+                                      * 0.5_wp * dy
 !
 !--                -1.0 for right-facing wall, 1.0 for left-facing wall
@@ -3410,5 +3461,5 @@
 
              km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
-                             * 0.5_wp * dy
+                                * 0.5_wp * dy
 !
 !--          -1.0 for right-facing wall, 1.0 for left-facing wall
@@ -3428,5 +3479,5 @@
 
              km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
-                             * 0.5_wp * dy
+                                * 0.5_wp * dy
 !
 !--          -1.0 for right-facing wall, 1.0 for left-facing wall
@@ -3446,5 +3497,5 @@
 
              km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
-                             * 0.5_wp * dy
+                                * 0.5_wp * dy
 !
 !--          -1.0 for right-facing wall, 1.0 for left-facing wall
@@ -3911,5 +3962,5 @@
                 flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) )
                 tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
-                              MERGE( vpt_reference, vpt(k,j,i),          &
+                              MERGE( vpt_reference, vpt(k,j,i),        &
                                      use_single_reference_value ) )
              ENDDO
@@ -3971,9 +4022,9 @@
     REAL(wp)     ::  duv2_dz2       !< squared vertical gradient of wind vector
     REAL(wp)     ::  dvar_dz        !< vertical gradient of var
-    REAL(wp)     ::  l              !< mixing length
+    REAL(wp)     ::  ml             !< mixing length
     REAL(wp)     ::  var_reference  !< reference temperature
 
-    REAL(wp), DIMENSION(nzb+1:nzt) ::  l_stable  !< mixing length according to stratification
-    REAL(wp), DIMENSION(nzb+1:nzt) ::  rif       !< Richardson flux number
+    REAL(wp), DIMENSION(nzb+1:nzt) ::  ml_stratification  !< mixing length according to stratification
+    REAL(wp), DIMENSION(nzb+1:nzt) ::  rif                !< Richardson flux number
 
     REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  var  !< temperature
@@ -3984,7 +4035,7 @@
 
        !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
-       !$ACC PRIVATE(l, l_stable, dvar_dz) &
+       !$ACC PRIVATE(ml, ml_stratification, dvar_dz) &
        !$ACC PRESENT(diss, e, var, wall_flags_total_0) &
-       !$ACC PRESENT(dd2zu, l_grid, l_wall)
+       !$ACC PRESENT(dd2zu, ml_wall_adjusted)
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -3995,12 +4046,12 @@
                 IF ( dvar_dz > 0.0_wp ) THEN
                    IF ( use_single_reference_value )  THEN
-                     l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )                  &
+                     ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) )         &
                                  / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
                    ELSE
-                     l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )               &
+                     ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) )         &
                                  / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
                    ENDIF
                 ELSE
-                   l_stable(k) = l_grid(k)
+                   ml_stratification(k) = ml_wall_adjusted(k,j,i)
                 ENDIF
 
@@ -4011,8 +4062,8 @@
              DO  k = nzb+1, nzt
 
-                l  = MIN( l_wall(k,j,i), l_stable(k) )
-
-                diss(k,j,i) = ( 0.19_wp + 0.74_wp * l / l_wall(k,j,i) )                &
-                              * e(k,j,i) * SQRT( e(k,j,i) ) / l                        &
+                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) )           &
+                              * 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 ) )
 
@@ -4025,7 +4076,7 @@
 
       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
-      !$ACC PRIVATE(l_stable, duv2_dz2, rif, dvar_dz) &
+      !$ACC PRIVATE(ml_stratification, duv2_dz2, rif, dvar_dz) &
       !$ACC PRESENT(diss, e, u, v, var, wall_flags_total_0) &
-      !$ACC PRESENT(dd2zu, l_black, l_wall)
+      !$ACC PRESENT(dd2zu, ml_blackadar, ml_wall_adjusted)
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -4062,7 +4113,7 @@
              DO  k = nzb+1, nzt
                 IF ( rif(k) >= 0.0_wp )  THEN
-                   l_stable(k) = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif(k) ), l_wall(k,j,i) )
+                   ml_stratification(k) = ml_blackadar(k) / ( 1.0_wp + 5.0_wp * rif(k) )
                 ELSE
-                   l_stable(k) = l_wall(k,j,i) * SQRT( 1.0_wp - 16.0_wp * rif(k) )
+                   ml_stratification(k) = ml_blackadar(k) * SQRT( 1.0_wp - 16.0_wp * rif(k) )
                 ENDIF
              ENDDO
@@ -4072,5 +4123,6 @@
              DO  k = nzb+1, nzt
 
-                diss(k,j,i) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / l_stable(k)     &
+                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) )         &
                             * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
 
@@ -4180,9 +4232,9 @@
     REAL(wp)     ::  duv2_dz2       !< squared vertical gradient of wind vector
     REAL(wp)     ::  dvar_dz        !< vertical gradient of var
-    REAL(wp)     ::  l              !< mixing length
+    REAL(wp)     ::  ml             !< mixing length
     REAL(wp)     ::  var_reference  !< reference temperature
 
-    REAL(wp), DIMENSION(nzb+1:nzt) ::  l_stable  !< mixing length according to stratification
-    REAL(wp), DIMENSION(nzb+1:nzt) ::  rif       !< Richardson flux number
+    REAL(wp), DIMENSION(nzb+1:nzt) ::  ml_stratification  !< mixing length according to stratification
+    REAL(wp), DIMENSION(nzb+1:nzt) ::  rif                !< Richardson flux number
 
     REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  var  !< temperature
@@ -4197,12 +4249,12 @@
           IF ( dvar_dz > 0.0_wp ) THEN
              IF ( use_single_reference_value )  THEN
-                l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )                  &
+                ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) )              &
                             / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
              ELSE
-                l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )               &
+                ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) )              &
                             / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
              ENDIF
           ELSE
-             l_stable(k) = l_grid(k)
+             ml_stratification(k) = ml_wall_adjusted(k,j,i)
           ENDIF
        ENDDO
@@ -4210,8 +4262,8 @@
        !DIR$ IVDEP
        DO  k = nzb+1, nzt
-          l  = MIN( l_wall(k,j,i), l_stable(k) )
-
-          diss(k,j,i) = ( 0.19_wp + 0.74_wp * l / l_wall(k,j,i) )              &
-                        * e(k,j,i) * SQRT( e(k,j,i) ) / l                      &
+          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) )           &
+                        * 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
@@ -4248,7 +4300,7 @@
        DO  k = nzb+1, nzt
           IF ( rif(k) >= 0.0_wp )  THEN
-             l_stable(k) = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif(k) ), l_wall(k,j,i) )
+             ml_stratification(k) = ml_blackadar(k) / ( 1.0_wp + 5.0_wp * rif(k) )
           ELSE
-             l_stable(k) = l_wall(k,j,i) * SQRT( 1.0_wp - 16.0_wp * rif(k) )
+             ml_stratification(k) = ml_blackadar(k) * SQRT( 1.0_wp - 16.0_wp * rif(k) )
           ENDIF
 
@@ -4257,5 +4309,6 @@
        !DIR$ IVDEP
        DO  k = nzb+1, nzt
-          diss(k,j,i) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / l_stable(k)     &
+          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) )         &
                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
        ENDDO
@@ -4652,11 +4705,11 @@
     REAL(wp)     ::  duv2_dz2            !< squared vertical gradient of wind vector
     REAL(wp)     ::  dvar_dz             !< vertical gradient of var
-    REAL(wp)     ::  l                   !< mixing length (single height)
+    REAL(wp)     ::  ml                  !< mixing length (single height)
     REAL(wp)     ::  var_reference       !< reference temperature
 
-    !DIR$ ATTRIBUTES ALIGN:64:: l_v, l_stable, rif
-    REAL(wp), DIMENSION(nzb+1:nzt) ::  l_v       !< mixing length (all heights)
-    REAL(wp), DIMENSION(nzb+1:nzt) ::  l_stable  !< mixing length according to stratification
-    REAL(wp), DIMENSION(nzb+1:nzt) ::  rif       !< Richardson flux number
+    !DIR$ ATTRIBUTES ALIGN:64:: ml_local_profile, ml_stratification, rif
+    REAL(wp), DIMENSION(nzb+1:nzt) ::  ml_local_profile   !< mixing length (all heights)
+    REAL(wp), DIMENSION(nzb+1:nzt) ::  ml_stratification  !< mixing length according to stratification
+    REAL(wp), DIMENSION(nzb+1:nzt) ::  rif                !< Richardson flux number
 
     REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  var  !< temperature
@@ -4674,5 +4727,5 @@
 !
 !-- Compute the turbulent diffusion coefficient for momentum
-    !$OMP PARALLEL PRIVATE (i,j,k,l,sr,tn)
+    !$OMP PARALLEL PRIVATE (i,j,k,ml,sr,tn)
 !$  tn = omp_get_thread_num()
 
@@ -4680,7 +4733,7 @@
        !$OMP DO
        !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
-       !$ACC PRIVATE(dvar_dz, l, l_stable, l_v) &
-       !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, l_wall, l_grid, rmask) &
-       !$ACC PRESENT(kh, km, sums_l_l)
+       !$ACC PRIVATE(dvar_dz, ml, ml_stratification, ml_local_profile) &
+       !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, ml_wall_adjusted) &
+       !$ACC PRESENT(kh, km, sums_l_l, rmask)
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
@@ -4696,12 +4749,12 @@
                 IF ( dvar_dz > 0.0_wp ) THEN
                    IF ( use_single_reference_value )  THEN
-                      l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )                  &
+                      ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) )         &
                                   / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
                    ELSE
-                      l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )               &
+                      ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) )         &
                                   / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
                    ENDIF
                 ELSE
-                   l_stable(k) = l_grid(k)
+                   ml_stratification(k) = ml_wall_adjusted(k,j,i)
                 ENDIF
 
@@ -4712,11 +4765,11 @@
              DO  k = nzb+1, nzt
 
-                l_v(k) = MIN( l_wall(k,j,i), l_stable(k) )                      &
-                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
-                l = l_v(k)
+                ml = MIN( ml_wall_adjusted(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
 !
 !--             Compute diffusion coefficients for momentum and heat
-                km(k,j,i) = c_0 * l * SQRT( e(k,j,i) )
-                kh(k,j,i) = ( 1.0_wp + 2.0_wp * l / l_wall(k,j,i) ) * km(k,j,i)
+                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)
 
              ENDDO
@@ -4726,5 +4779,5 @@
              DO  sr = 0, statistic_regions
                 DO  k = nzb+1, nzt
-                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l_v(k) * rmask(j,i,sr)
+                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + ml_local_profile(k) * rmask(j,i,sr)
                 ENDDO
              ENDDO
@@ -4737,7 +4790,7 @@
        !$OMP DO
        !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
-       !$ACC PRIVATE(dvar_dz, duv2_dz2, l_stable, l_v, rif) &
-       !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, u, v, l_wall, l_black, rmask) &
-       !$ACC PRESENT(kh, km, sums_l_l)
+       !$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(kh, km, sums_l_l, rmask)
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
@@ -4775,7 +4828,7 @@
              DO  k = nzb+1, nzt
                 IF ( rif(k) >= 0.0_wp )  THEN
-                   l_stable(k)  = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif(k) ), l_wall(k,j,i) )
+                   ml_stratification(k) = ml_blackadar(k) / ( 1.0_wp + 5.0_wp * rif(k) )
                 ELSE
-                   l_stable(k)  = l_wall(k,j,i)
+                   ml_stratification(k) = ml_blackadar(k)
                 ENDIF
 
@@ -4786,6 +4839,7 @@
              !DIR$ IVDEP
              DO  k = nzb+1, nzt
-                l_v(k)    = l_stable(k) * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
-                km(k,j,i) = c_0 * l_v(k) * SQRT( e(k,j,i) )
+                ml_local_profile(k) = MIN( ml_stratification(k), ml_wall_adjusted(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) )
                 kh(k,j,i) = km(k,j,i) / prandtl_number
              ENDDO
@@ -4795,5 +4849,5 @@
              DO  sr = 0, statistic_regions
                 DO  k = nzb+1, nzt
-                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l_v(k) * rmask(j,i,sr)
+                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + ml_local_profile(k) * rmask(j,i,sr)
                 ENDDO
              ENDDO
@@ -4806,5 +4860,5 @@
        !$OMP DO
        !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
-       !$ACC PRIVATE(l_v) &
+       !$ACC PRIVATE(ml_local_profile) &
        !$ACC PRESENT(wall_flags_total_0, e, diss, rmask) &
        !$ACC PRESENT(kh, km, sums_l_l)
@@ -4817,8 +4871,9 @@
              DO  k = nzb+1, nzt
 
-                l_v(k) = c_0**3 * e(k,j,i) * SQRT(e(k,j,i)) / ( diss(k,j,i) + 1.0E-30_wp ) &
-                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
-
-                km(k,j,i) = c_0 * SQRT( e(k,j,i) ) * l_v(k)
+                ml_local_profile(k) = c_0**3 * e(k,j,i) * SQRT(e(k,j,i))                           &
+                                    / ( diss(k,j,i) + 1.0E-30_wp )                                 &
+                                    * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+
+                km(k,j,i) = c_0 * SQRT( e(k,j,i) ) * ml_local_profile(k)
                 kh(k,j,i) = km(k,j,i) / prandtl_number
 
@@ -4829,5 +4884,5 @@
              DO  sr = 0, statistic_regions
                 DO  k = nzb+1, nzt
-                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l_v(k) * rmask(j,i,sr)
+                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + ml_local_profile(k) * rmask(j,i,sr)
                 ENDDO
              ENDDO
