Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 4169)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 4170)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! copy diss, diss_p, tdiss_m to GPU
+! 
+! 4144 2019-08-06 09:11:47Z raasch
 ! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
 ! 
@@ -718,6 +721,6 @@
         ONLY:  d, dd2zu, ddzu, ddzw, drho_air, drho_air_zw, dzw, heatflux_output_conversion, kh,   &
                km, momentumflux_output_conversion, p, ptdf_x, ptdf_y, rdf, rdf_sc, rho_air,        &
-               rho_air_zw, te_m, tpt_m, tu_m, tv_m, tw_m, ug, u_init, u_stokes_zu, vg, v_init,     &
-               v_stokes_zu, zu
+               rho_air_zw, tdiss_m, te_m, tpt_m, tu_m, tv_m, tw_m, ug, u_init, u_stokes_zu, vg,    &
+               v_init, v_stokes_zu, zu
 
     USE control_parameters,                                                                        &
@@ -761,4 +764,5 @@
 !$ACC DATA & 
 !$ACC COPY(d(nzb+1:nzt,nys:nyn,nxl:nxr)) &
+!$ACC COPY(diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
 !$ACC COPY(e(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
 !$ACC COPY(u(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
@@ -771,4 +775,5 @@
 
 !$ACC DATA &
+!$ACC COPY(diss_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
 !$ACC COPY(e_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
 !$ACC COPY(u_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
@@ -777,4 +782,5 @@
 !$ACC COPY(pt_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
 !$ACC COPY(tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+!$ACC COPY(tdiss_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
 !$ACC COPY(te_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
 !$ACC COPY(tu_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
Index: palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 4169)
+++ palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 4170)
@@ -25,4 +25,11 @@
 ! -----------------
 ! $Id$
+! - add performance optimizations according to K. Ketelsen
+!   to diffusion_e and tcm_diffusivities_default
+! - bugfix in calculating l_wall for vertical walls
+! - bugfix in using l_wall in initialization (consider wall_adjustment_factor)
+! - always initialize diss and save the dissipation to that array
+! 
+! 4168 2019-08-16 13:50:17Z suehring
 ! Replace function get_topography_top_index by topo_top_ind
 ! 
@@ -205,5 +212,5 @@
 !>
 !> @todo test initialization for all possibilities
-!>       add OpenMP directives whereever possible
+!> @todo add OpenMP directives whereever possible
 !> @todo Check for random disturbances
 !> @note <Enter notes on the module>
@@ -657,8 +664,8 @@
        ENDDO
 
-       
        diss_p(nzt+1,:,:) = diss_p(nzt,:,:)
-     
+
     ENDIF
+
  END SUBROUTINE tcm_boundary_conds 
  
@@ -1160,7 +1167,4 @@
         ONLY:  collision_turbulence
 
-    USE particle_attributes,                                                   &
-        ONLY:  use_sgs_for_particles, wang_kernel
-
     USE pmc_interface,                                                         &
         ONLY:  nested_run
@@ -1180,13 +1184,9 @@
 !-- they do not necessarily need to be transferred, which is attributed to
 !-- the design of the model coupler which allocates memory for each variable.
-    IF ( rans_mode  .OR.  use_sgs_for_particles  .OR.  wang_kernel  .OR.       &
-         collision_turbulence  .OR.  nested_run )  THEN
-
-       ALLOCATE( diss_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-       IF ( rans_tke_e  .OR.  nested_run )  THEN
-          ALLOCATE( diss_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-          ALLOCATE( diss_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-       ENDIF
-
+    ALLOCATE( diss_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+
+    IF ( rans_tke_e  .OR.  nested_run )  THEN
+       ALLOCATE( diss_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+       ALLOCATE( diss_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     ENDIF
 
@@ -1195,10 +1195,7 @@
     e  => e_1;   e_p  => e_2;   te_m  => e_3
 
-    IF ( rans_mode  .OR.  use_sgs_for_particles  .OR.     &
-         wang_kernel  .OR.  collision_turbulence  .OR.  nested_run )  THEN
-       diss => diss_1
-       IF ( rans_tke_e  .OR.  nested_run )  THEN
+    diss => diss_1
+    IF ( rans_tke_e  .OR.  nested_run )  THEN
        diss_p => diss_2; tdiss_m => diss_3
-       ENDIF
     ENDIF
 
@@ -1252,4 +1249,6 @@
                 e = 0.0_wp
              ENDIF
+
+             diss = 0.0_wp
 
           ELSE
@@ -1298,5 +1297,5 @@
           ENDIF
 
-       ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 .OR. &
+       ELSEIF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 .OR. &
                 INDEX( initializing_actions, 'inifor' ) /= 0 )  THEN
 
@@ -1340,4 +1339,6 @@
              diss(nzb,:,:) = diss(nzb+1,:,:)
              diss(nzt+1,:,:) = diss(nzt,:,:)
+          ELSE
+             diss = 0.0_wp
           ENDIF
 
@@ -1383,4 +1384,6 @@
                 ENDDO
              ENDDO
+          ELSE
+             diss = 0.0_wp
           ENDIF
        ENDIF
@@ -1469,6 +1472,7 @@
 
 
+!------------------------------------------------------------------------------!
 ! Description:
-! -----------------------------------------------------------------------------!
+! ------------
 !> Pre-computation of grid-dependent and near-wall mixing length.
 !> @todo consider walls in horizontal direction at a distance further than a
@@ -1481,5 +1485,5 @@
 
     USE control_parameters,                                                    &
-        ONLY:  f, message_string, wall_adjustment_factor
+        ONLY:  f, message_string, wall_adjustment, wall_adjustment_factor
 
     USE grid_variables,                                                        &
@@ -1495,20 +1499,20 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) :: dist_dx        !< found distance devided by dx
-    INTEGER(iwp) :: i              !< index variable along x
-    INTEGER(iwp) :: ii             !< index variable along x
-    INTEGER(iwp) :: j              !< index variable along y
-    INTEGER(iwp) :: k              !< index variable along z
-    INTEGER(iwp) :: k_max_topo = 0 !< index of maximum topography height
-    INTEGER(iwp) :: kk             !< index variable along z
-    INTEGER(iwp) :: rad_i          !< search radius in grid points along x
-    INTEGER(iwp) :: rad_i_l        !< possible search radius to the left
-    INTEGER(iwp) :: rad_i_r        !< possible search radius to the right
-    INTEGER(iwp) :: rad_j          !< search radius in grid points along y
-    INTEGER(iwp) :: rad_j_n        !< possible search radius to north
-    INTEGER(iwp) :: rad_j_s        !< possible search radius to south
-    INTEGER(iwp) :: rad_k          !< search radius in grid points along z
-    INTEGER(iwp) :: rad_k_b        !< search radius in grid points along negative z
-    INTEGER(iwp) :: rad_k_t        !< search radius in grid points along positive z
+    INTEGER(iwp) :: dist_dx     !< found distance devided by dx
+    INTEGER(iwp) :: i           !< index variable along x
+    INTEGER(iwp) :: ii          !< index variable along x
+    INTEGER(iwp) :: j           !< index variable along y
+    INTEGER(iwp) :: k           !< index variable along z
+    INTEGER(iwp) :: k_max_topo  !< index of maximum topography height
+    INTEGER(iwp) :: kk          !< index variable along z
+    INTEGER(iwp) :: rad_i       !< search radius in grid points along x
+    INTEGER(iwp) :: rad_i_l     !< possible search radius to the left
+    INTEGER(iwp) :: rad_i_r     !< possible search radius to the right
+    INTEGER(iwp) :: rad_j       !< search radius in grid points along y
+    INTEGER(iwp) :: rad_j_n     !< possible search radius to north
+    INTEGER(iwp) :: rad_j_s     !< possible search radius to south
+    INTEGER(iwp) :: rad_k       !< search radius in grid points along z
+    INTEGER(iwp) :: rad_k_b     !< search radius in grid points along negative z
+    INTEGER(iwp) :: rad_k_t     !< search radius in grid points along positive z
 
     INTEGER(KIND=1), DIMENSION(:,:), ALLOCATABLE :: vic_yz !< contains a quarter of a single yz-slice of vicinity
@@ -1516,5 +1520,5 @@
     INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE :: vicinity !< contains topography information of the vicinity of (i/j/k)
 
-    REAL(wp) :: radius           !< search radius in meter
+    REAL(wp) :: radius          !< search radius in meter
 
     ALLOCATE( l_grid(1:nzt) )
@@ -1537,110 +1541,118 @@
        l_wall(nzt+1,:,:) = l_grid(nzt)
 
-       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?
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             DO  k = nzb+1, nzt
-!
-!--             Check if current gridpoint belongs to the atmosphere
-                IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
-!
-!--                Check for neighbouring grid-points.
-!--                Vertical distance, down
-                   IF ( .NOT. BTEST( wall_flags_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_0(k+1,j,i), 0 ) )              &
-                      l_wall(k,j,i) = MIN( l_grid(k), zw(k) - zu(k) )
-!
-!--                y-distance
-                   IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 )  .OR.          &
-                        .NOT. BTEST( wall_flags_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_0(k,j,i-1), 0 )  .OR.          &
-                        .NOT. BTEST( wall_flags_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_0(k-1,j-1,i), 0 )  .OR.       &
-                         .NOT. BTEST( wall_flags_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_0(k+1,j-1,i), 0 )  .OR.       &
-                         .NOT. BTEST( wall_flags_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_0(k-1,j,i-1), 0 )  .OR.       &
-                         .NOT. BTEST( wall_flags_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_0(k+1,j,i-1), 0 )  .OR.       &
-                         .NOT. BTEST( wall_flags_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_0(k,j-1,i-1), 0 )  .OR.        &
-                        .NOT. BTEST( wall_flags_0(k,j+1,i-1), 0 )  .OR.        &
-                        .NOT. BTEST( wall_flags_0(k,j-1,i+1), 0 )  .OR.        &
-                        .NOT. BTEST( wall_flags_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_0(k-1,j-1,i-1), 0 )  .OR.      &
-                        .NOT. BTEST( wall_flags_0(k-1,j+1,i-1), 0 )  .OR.      &
-                        .NOT. BTEST( wall_flags_0(k-1,j-1,i+1), 0 )  .OR.      &
-                        .NOT. BTEST( wall_flags_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_0(k+1,j-1,i-1), 0 )  .OR.      &
-                        .NOT. BTEST( wall_flags_0(k+1,j+1,i-1), 0 )  .OR.      &
-                        .NOT. BTEST( wall_flags_0(k+1,j-1,i+1), 0 )  .OR.      &
-                        .NOT. BTEST( wall_flags_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
-             ENDDO
+       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
-       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?
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt
+!
+!--                Check if current gridpoint belongs to the atmosphere
+                   IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
+!
+!--                   Check for neighbouring grid-points.
+!--                   Vertical distance, down
+                      IF ( .NOT. BTEST( wall_flags_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_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_0(k,j-1,i), 0 )  .OR.         &
+                           .NOT. BTEST( wall_flags_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_0(k,j,i-1), 0 )  .OR.         &
+                           .NOT. BTEST( wall_flags_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_0(k-1,j-1,i), 0 )  .OR.       &
+                           .NOT. BTEST( wall_flags_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_0(k+1,j-1,i), 0 )  .OR.       &
+                           .NOT. BTEST( wall_flags_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_0(k-1,j,i-1), 0 )  .OR.       &
+                           .NOT. BTEST( wall_flags_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_0(k+1,j,i-1), 0 )  .OR.       &
+                           .NOT. BTEST( wall_flags_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_0(k,j-1,i-1), 0 )  .OR.        &
+                           .NOT. BTEST( wall_flags_0(k,j+1,i-1), 0 )  .OR.        &
+                           .NOT. BTEST( wall_flags_0(k,j-1,i+1), 0 )  .OR.        &
+                           .NOT. BTEST( wall_flags_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_0(k-1,j-1,i-1), 0 )  .OR.      &
+                           .NOT. BTEST( wall_flags_0(k-1,j+1,i-1), 0 )  .OR.      &
+                           .NOT. BTEST( wall_flags_0(k-1,j-1,i+1), 0 )  .OR.      &
+                           .NOT. BTEST( wall_flags_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_0(k+1,j-1,i-1), 0 )  .OR.      &
+                           .NOT. BTEST( wall_flags_0(k+1,j+1,i-1), 0 )  .OR.      &
+                           .NOT. BTEST( wall_flags_0(k+1,j-1,i+1), 0 )  .OR.      &
+                           .NOT. BTEST( wall_flags_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) )
+
+                ENDDO  !k loop
+             ENDDO  !j loop
+          ENDDO  !i loop
+
+       ENDIF  !if wall_adjustment
 
     ELSE
 !
-!-- Initialize the mixing length in case of a RANS simulation
+!--    Initialize the mixing length in case of a RANS simulation
        ALLOCATE( l_black(nzb:nzt+1) )
 
@@ -1662,4 +1674,5 @@
 !
 !--    Get height level of highest topography within local subdomain
+       k_max_topo = 0
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
@@ -4043,15 +4056,10 @@
  SUBROUTINE diffusion_e( var, var_reference )
 
-#if defined( _OPENACC )
     USE arrays_3d,                                                             &
-        ONLY:  ddzu, dd2zu, ddzw, drho_air, rho_air_zw
+        ONLY:  dd2zu, ddzu, ddzw, drho_air, rho_air_zw
 
     USE control_parameters,                                                    &
-        ONLY:  atmos_ocean_sign, use_single_reference_value,                   &
-               wall_adjustment, wall_adjustment_factor
-#else
-    USE arrays_3d,                                                             &
-        ONLY:  ddzu, ddzw, drho_air, rho_air_zw
-#endif
+        ONLY:  atmos_ocean_sign, use_single_reference_value
+
     USE grid_variables,                                                        &
         ONLY:  ddx2, ddy2
@@ -4070,150 +4078,145 @@
     INTEGER(iwp) ::  m              !< running index surface elements
 
-    REAL(wp)     ::  flag           !< flag to mask topography
+    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)     ::  ll             !< adjusted l
     REAL(wp)     ::  var_reference  !< reference temperature
-#if defined( _OPENACC )
-!
-!-- From mixing_length_les:
-    REAL(wp)     :: l_stable        !< mixing length according to stratification
-!
-!-- From mixing_length_rans:
-    REAL(wp)     :: duv2_dz2        !< squared vertical gradient of wind vector
-    REAL(wp)     :: l_diss          !< mixing length for dissipation
-    REAL(wp)     :: rif             !< Richardson flux number
-!
-!-- From both:
-    REAL(wp)     :: dvar_dz         !< vertical gradient of var
-#endif
-
-    REAL(wp)     ::  dissipation  !< TKE dissipation
-
-    REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< temperature
-
-
-!
-!-- Calculate the tendency terms
+
+    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:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  var  !< temperature
+
+!
+!-- Calculate the dissipation
+    IF ( les_dynamic .OR. les_mw )  THEN
+
+       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
+       !$ACC PRIVATE(l, l_stable, dvar_dz) &
+       !$ACC PRESENT(diss, e, var, wall_flags_0) &
+       !$ACC PRESENT(dd2zu, l_grid, l_wall)
+       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
+                   IF ( use_single_reference_value )  THEN
+                     l_stable(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) )               &
+                                 / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
+                   ENDIF
+                ELSE
+                   l_stable(k) = l_grid(k)
+                ENDIF
+
+             ENDDO
+
+             !$ACC LOOP PRIVATE(k)
+             !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                        &
+                              * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+
+             ENDDO
+
+          ENDDO
+       ENDDO
+
+    ELSEIF ( rans_tke_l )  THEN
+
+      !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
+      !$ACC PRIVATE(l_stable, duv2_dz2, rif, dvar_dz) &
+      !$ACC PRESENT(diss, e, u, v, var, wall_flags_0) &
+      !$ACC PRESENT(dd2zu, l_black, l_wall)
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+!
+!--          Calculate Richardson-flux number
+             IF ( use_single_reference_value )  THEN
+                !$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)
+
+                   duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2  &
+                              + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2  &
+                              + 1E-30_wp
+
+                   rif(k) = MIN( MAX( g / var_reference * dvar_dz / duv2_dz2, -5.0_wp ),  1.0_wp )
+                ENDDO
+             ELSE
+                !$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)
+
+                   duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2  &
+                              + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2  &
+                              + 1E-30_wp
+
+                   rif(k) = MIN( MAX( g / var(k,j,i) * dvar_dz / duv2_dz2, -5.0_wp ),  1.0_wp )
+                ENDDO
+             ENDIF
+!
+!--          Calculate diabatic mixing length using Dyer-profile functions
+             !$ACC LOOP PRIVATE(k)
+             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) )
+                ELSE
+                   l_stable(k) = l_wall(k,j,i) * SQRT( 1.0_wp - 16.0_wp * rif(k) )
+                ENDIF
+             ENDDO
+
+             !$ACC LOOP PRIVATE(k)
+             !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)     &
+                            * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+
+             ENDDO
+
+          ENDDO
+       ENDDO
+
+!-- Note, in case of rans_tke_e, the dissipation is already calculated
+!-- in prognostic_equations
+    ENDIF
+
     !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
-    !$ACC PRIVATE(flag, l, ll, dissipation) &
-    !$ACC PRIVATE(l_stable, duv2_dz2, l_diss, rif, dvar_dz) &
-    !$ACC PRESENT(e, u, v, km, var, wall_flags_0) &
-    !$ACC PRESENT(ddzu, dd2zu, ddzw, rho_air_zw, drho_air) &
-    !$ACC PRESENT(l_black, l_grid, l_wall) &
-    !$ACC PRESENT(diss, tend)
+    !$ACC PRESENT(diss, e, km, tend, wall_flags_0) &
+    !$ACC PRESENT(ddzu, ddzw, rho_air_zw, drho_air)
     DO  i = nxl, nxr
        DO  j = nys, nyn
           DO  k = nzb+1, nzt
-!
-!--          Predetermine flag to mask topography
-             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-
-!
-!--          Calculate dissipation
-             IF ( les_dynamic .OR. les_mw )  THEN
-
-#if defined( _OPENACC )
-!
-!--             Cannot call subroutine mixing_length_les because after adding all required
-!--             OpenACC directives the execution crashed reliably when inside the called
-!--             subroutine. I'm not sure why that is...
-                dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
-                IF ( dvar_dz > 0.0_wp ) THEN
-                   IF ( use_single_reference_value )  THEN
-                      l_stable = 0.76_wp * SQRT( e(k,j,i) )                                &
-                                         / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
-                   ELSE
-                      l_stable = 0.76_wp * SQRT( e(k,j,i) )                                &
-                                         / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
-                   ENDIF
-                ELSE
-                   l_stable = l_grid(k)
-                ENDIF
-!
-!--             Adjustment of the mixing length
-                IF ( wall_adjustment )  THEN
-                   l  = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k), l_stable )
-                   ll = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) )
-                ELSE
-                   l  = MIN( l_grid(k), l_stable )
-                   ll = l_grid(k)
-                ENDIF
-#else
-                CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
-#endif
-
-                dissipation = ( 0.19_wp + 0.74_wp * l / ll )                   &
-                              * e(k,j,i) * SQRT( e(k,j,i) ) / l
-
-             ELSEIF ( rans_tke_l )  THEN
-
-#if defined( _OPENACC )
-!
-!--             Same reason as above...
-                dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
-
-                duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2                 &
-                           + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2                 &
-                           + 1E-30_wp
-
-                IF ( use_single_reference_value )  THEN
-                   rif = g / var_reference * dvar_dz / duv2_dz2
-                ELSE
-                   rif = g / var(k,j,i) * dvar_dz / duv2_dz2
-                ENDIF
-
-                rif = MAX( rif, -5.0_wp )
-                rif = MIN( rif,  1.0_wp )
-
-!
-!--             Calculate diabatic mixing length using Dyer-profile functions
-                IF ( rif >= 0.0_wp )  THEN
-                   l      = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif ), l_wall(k,j,i) )
-                   l_diss = l
-                ELSE
-!
-!--                In case of unstable stratification, use mixing length of neutral case
-!--                for l, but consider profile functions for l_diss
-                   l      = l_wall(k,j,i)
-                   l_diss = l * SQRT( 1.0_wp - 16.0_wp * rif )
-                ENDIF
-#else
-                CALL mixing_length_rans( i, j, k, l, ll, var, var_reference )
-#endif
-
-                dissipation = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll
-
-                diss(k,j,i) = dissipation * flag
-
-             ELSEIF ( rans_tke_e )  THEN
-
-                dissipation = diss(k,j,i)
-
-             ENDIF
-
-             tend(k,j,i) = tend(k,j,i) + (                                     &
-                                           (                                   &
-                       ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )     &
-                     - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )     &
-                                           ) * ddx2  * flag                    &
-                                         + (                                   &
-                       ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )     &
-                     - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )     &
-                                           ) * ddy2  * flag                    &
-                                         + (                                   &
-            ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1)    &
-                                                          * rho_air_zw(k)      &
-          - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)      &
-                                                          * rho_air_zw(k-1)    &
-                                           ) * ddzw(k) * drho_air(k)           &
-                                         ) * flag * dsig_e                     &
-                          - dissipation * flag
-!
-!--          Store dissipation if needed for calculating the sgs particle
-!--          velocities
-             IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles  .OR.         &
-                  wang_kernel  .OR.  collision_turbulence  ) )  THEN
-                diss(k,j,i) = dissipation * flag
-             ENDIF
+
+             tend(k,j,i) = tend(k,j,i) + (                                            &
+                                           (                                          &
+                       ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )            &
+                     - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )            &
+                                           ) * ddx2                                   &
+                                         + (                                          &
+                       ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )            &
+                     - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )            &
+                                           ) * ddy2                                   &
+                                         + (                                          &
+            ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1)           &
+                                                          * rho_air_zw(k)             &
+          - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)             &
+                                                          * rho_air_zw(k-1)           &
+                                           ) * ddzw(k) * drho_air(k)                  &
+                                         ) * dsig_e                                   &
+                                           * MERGE( 1.0_wp, 0.0_wp,                   &
+                                                    BTEST( wall_flags_0(k,j,i), 0 ) ) &
+                          - diss(k,j,i)
 
           ENDDO
