Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 2912)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 2918)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Add check for 1D model
+! 
+! 2883 2018-03-14 08:29:10Z Giersch
 ! dt_dopr_listing is not set to the default value zero anymore
 ! 
@@ -3993,4 +3996,10 @@
        CALL message( 'check_parameters', 'PA0138', 1, 2, 0, 6, 0 )
     ENDIF
+    IF ( TRIM( mixing_length_1d ) /= 'as_in_3d_model'  .AND.                   &
+         TRIM( dissipation_1d ) == 'as_in_3d_model' )  THEN
+       message_string = 'dissipation_1d = "' // TRIM( dissipation_1d ) //      &
+                        '" requires mixing_length_1d = "as_in_3d_model"'
+       CALL message( 'check_parameters', 'PA0485', 1, 2, 0, 6, 0 )
+    ENDIF
 
 !
Index: palm/trunk/SOURCE/init_grid.f90
===================================================================
--- palm/trunk/SOURCE/init_grid.f90	(revision 2912)
+++ palm/trunk/SOURCE/init_grid.f90	(revision 2918)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Moved init_mixing_length to turbulence_closure_mod.f90
+! 
+! 2897 2018-03-15 11:47:16Z suehring
 ! Relax restrictions for topography input, terrain and building heights can be
 ! input separately and are not mandatory any more. 
@@ -284,5 +287,4 @@
 ! -----------------------------------------------------------------------------!
 !> Creating grid depending constants
-!> @todo: Move initialization mixing length
 !> @todo: Rearrange topo flag list
 !> @todo: reference 3D buildings on top of orography is not tested and may need
@@ -517,8 +519,4 @@
 
 !
-!-- Calculated grid-dependent as well as near-wall mixing length 
-    CALL init_mixing_length
-
-!
 !-- Determine the maximum level of topography. It is used for
 !-- steering the degradation of order of the applied advection scheme, 
@@ -1342,154 +1340,4 @@
  END SUBROUTINE filter_topography
 
-! Description:
-! -----------------------------------------------------------------------------!
-!> Pre-computation of grid-dependent and near-wall mixing length.
-!> @todo: move subroutine to turbulence module developed by M1 (Tobi).
-!------------------------------------------------------------------------------!
- SUBROUTINE init_mixing_length
-
-    USE arrays_3d,                                                             &
-        ONLY:  dzw, l_grid, l_wall, zu, zw
-
-    USE control_parameters,                                                    &
-        ONLY:  message_string, wall_adjustment_factor
-
-    USE grid_variables,                                                        &
-        ONLY:  dx, dy
-
-    USE indices,                                                               &
-        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,     &
-               wall_flags_0
-    
-    USE kinds
-
-    IMPLICIT NONE
-
-    INTEGER(iwp) ::  i             !< index variable along x
-    INTEGER(iwp) ::  j             !< index variable along y
-    INTEGER(iwp) ::  k             !< index variable along z
-
-    ALLOCATE( l_grid(1:nzt) )
-    ALLOCATE( l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-!
-!-- Compute the grid-dependent mixing length.
-    DO  k = 1, nzt
-       l_grid(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)
-    ENDDO
-    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
-       ENDDO
-    ENDDO
-!
-!-- Set lateral boundary conditions for l_wall
-    CALL exchange_horiz( l_wall, nbgp )     
-
- END SUBROUTINE init_mixing_length
 
 ! Description:
Index: palm/trunk/SOURCE/model_1d_mod.f90
===================================================================
--- palm/trunk/SOURCE/model_1d_mod.f90	(revision 2912)
+++ palm/trunk/SOURCE/model_1d_mod.f90	(revision 2918)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! - rename l_black into l1d_init
+! - calculate l_grid within init_1d_model and save it as l1d_init
+! 
+! 2718 2018-01-02 08:49:38Z maronga
 ! Corrected "Former revisions" section
 ! 
@@ -128,9 +132,12 @@
 !> @todo harmonize code with new surface_layer_fluxes module
 !> @bug 1D model crashes when using small grid spacings in the order of 1 m
+!> @fixme option "as_in_3d_model" seems to eb an inappropriate option because
+!>        the 1D model uses different turbulence closure approaches at least if
+!>        the 3D model is set to LES-mode.
 !------------------------------------------------------------------------------!
  MODULE model_1d_mod
 
     USE arrays_3d,                                                             &
-        ONLY:  dd2zu, ddzu, ddzw, dzu, l_grid, pt_init, q_init, ug, u_init,    &
+        ONLY:  dd2zu, ddzu, ddzw, dzu, dzw, pt_init, q_init, ug, u_init,       &
                vg, v_init, zu
     
@@ -190,6 +197,6 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  kh1d     !< turbulent diffusion coefficient for heat (1d-model)
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  km1d     !< turbulent diffusion coefficient for momentum (1d-model)
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l_black  !< mixing length Blackadar (1d-model)
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  l1d      !< mixing length for turbulent diffusion coefficients (1d-model)
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l1d_init !< initial mixing length (1d-model)
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  l1d_diss !< mixing length for dissipation (1d-model)
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  rif1d    !< Richardson flux number (1d-model)
@@ -255,4 +262,7 @@
  SUBROUTINE init_1d_model
  
+    USE grid_variables,                                                        &
+        ONLY:  dx, dy
+
     IMPLICIT NONE
 
@@ -267,5 +277,5 @@
     ALLOCATE( diss1d(nzb:nzt+1), diss1d_p(nzb:nzt+1),                          &
               e1d(nzb:nzt+1), e1d_p(nzb:nzt+1), kh1d(nzb:nzt+1),               &
-              km1d(nzb:nzt+1), l_black(nzb:nzt+1), l1d(nzb:nzt+1),             &
+              km1d(nzb:nzt+1), l1d(nzb:nzt+1), l1d_init(nzb:nzt+1),            &
               l1d_diss(nzb:nzt+1), rif1d(nzb:nzt+1), te_diss(nzb:nzt+1),       &
               te_dissm(nzb:nzt+1), te_e(nzb:nzt+1),                            &
@@ -286,5 +296,5 @@
 !
 !--    Compute the mixing length
-       l_black(nzb) = 0.0_wp
+       l1d_init(nzb) = 0.0_wp
 
        IF ( TRIM( mixing_length_1d ) == 'blackadar' )  THEN
@@ -299,17 +309,26 @@
 
           DO  k = nzb+1, nzt+1
-             l_black(k) = kappa * zu(k) / ( 1.0_wp + kappa * zu(k) / lambda )
+             l1d_init(k) = kappa * zu(k) / ( 1.0_wp + kappa * zu(k) / lambda )
           ENDDO
 
        ELSEIF ( TRIM( mixing_length_1d ) == 'as_in_3d_model' )  THEN
 !
-!--       Use the same mixing length as in 3D model
-          l_black(1:nzt) = l_grid
-          l_black(nzt+1) = l_black(nzt)
+!--       Use the same mixing length as in 3D model (LES-mode)
+          !### TODO: rename (delete?) this option
+          !###  As the mixing length is different between RANS and LES mode, it
+          !###  must be distinguished here between these modes. For RANS mode,
+          !###  the mixing length is calculated accoding to Blackadar, which is
+          !###  the other option at this point.
+          !###  Maybe delete this option entirely (not appropriate in LES case)
+          !###  2018-03-20, gronemeier
+          DO  k = nzb+1, nzt
+             l1d_init(k)  = ( dx * dy * dzw(k) )**0.33333333333333_wp
+          ENDDO
+          l1d_init(nzt+1) = l1d_init(nzt)
 
        ENDIF
     ENDIF
-    l1d      = l_black
-    l1d_diss = l_black
+    l1d      = l1d_init
+    l1d_diss = l1d_init
     u1d      = u_init
     u1d_p    = u_init
@@ -457,5 +476,5 @@
                    diss1d(k) = c_m**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
                 ELSEIF ( dissipation_1d == 'as_in_3d_model' )  THEN
-                   diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) /  l_grid(k)  &
+                   diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) &
                                ) * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
                 ENDIF