@@ -4258,5 +4261,8 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  ddzu, ddzw, drho_air, rho_air_zw
+        ONLY:  dd2zu, ddzu, ddzw, drho_air, rho_air_zw
+
+    USE control_parameters,                                                    &
+        ONLY:  atmos_ocean_sign, use_single_reference_value
 
     USE grid_variables,                                                        &
@@ -4278,48 +4284,96 @@
     INTEGER(iwp) ::  surf_s         !< Start index of surface elements at (j,i)-gridpoint
 
-    REAL(wp)     ::  flag           !< flag to mask topography
+    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)     ::  ll             !< adjusted l
     REAL(wp)     ::  var_reference  !< reference temperature
 
-    REAL(wp), DIMENSION(nzb+1:nzt) ::  dissipation  !< dissipation of TKE
-
-    REAL(wp), DIMENSION(:,:,:), POINTER ::  var     !< temperature
-
-!
-!-- Calculate the mixing length (for dissipation)
+    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:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  var  !< temperature
+
+!
+!-- Calculate the mixing length and dissipation...
+!-- ...in case of LES
+    IF ( les_dynamic .OR. les_mw )  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
+             IF ( use_single_reference_value )  THEN
+                l_stable(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) )               &
+                            / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
+             ENDIF
+          ELSE
+             l_stable(k) = l_grid(k)
+          ENDIF
+       ENDDO
+
+       !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                      &
+                        * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+       ENDDO
+
+!
+!-- ...in case of RANS
+    ELSEIF ( rans_tke_l )  THEN
+
+!
+!--    Calculate Richardson-flux number
+       IF ( use_single_reference_value )  THEN
+          DO  k = nzb+1, nzt
+             dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+
+             duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2           &
+                        + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2           &
+                        + 1E-30_wp
+
+             rif(k) = MIN( MAX( g / var_reference * dvar_dz / duv2_dz2, -5.0_wp ),  1.0_wp )
+          ENDDO
+       ELSE
+          DO  k = nzb+1, nzt
+             dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+
+             duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2           &
+                        + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2           &
+                        + 1E-30_wp
+
+             rif(k) = MIN( MAX( g / var(k,j,i) * dvar_dz / duv2_dz2, -5.0_wp ), 1.0_wp )
+          ENDDO
+       ENDIF
+!
+!--    Calculate diabatic mixing length using Dyer-profile functions
+       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) )
+          ELSE
+             l_stable(k) = l_wall(k,j,i) * SQRT( 1.0_wp - 16.0_wp * rif(k) )
+          ENDIF
+
+       ENDDO
+
+       !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)     &
+                      * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+       ENDDO
+
+!-- Note, in case of rans_tke_e, the dissipation is already calculated
+!-- in prognostic_equations
+    ENDIF
+
+!
+!-- Calculate the tendency term
+    !DIR$ IVDEP
     DO  k = nzb+1, nzt
-!
-!--    Predetermine flag to mask topography
-       flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-
-!
-!--    Calculate dissipation...
-!--    ...in case of LES
-       IF ( les_dynamic .OR. les_mw )  THEN
-
-          CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
-
-          dissipation(k) = ( 0.19_wp + 0.74_wp * l / ll )                      &
-                           * e(k,j,i) * SQRT( e(k,j,i) ) / l
-
-!
-!--    ...in case of RANS
-       ELSEIF ( rans_tke_l )  THEN
-
-          CALL mixing_length_rans( i, j, k, l, ll, var, var_reference  )
-
-          dissipation(k) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll
-
-          diss(k,j,i) = dissipation(k) * flag
-
-       ELSEIF ( rans_tke_e )  THEN
-
-          dissipation(k) = diss(k,j,i)
-
-       ENDIF
-
-!
-!--    Calculate the tendency term
+
        tend(k,j,i) = tend(k,j,i) + (                                           &
                                       (                                        &
@@ -4337,23 +4391,21 @@
                                                          * rho_air_zw(k-1)     &
                                       ) * ddzw(k) * drho_air(k)                &
-                                   ) * flag * dsig_e                           &
-                                 - dissipation(k) * flag
+                                   ) * dsig_e                                  &
+                                     * MERGE( 1.0_wp, 0.0_wp,                  &
+                                              BTEST( wall_flags_0(k,j,i), 0 ) )&
+                                 - diss(k,j,i)
 
     ENDDO
 
 !