@@ -520,5 +539,5 @@
                 diss1d(k) = c_m**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
              ELSEIF ( dissipation_1d == 'as_in_3d_model' )  THEN
-                diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l_grid(k) )    &
+                diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) )  &
                             * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
              ENDIF
@@ -804,9 +823,9 @@
                 DO  k = nzb+1, nzt
                    IF ( rif1d(k) >= 0.0_wp )  THEN
-                      l1d(k) = l_black(k) / ( 1.0_wp + 5.0_wp * rif1d(k) )
+                      l1d(k) = l1d_init(k) / ( 1.0_wp + 5.0_wp * rif1d(k) )
                       l1d_diss(k) = l1d(k)
                    ELSE
-                      l1d(k) = l_black(k)
-                      l1d_diss(k) = l_black(k) *                               &
+                      l1d(k) = l1d_init(k)
+                      l1d_diss(k) = l1d_init(k) *                              &
                                     SQRT( 1.0_wp - 16.0_wp * rif1d(k) )
                    ENDIF
@@ -819,7 +838,7 @@
                                      SQRT( g / pt_init(k) * dpt_dz ) + 1E-5_wp
                    ELSE
-                      l_stable = l_grid(k)
+                      l_stable = l1d_init(k)
                    ENDIF
-                   l1d(k) = MIN( l_grid(k), l_stable )
+                   l1d(k) = MIN( l1d_init(k), l_stable )
                    l1d_diss(k) = l1d(k)
                 ENDDO
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 2912)
+++ palm/trunk/SOURCE/modules.f90	(revision 2918)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! -l_grid, -l_wall
+! 
+! 2906 2018-03-19 08:56:40Z Giersch
 ! Module control_parameters has been extended with ENVIRONMENT variables 
 ! read/write_svf 
@@ -621,5 +624,4 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  hyp                    !< hydrostatic pressure 
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  inflow_damping_factor  !< used for turbulent inflow (non-cyclic boundary conditions) 
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l_grid                 !< geometric mean of grid sizes dx, dy, dz 
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  ptdf_x                 !< damping factor for potential temperature in x-direction 
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  ptdf_y                 !< damping factor for potential temperature in y-direction 
@@ -715,5 +717,4 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  kh          !< eddy diffusivity for heat 
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  km          !< eddy diffusivity for momentum
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  l_wall      !< near-wall mixing length
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  prr         !< rain rate
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  p_loc       !< local array in multigrid/sor solver containing the pressure which is iteratively advanced in each iteration step
Index: palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 2912)
+++ palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 2918)
@@ -25,4 +25,15 @@
 ! -----------------
 ! $Id$
+! - defined l_grid only within this module
+! - Moved l_wall definition from modules.f90
+! - Get level of highest topography, used to limit upward distance calculation
+! - Consider cyclic boundary conditions for mixing length calculation
+! - Moved copy of wall_flags into subarray to subroutine
+! - Implemented l_wall calculation in case of RANS simulation
+! - Moved init of l_black to tcm_init_mixing_length
+! - Moved init_mixing_length from init_grid.f90 and
+!   renamed it to tcm_init_mixing_length
+! 
+! 2764 2018-01-22 09:25:36Z gronemeier
 ! Bugfix: remove duplicate SAVE statements
 ! 
@@ -66,10 +77,10 @@
 #if defined( __nopointer )
     USE arrays_3d,                                                             &
-        ONLY:  diss, diss_p, dzu, e, e_p, kh, km, l_grid,                      &
+        ONLY:  diss, diss_p, dzu, e, e_p, kh, km,                              &
                mean_inflow_profiles, prho, pt, tdiss_m, te_m, tend, u, v, vpt, w
 #else
     USE arrays_3d,                                                             &
         ONLY:  diss, diss_1, diss_2, diss_3, diss_p, dzu, e, e_1, e_2, e_3,    &
-               e_p, kh, km, l_grid, mean_inflow_profiles, prho, pt, tdiss_m,   &
+               e_p, kh, km, mean_inflow_profiles, prho, pt, tdiss_m,           &
                te_m, tend, u, v, vpt, w
 #endif