-!-- Store dissipation if needed for calculating the sgs particle velocities
+!-- Set boundary conditions of dissipation if needed for calculating the sgs
+!-- particle velocities.
+!-- Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:)
+!-- For each surface type determine start and end index (in case of elevated
+!-- topography several up/downward facing surfaces may exist.
+!-- Note, bc cannot be set in tcm_boundary conditions as the dissipation 
+!-- in LES mode is only a diagnostic quantity. 
     IF ( .NOT. rans_tke_e .AND.  ( use_sgs_for_particles  .OR.  wang_kernel    &
           .OR.  collision_turbulence ) )  THEN
-       DO  k = nzb+1, nzt
-          diss(k,j,i) = dissipation(k) * MERGE( 1.0_wp, 0.0_wp,                &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-       ENDDO
-!
-!--    Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:)
-!--    For each surface type determine start and end index (in case of elevated
-!--    topography several up/downward facing surfaces may exist.
-!--    Note, bc cannot be set in tcm_boundary conditions as the dissipation 
-!--    in LES mode is only a diagnostic quantity. 
        surf_s = bc_h(0)%start_index(j,i)
        surf_e = bc_h(0)%end_index(j,i)
@@ -4490,118 +4542,4 @@
 ! Description:
 ! ------------
-!> Calculate mixing length for LES mode.
-!------------------------------------------------------------------------------!
- SUBROUTINE mixing_length_les( i, j, k, l, ll, var, var_reference )
-
-    USE arrays_3d,                                                             &
-        ONLY:  dd2zu
-
-    USE control_parameters,                                                    &
-        ONLY:  atmos_ocean_sign, use_single_reference_value,                   &
-               wall_adjustment, wall_adjustment_factor
-
-    IMPLICIT NONE
-
-    INTEGER(iwp) :: i   !< loop index
-    INTEGER(iwp) :: j   !< loop index
-    INTEGER(iwp) :: k   !< loop index
-
-    REAL(wp)     :: dvar_dz         !< vertical gradient of var
-    REAL(wp)     :: l               !< mixing length
-    REAL(wp)     :: l_stable        !< mixing length according to stratification
-    REAL(wp)     :: ll              !< adjusted l_grid
-    REAL(wp)     :: var_reference   !< var at reference height
-
-    REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< temperature
-
-
-    dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
-    IF ( dvar_dz > 0.0_wp ) THEN
-       IF ( use_single_reference_value )  THEN
-          l_stable = 0.76_wp * SQRT( e(k,j,i) )                                &
-                             / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
-       ELSE
-          l_stable = 0.76_wp * SQRT( e(k,j,i) )                                &
-                             / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
-       ENDIF
-    ELSE
-       l_stable = l_grid(k)
-    ENDIF
-!
-!-- Adjustment of the mixing length
-    IF ( wall_adjustment )  THEN
-       l  = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k), l_stable )
-       ll = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) )
-    ELSE
-       l  = MIN( l_grid(k), l_stable )
-       ll = l_grid(k)
-    ENDIF
-
- END SUBROUTINE mixing_length_les
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Calculate mixing length for RANS mode.
-!------------------------------------------------------------------------------!
- SUBROUTINE mixing_length_rans( i, j, k, l, l_diss, var, var_reference  )
-
-    USE arrays_3d,                                                             &
-        ONLY:  dd2zu
-
-    USE control_parameters,                                                    &
-        ONLY:  atmos_ocean_sign, use_single_reference_value
-
-    IMPLICIT NONE
-
-    INTEGER(iwp) :: i   !< loop index
-    INTEGER(iwp) :: j   !< loop index
-    INTEGER(iwp) :: k   !< loop index
-
-    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)     :: l_diss          !< mixing length for dissipation
-    REAL(wp)     :: rif             !< Richardson flux number
-    REAL(wp)     :: var_reference   !< var at reference height
-
-    REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< temperature
-
-
-    dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
-
-    duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2                 &
-               + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2                 &
-               + 1E-30_wp
-
-    IF ( use_single_reference_value )  THEN
-       rif = g / var_reference * dvar_dz / duv2_dz2
-    ELSE
-       rif = g / var(k,j,i) * dvar_dz / duv2_dz2
-    ENDIF
-
-    rif = MAX( rif, -5.0_wp )
-    rif = MIN( rif,  1.0_wp )
-
-!
-!-- Calculate diabatic mixing length using Dyer-profile functions
-    IF ( rif >= 0.0_wp )  THEN
-       l      = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif ), l_wall(k,j,i) )
-       l_diss = l
-    ELSE
-!
-!--    In case of unstable stratification, use mixing length of neutral case
-!--    for l, but consider profile functions for l_diss
-       l      = l_wall(k,j,i)
-       l_diss = l * SQRT( 1.0_wp - 16.0_wp * rif )
-    ENDIF
-
- END SUBROUTINE mixing_length_rans
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
 !> Computation of the turbulent diffusion coefficients for momentum and heat.
 !> @bug unstable stratification is not properly considered for kh in rans mode.
@@ -4624,5 +4562,5 @@
     REAL(wp) ::  var_reference  !< reference temperature
 
-    REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< temperature
+    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  var  !< temperature
 
 
@@ -4797,12 +4735,9 @@
  SUBROUTINE tcm_diffusivities_default( var, var_reference )
 
-#if defined( _OPENACC )
     USE arrays_3d,                                                             &
         ONLY:  dd2zu
 
     USE control_parameters,                                                    &
-        ONLY:  atmos_ocean_sign, use_single_reference_value,                   &
-               wall_adjustment, wall_adjustment_factor
-#endif
+        ONLY:  atmos_ocean_sign, use_single_reference_value
 
     USE statistics,                                                            &
@@ -4818,17 +4753,15 @@
     INTEGER(iwp) ::  tn                  !< thread number
 
-    REAL(wp)     ::  flag                !< topography flag
-    REAL(wp)     ::  l                   !< mixing length
-    REAL(wp)     ::  ll                  !< adjusted mixing length
+    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)     ::  var_reference       !< reference temperature
-#if defined( _OPENACC )
-!
-!-- From mixing_length_les:
-    REAL(wp)     :: l_stable        !< mixing length according to stratification
-    REAL(wp)     :: dvar_dz         !< vertical gradient of var
-#endif
-
-    REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< 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
+
+    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  var  !< temperature
 
 !