@@ -133,5 +144,8 @@
     INTEGER(iwp) ::  surf_s        !< start index of surface elements at given i-j position
 
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l_black      !< mixing length according to Blackadar
+    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(:,:,:), ALLOCATABLE :: dummy1 !< debug output variable
@@ -368,5 +382,5 @@
 
        CASE ( 'dummy2', 'dummy3', 'dummy1' )
-          unit = '?'
+          unit = 'mixing length'
 
        CASE ( 'kh', 'km' )
@@ -875,9 +889,6 @@
  SUBROUTINE tcm_init
 
-    USE arrays_3d,                                                             &
-        ONLY:  ug, vg, zu
-
     USE control_parameters,                                                    &
-        ONLY:  complex_terrain, f, kappa, dissipation_1d, topography
+        ONLY:  complex_terrain, dissipation_1d, topography
 
     USE model_1d_mod,                                                          &
@@ -894,4 +905,8 @@
     INTEGER(iwp) :: nz_s_shift   !<
     INTEGER(iwp) :: nz_s_shift_l !< 
+
+!
+!-- Initialize mixing length
+    CALL tcm_init_mixing_length
 
 !
@@ -1061,10 +1076,188 @@
     ENDIF
 
-!
-!-- Calculate mixing length according to Blackadar (1962)
-    IF ( rans_mode )  THEN
-
+ END SUBROUTINE tcm_init
+
+
+! Description:
+! -----------------------------------------------------------------------------!
+!> Pre-computation of grid-dependent and near-wall mixing length.
+!------------------------------------------------------------------------------!
+ SUBROUTINE tcm_init_mixing_length
+
+    USE arrays_3d,                                                             &
+        ONLY:  dzw, ug, vg, zu, zw
+
+    USE control_parameters,                                                    &
+        ONLY:  bc_lr_cyc, bc_ns_cyc, f, kappa, message_string,                 &
+               wall_adjustment_factor
+
+    USE grid_variables,                                                        &
+        ONLY:  dx, dy
+
+    USE indices,                                                               &
+        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb,  &
+               nzt, wall_flags_0
+
+    USE kinds
+
+
+    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) :: jj             !< 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(KIND=1), DIMENSION(:,:), ALLOCATABLE :: vic_yz !< contains a quarter of a single yz-slice of vicinity
+
+    INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE :: vicinity !< contains topography information of the vicinity of (i/j/k)
+
+    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: wall_flags_0_global !< wall_flags_0 of whole domain
+    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: wall_flags_dummy    !< dummy array required for MPI_ALLREDUCE command
+
+    REAL(wp) :: radius           !< search radius in meter
+
+    ALLOCATE( l_grid(1:nzt) )
+    ALLOCATE( l_wall(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
+       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)
+       ENDDO
+       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
+          ENDDO
+       ENDDO
+
+    ELSE
+!
+!-- Initialize the mixing length in case of a RANS simulation
+
+!
+!--    Calculate mixing length according to Blackadar (1962)
        IF ( f /= 0.0_wp )  THEN
-          l_max = 2.7E-4 * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) /        &
+          l_max = 2.7E-4 * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) /               &
                   ABS( f ) + 1E-10_wp
        ELSE
@@ -1078,7 +1271,391 @@
        l_black(nzt+1) = l_black(nzt)
 