@@ -4844,63 +4777,60 @@
 !
 !-- Compute the turbulent diffusion coefficient for momentum
-    !$OMP PARALLEL PRIVATE (i,j,k,l,ll,sr,tn,flag)
+    !$OMP PARALLEL PRIVATE (i,j,k,l,sr,tn)
 !$  tn = omp_get_thread_num()
 
     IF ( les_dynamic .OR. les_mw )  THEN
        !$OMP DO
-       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(sr) &
-       !$ACC PRIVATE(flag, dvar_dz, l_stable, l, ll) &
+       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
+       !$ACC PRIVATE(dvar_dz, l, l_stable, l_v) &
        !$ACC PRESENT(wall_flags_0, var, dd2zu, e, l_wall, l_grid, rmask) &
        !$ACC PRESENT(kh, km, sums_l_l)
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
+             !$ACC LOOP PRIVATE(k)
              DO  k = nzb+1, nzt
-
-                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-
-!
-!--             Determine the mixing length for LES closure
-#if defined( _OPENACC )
-!
-!--             Cannot call subroutine mixing_length_les because after adding all required
-!--             OpenACC directives the execution crashed reliably when inside the called
-!--             subroutine. I'm not sure why that is...
+!
+!--             Determine the mixing length
+!--             @note The following code cannot be transferred to a subroutine
+!--             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).
                 dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
                 IF ( dvar_dz > 0.0_wp ) THEN
                    IF ( use_single_reference_value )  THEN
-                      l_stable = 0.76_wp * SQRT( e(k,j,i) )                                &
-                                         / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
+                      l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )                  &
+                                  / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
                    ELSE
-                      l_stable = 0.76_wp * SQRT( e(k,j,i) )                                &
-                                         / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
+                      l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )               &
+                                  / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
                    ENDIF
                 ELSE
-                   l_stable = l_grid(k)
+                   l_stable(k) = l_grid(k)
                 ENDIF
-!
-!--             Adjustment of the mixing length
-                IF ( wall_adjustment )  THEN
-                   l  = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k), l_stable )
-                   ll = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) )
-                ELSE
-                   l  = MIN( l_grid(k), l_stable )
-                   ll = l_grid(k)
-                ENDIF
-#else
-                CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
-#endif
-
+
+             ENDDO
+
+             !$ACC LOOP PRIVATE(k)
+             !DIR$ IVDEP
+             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_0(k,j,i), 0 ) )
+                l = l_v(k)
 !
 !--             Compute diffusion coefficients for momentum and heat
-                km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag
-                kh(k,j,i) = ( 1.0_wp + 2.0_wp * l / ll ) * km(k,j,i) * flag
-!
-!--             Summation for averaged profile (cf. flow_statistics)
-                DO  sr = 0, statistic_regions
-                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr)   &
-                                                             * flag
+                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)
+
+             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) + l_v(k) * rmask(j,i,sr)
                 ENDDO