-    ENDIF
-
- END SUBROUTINE tcm_init
+!
+!--    Gather topography information of whole domain
+       !## TODO: reduce amount of data sent by MPI call
+       !##  By now, a whole global 3D-array is sent and received with
+       !##  MPI_ALLREDUCE although most of the array is 0. This can be
+       !##  drastically reduced if only the local subarray is sent and stored
+       !##  in a global array. For that, an MPI data type or subarray must be
+       !##  defined.
+       !##  2018-03-19, gronemeier
+       ALLOCATE( wall_flags_0_global(nzb:nzt+1,0:ny,0:nx) )
+
+#if defined ( __parallel )
+       ALLOCATE( wall_flags_dummy(nzb:nzt+1,0:ny,0:nx) )
+       wall_flags_dummy = 0
+       wall_flags_dummy(nzb:nzt+1,nys:nyn,nxl:nxr) =  &
+           wall_flags_0(nzb:nzt+1,nys:nyn,nxl:nxr)
+
+       CALL MPI_ALLREDUCE( wall_flags_dummy,                  &
+                           wall_flags_0_global,               &
+                           (nzt-nzb+2)*(ny+1)*(nx+1),         &
+                           MPI_INTEGER, MPI_SUM, comm2d, ierr )
+       DEALLOCATE( wall_flags_dummy )
+#else
+       wall_flags_0_global(nzb:nzt+1,nys:nyn,nxl:nxr) =  &
+              wall_flags_0(nzb:nzt+1,nys:nyn,nxl:nxr)
+#endif
+!
+!--    Get height level of highest topography
+       DO  i = 0, nx
+          DO  j = 0, ny
+             DO  k = nzb+1, nzt-1
+                IF ( .NOT. BTEST( wall_flags_0_global(k,j,i), 0 ) .AND.  &
+                     k > k_max_topo )  &
+                   k_max_topo = k
+             ENDDO
+          ENDDO
+       ENDDO
+!
+!--    Limit mixing length to either nearest wall or Blackadar mixing length.
+!--    For that, analyze each grid point (i/j/k) ("analysed grid point") and 
+!--    search within its vicinity for the shortest distance to a wall by cal-
+!--    culating the distance between the analysed grid point and the "viewed
+!--    grid point" if it contains a wall (belongs to topography).
+       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
+
+!
+!--       Compute search radius as number of grid points in all directions
+          rad_i = CEILING( radius / dx )
+          rad_j = CEILING( radius / dy )
+
+          DO  kk = 0, nzt-k
+             rad_k_t = kk
+!
+!--          Limit upward search radius to height of maximum topography
+             IF ( zu(k+kk)-zu(k) >= radius .OR. k+kk >= k_max_topo )  EXIT
+          ENDDO
+
+          DO  kk = 0, k
+             rad_k_b = kk
+             IF ( zu(k)-zu(k-kk) >= radius )  EXIT
+          ENDDO
+
+!
+!--       Get maximum vertical radius; necessary for defining arrays
+          rad_k = MAX( rad_k_b, rad_k_t )
+!
+!--       When analysed grid point lies above maximum topography, set search
+!--       radius to 0 if the distance between the analysed grid point and max 
+!--       topography height is larger than the maximum search radius
+          IF ( zu(k-rad_k_b) > zu(k_max_topo) )  rad_k_b = 0
+!
+!--       Search within vicinity only if the vertical search radius is >0
+          IF ( rad_k_b /= 0 .OR. rad_k_t /= 0 )  THEN
+
+             !## NOTE: shape of vicinity is larger in z direction
+             !##  Shape of vicinity is two grid points larger than actual search
+             !##  radius in vertical direction. The first and last grid point is
+             !##   always set to 1 to asure correct detection of topography. See
+             !##  function "shortest_distance" for details.
+             !##  2018-03-16, gronemeier
+             ALLOCATE( vicinity(-rad_k-1:rad_k+1,-rad_j:rad_j,-rad_i:rad_i) )
+             ALLOCATE( vic_yz(0:rad_k+1,0:rad_j) )
+
+             vicinity = 1
+
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+!
+!--                Start search only if (i/j/k) belongs to atmosphere
+                   IF ( BTEST( wall_flags_0(k,j,i), 0 )  )  THEN
+!
+!--                   Reset topography within vicinity
+                      vicinity(-rad_k:rad_k,:,:) = 0
+!
+!--                   Copy area surrounding analysed grid point into vicinity.
+!--                   First, limit size of data copied to vicinity by the domain
+!--                   border
+                      rad_i_l = MIN( rad_i, i )
+                      rad_i_r = MIN( rad_i, nx-i )
+
+                      rad_j_s = MIN( rad_j, j )
+                      rad_j_n = MIN( rad_j, ny-j )
+
+                      CALL copy_into_vicinity( k, j, i,           &
+                                               -rad_k_b, rad_k_t, &
+                                               -rad_j_s, rad_j_n, &
+                                               -rad_i_l, rad_i_r  )
+!
+!--                   In case of cyclic boundaries, copy parts into vicinity
+!--                   where vicinity reaches over the domain borders.
+                      IF ( bc_lr_cyc )  THEN
+!
+!--                      Vicinity reaches over left domain boundary
+                         IF ( rad_i > rad_i_l )  THEN
+                            CALL copy_into_vicinity( k, j, nx+rad_i_l+1, &
+                                                     -rad_k_b, rad_k_t,  &
+                                                     -rad_j_s, rad_j_n,  &
+                                                     -rad_i, -rad_i_l-1  )
+!
+!--                         ...and over southern domain boundary
+                            IF ( bc_ns_cyc .AND. rad_j > rad_j_s )  &
+                               CALL copy_into_vicinity( k, ny+rad_j_s+1,    &
+                                                        nx+rad_i_l+1,       &
+                                                        -rad_k_b, rad_k_t,  &
+                                                        -rad_j, -rad_j_s-1, &
+                                                        -rad_i, -rad_i_l-1  )
+!
+!--                         ...and over northern domain boundary
+                            IF ( bc_ns_cyc .AND. rad_j > rad_j_n )  &
+                               CALL copy_into_vicinity( k, 0-rad_j_n-1,    &
+                                                        nx+rad_i_l+1,      &
+                                                        -rad_k_b, rad_k_t, &
+                                                         rad_j_n+1, rad_j, &
+                                                        -rad_i, -rad_i_l-1 )
+                         ENDIF
+!
+!--                      Vicinity reaches over right domain boundary
+                         IF ( rad_i > rad_i_r )  THEN
+                            CALL copy_into_vicinity( k, j, 0-rad_i_r-1, &
+                                                     -rad_k_b, rad_k_t, &
+                                                     -rad_j_s, rad_j_n, &
+                                                      rad_i_r+1, rad_i  )
+!
+!--                         ...and over southern domain boundary
+                            IF ( bc_ns_cyc .AND. rad_j > rad_j_s )  &
+                               CALL copy_into_vicinity( k, ny+rad_j_s+1,    &
+                                                        0-rad_i_r-1,        &
+                                                        -rad_k_b, rad_k_t,  &
+                                                        -rad_j, -rad_j_s-1, &
+                                                         rad_i_r+1, rad_i   )
+!
+!--                         ...and over northern domain boundary
+                            IF ( bc_ns_cyc .AND. rad_j > rad_j_n )  &
+                               CALL copy_into_vicinity( k, 0-rad_j_n-1,    &
+                                                        0-rad_i_r-1,       &
+                                                        -rad_k_b, rad_k_t, &
+                                                         rad_j_n+1, rad_j, &
+                                                         rad_i_r+1, rad_i  )
+                         ENDIF
+                      ENDIF
+
+                      IF ( bc_ns_cyc )  THEN
+!
+!--                      Vicinity reaches over southern domain boundary
+                         IF ( rad_j > rad_j_s )  &
+                            CALL copy_into_vicinity( k, ny+rad_j_s+1, i, &
+                                                     -rad_k_b, rad_k_t,  &
+                                                     -rad_j, -rad_j_s-1, &
+                                                     -rad_i_l, rad_i_r   )
+!
+!--                      Vicinity reaches over northern domain boundary
+                         IF ( rad_j > rad_j_n )  &
+                            CALL copy_into_vicinity( k, 0-rad_j_n-1, i, &
+                                                     -rad_k_b, rad_k_t, &
+                                                      rad_j_n+1, rad_j, &
+                                                      rad_i_l, rad_i_r  )
+                      ENDIF
+!
+!--                   Search for walls only if there is any within vicinity
+                      IF ( MAXVAL( vicinity(-rad_k:rad_k,:,:) ) /= 0 )  THEN
+!
+!--                      Search within first half (positive x)
+                         dist_dx = rad_i
+                         DO  ii = 0, dist_dx
+!
+!--                         Search along vertical direction only if below
+!--                         maximum topography
+                            IF ( rad_k_t > 0 ) THEN
+!
+!--                            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),             &
+                                       shortest_distance( vic_yz, .TRUE., ii ) )
+!
+!--                            Search for walls within octant (+-+)
+!--                            Switch order of array so that the analysed grid
+!--                            point is always located at (0/0) (required by
+!--                            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),             &
+                                       shortest_distance( vic_yz, .TRUE., ii ) )
+
+                            ENDIF
+!
+!--                         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),                &
+                                      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),                &
+                                      shortest_distance( vic_yz, .FALSE., ii ) )
+!
+!--                         Reduce search along x by already found distance
+                            dist_dx = CEILING( l_wall(k,j,i) / dx )
+
+                         ENDDO
+!
+!-                       Search within second half (negative x)
+                         DO  ii = 0, -dist_dx, -1
+!
+!--                         Search along vertical direction only if below
+!--                         maximum topography
+                            IF ( rad_k_t > 0 ) THEN
+!
+!--                            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),             &
+                                      shortest_distance( vic_yz, .TRUE., -ii ) )
+!
+!--                            Search for walls within octant (--+)
+!--                            Switch order of array so that the analysed grid
+!--                            point is always located at (0/0) (required by
+!--                            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),             &
+                                      shortest_distance( vic_yz, .TRUE., -ii ) )
+
+                            ENDIF
+!
+!--                         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),                &
+                                     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),                &
+                                     shortest_distance( vic_yz, .FALSE., -ii ) )
+!
+!--                         Reduce search along x by already found distance
+                            dist_dx = CEILING( l_wall(k,j,i) / dx )
+
+                         ENDDO
+
+                      ENDIF  !Check for any walls within vicinity
+
+                   ELSE  !Check if (i,j,k) belongs to atmosphere
+
+                      l_wall(k,j,i) = -999.0
+
+                   ENDIF
+
+                ENDDO  !j loop
+             ENDDO  !i loop
+
+             DEALLOCATE( vicinity )
+             DEALLOCATE( vic_yz )
+
+          ENDIF  !check vertical size of vicinity
+
+       ENDDO  !k loop
+
+       DEALLOCATE( wall_flags_0_global )
+
+    ENDIF  !LES or RANS mode
+
+!
+!-- Set lateral boundary conditions for l_wall
+    CALL exchange_horiz( l_wall, nbgp )
+
+    CONTAINS
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculate the shortest distance between position (i/j/k)=(0/0/0) and 
+!> (pos_i/jj/kk), where (jj/kk) is the position of the maximum of 'array'
+!> closest to the origin (0/0) of 'array'.
+!------------------------------------------------------------------------------!
+    REAL FUNCTION shortest_distance( array, orientation, pos_i )
+
+       IMPLICIT NONE
+
+       LOGICAL, INTENT(IN) :: orientation    !< flag if array represents an array oriented upwards (true) or downwards (false)
+
+       INTEGER(iwp), INTENT(IN) :: pos_i     !< x position of the yz-plane 'array'
+
+       INTEGER(iwp) :: jj                    !< loop index
+
+       INTEGER(iwp), DIMENSION(0:rad_j) :: loc_k !< location of closest wall along vertical dimension
+
+       INTEGER(KIND=1), DIMENSION(0:rad_k+1,0:rad_j), INTENT(IN) :: array !< array containing a yz-plane at position pos_i
+
+!
+!--    Get coordinate of first maximum along vertical dimension
+!--    at each y position of array.
+!--    Substract 1 because indices count from 1 instead of 0 by MAXLOC
+       loc_k = MAXLOC( array, DIM = 1) - 1
+
+!
+!--    Set distance to the default maximum value (=search radius)
+       shortest_distance = radius
+!
+!--    Calculate distance between position (0/0/0) and
+!--    position (pos_i/jj/loc(jj)) and only save the shortest distance.
+       IF ( orientation ) THEN  !if array is oriented upwards
+          DO  jj = 0, rad_j
+             shortest_distance = MIN( shortest_distance,                        &
+                                      SQRT( MAX(pos_i*dx-0.5*dx,0.0)**2         &
+                                          + MAX(jj*dy-0.5*dy,0.0)**2            &
+                                          + MAX(zw(loc_k(jj)+k-1)-zu(k),0.0)**2 &
+                                          )                                     &
+                                    )
+          ENDDO
+       ELSE  !if array is oriented downwards
+          !## NOTE: MAX within zw required to circumvent error at domain border
+          !##  At the domain border, if non-cyclic boundary is present, the
+          !##  index for zw could be -1, which will be errorneous (zw(-1) does
+          !##  not exist). The MAX function limits the index to be at least 0.
+          DO  jj = 0, rad_j
+             shortest_distance = MIN( shortest_distance,                       &
+                                      SQRT( MAX(pos_i*dx-0.5*dx,0.0)**2        &
+                                          + MAX(jj*dy-0.5*dy,0.0)**2           &
+                                          + MAX(zu(k)-zw(MAX(k-loc_k(jj),      &
+                                                             0_iwp)),          &
+                                                0.0)**2                        &
+                                          )                                    &
+                                    )
+          ENDDO
+       ENDIF
+
+    END FUNCTION
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Copy a subarray of size (kb:kt,js:jn,il:ir) centered around grid point 
+!> (kp,jp,ip) containing the first bit of wall_flags_0 into the array
+!> 'vicinity'. Only copy first bit as this indicates the presence of topography.
+!------------------------------------------------------------------------------!
+    SUBROUTINE copy_into_vicinity( kp, jp, ip, kb, kt, js, jn, il, ir )
+
+       IMPLICIT NONE
+
+       INTEGER(iwp), INTENT(IN) :: il !< left loop boundary
+       INTEGER(iwp), INTENT(IN) :: ip !< center position in x-direction
+       INTEGER(iwp), INTENT(IN) :: ir !< right loop boundary
+       INTEGER(iwp), INTENT(IN) :: jn !< northern loop boundary
+       INTEGER(iwp), INTENT(IN) :: jp !< center position in y-direction
+       INTEGER(iwp), INTENT(IN) :: js !< southern loop boundary
+       INTEGER(iwp), INTENT(IN) :: kb !< bottom loop boundary
+       INTEGER(iwp), INTENT(IN) :: kp !< center position in z-direction
+       INTEGER(iwp), INTENT(IN) :: kt !< top loop boundary
+
+       INTEGER(iwp) :: i   !< loop index
+       INTEGER(iwp) :: j   !< loop index
+       INTEGER(iwp) :: k   !< loop index
+
+
+       DO  i = il, ir
+          DO  j = js, jn
+             DO  k = kb, kt
+                vicinity(k,j,i) = MERGE( 0, 1,               &
+                       BTEST( wall_flags_0_global(kp+k,jp+j,ip+i), 0 ) )
+             ENDDO
+          ENDDO
+       ENDDO
+
+    END SUBROUTINE copy_into_vicinity
+
+ END SUBROUTINE tcm_init_mixing_length
 
 
@@ -3531,5 +4108,5 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  dd2zu, l_grid, l_wall
+        ONLY:  dd2zu
 
     USE control_parameters,                                                    &
@@ -3630,5 +4207,5 @@
 !-- Calculate diabatic mixing length using Dyer-profile functions
     IF ( rif >= 0.0_wp )  THEN
-       l      = l_black(k) / ( 1.0_wp + 5.0_wp * rif )
+       l      = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif ), l_wall(k,j,i) )
        l_diss = l
     ELSE
@@ -3636,6 +4213,6 @@
 !--    In case of unstable stratification, use mixing length of neutral case
 !--    for l, but consider profile functions for l_diss
-       l      = l_black(k)
-       l_diss = l_black(k) * SQRT( 1.0_wp - 16.0_wp * rif )
+       l      = l_wall(k,j,i)
+       l_diss = l * SQRT( 1.0_wp - 16.0_wp * rif )
     ENDIF
 