-
-             ENDDO
+             ENDDO
+
           ENDDO
        ENDDO
@@ -4909,24 +4839,67 @@
 
        !$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_0, var, dd2zu, e, u, v, l_wall, l_black, rmask) &
+       !$ACC PRESENT(kh, km, sums_l_l)
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
+!
+!--          Calculate Richardson-flux number
+             IF ( use_single_reference_value )  THEN
+                !$ACC LOOP PRIVATE(k)
+                !DIR$ IVDEP
+                DO  k = nzb+1, nzt
+                   dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+
+                   duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2  &
+                              + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2  &
+                              + 1E-30_wp
+
+                   rif(k) = MIN( MAX( g / var_reference * dvar_dz / duv2_dz2, -5.0_wp ),  1.0_wp )
+                ENDDO
+             ELSE
+                !$ACC LOOP PRIVATE(k)
+                !DIR$ IVDEP
+                DO  k = nzb+1, nzt
+                   dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+
+                   duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2  &
+                              + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2  &
+                              + 1E-30_wp
+
+                   rif(k) = MIN( MAX( g / var(k,j,i) * dvar_dz / duv2_dz2, -5.0_wp ),  1.0_wp )
+                ENDDO
+             ENDIF
+!
+!--          Calculate diabatic mixing length using Dyer-profile functions
+!--          In case of unstable stratification, use mixing length of neutral case
+             !$ACC LOOP PRIVATE(k)
              DO  k = nzb+1, nzt
-
-                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-!
-!--             Mixing length for RANS mode with TKE-l closure
-                CALL mixing_length_rans( i, j, k, l, ll, var, var_reference )
-!
-!--             Compute diffusion coefficients for momentum and heat
-                km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag
-                kh(k,j,i) = km(k,j,i) / prandtl_number * flag
-!
-!--             Summation for averaged profile (cf. flow_statistics)
-                DO  sr = 0, statistic_regions
-                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr)   &
-                                                             * flag
+                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) )
+                ELSE
+                   l_stable(k)  = l_wall(k,j,i)
+                ENDIF
+
+             ENDDO
+!
+!--          Compute diffusion coefficients for momentum and heat
+             !$ACC LOOP PRIVATE(k)
+             !DIR$ IVDEP
+             DO  k = nzb+1, nzt
+                l_v(k)    = l_stable(k) * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+                km(k,j,i) = c_0 * l_v(k) * SQRT( e(k,j,i) )
+                kh(k,j,i) = km(k,j,i) / prandtl_number
+             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) + l_v(k) * rmask(j,i,sr)
                 ENDDO
-
-             ENDDO
+             ENDDO
+
           ENDDO
        ENDDO
@@ -4935,22 +4908,32 @@
 
        !$OMP DO
+       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
+       !$ACC PRIVATE(l_v) &
+       !$ACC PRESENT(wall_flags_0, e, diss, rmask) &
+       !$ACC PRESENT(kh, km, sums_l_l)
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
+!
+!--          Compute diffusion coefficients for momentum and heat
+             !$ACC LOOP PRIVATE(k)
+             !DIR$ IVDEP
              DO  k = nzb+1, nzt
 
-                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-!
-!--             Compute diffusion coefficients for momentum and heat
-                km(k,j,i) = c_0**4 * e(k,j,i)**2 / ( diss(k,j,i) + 1.0E-30_wp ) * flag
-                kh(k,j,i) = km(k,j,i) / prandtl_number * flag
-!
-!--             Summation for averaged profile of mixing length (cf. flow_statistics)
-                DO  sr = 0, statistic_regions
-                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) +                     &
-                      c_0**3 * e(k,j,i) * SQRT(e(k,j,i)) /                     &
-                      ( diss(k,j,i) + 1.0E-30_wp ) * rmask(j,i,sr) * flag
+                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_0(k,j,i), 0 ) )
+
+                km(k,j,i) = c_0 * SQRT( e(k,j,i) ) * l_v(k)
+                kh(k,j,i) = km(k,j,i) / prandtl_number
+
+             ENDDO
+!
+!--          Summation for averaged profile of mixing length (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) + l_v(k) * rmask(j,i,sr)
                 ENDDO
-
-             ENDDO
+             ENDDO
+
           ENDDO
        ENDDO
