Index: /palm/trunk/SOURCE/lpm.f90
===================================================================
--- /palm/trunk/SOURCE/lpm.f90	(revision 2416)
+++ /palm/trunk/SOURCE/lpm.f90	(revision 2417)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! Major bugfixes in modeling SGS particle speeds (since revision 1359).
+! Particle sorting added to distinguish between already completed and 
+! non-completed particles.
+! 
+! 2263 2017-06-08 14:59:01Z schwenkel
 ! Implemented splitting and merging algorithm
 ! 
@@ -140,5 +145,5 @@
 
     USE lpm_pack_arrays_mod,                                                   &
-        ONLY:  lpm_pack_all_arrays
+        ONLY:  lpm_pack_all_arrays, lpm_sort
 
     USE particle_attributes,                                                   &
@@ -259,5 +264,4 @@
        CALL cpu_log( log_point_s(44), 'lpm_advec', 'pause' )
        
-       grid_particles(:,:,:)%time_loop_done = .TRUE.
 !
 !--    If particle advection includes SGS velocity components, calculate the
@@ -265,8 +269,17 @@
 !--    horizontally averaged profiles of the SGS TKE and the resolved-scale 
 !--    velocity variances)
-
        IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
           CALL lpm_init_sgs_tke
        ENDIF
+
+!
+!--    In case SGS-particle speed is considered, particles may carry out 
+!--    several particle timesteps. In order to prevent unnecessary 
+!--    treatment of particles that already reached the final time level, 
+!--    particles are sorted into contiguous blocks of finished and 
+!--    not-finished particles, in addition to their already sorting
+!--    according to their sub-boxes. 
+       IF ( .NOT. first_loop_stride  .AND.  use_sgs_for_particles )            &
+          CALL lpm_sort
 
        DO  i = nxl, nxr
@@ -323,5 +336,5 @@
 !--             optimization is still possible.)
                 IF ( topography /= 'flat' .AND. k < nzb_max + 2 )  THEN
-                   CALL lpm_boundary_conds( 'walls' )
+                   CALL lpm_boundary_conds( 'walls' )!, i, j, k )
                 ENDIF
 !
@@ -333,18 +346,22 @@
 !--             the top or bottom boundary and delete those particles, which are
 !--             older than allowed
-                CALL lpm_boundary_conds( 'bottom/top' )
+                CALL lpm_boundary_conds( 'bottom/top' ) !, i, j, k )
 !
 !---            If not all particles of the actual grid cell have reached the 
-!--             LES timestep, this cell has to to another loop iteration. Due to
-!--             the fact that particles can move into neighboring grid cell, 
-!--             these neighbor cells also have to perform another loop iteration
+!--             LES timestep, this cell has to do another loop iteration. Due to
+!--             the fact that particles can move into neighboring grid cells, 
+!--             these neighbor cells also have to perform another loop iteration.
+!--             Please note, this realization does not work properly if 
+!--             particles move into another subdomain. 
                 IF ( .NOT. dt_3d_reached_l )  THEN
-                   ks = MAX(nzb+1,k)
-                   ke = MIN(nzt,k)
-                   js = MAX(nys,j)
-                   je = MIN(nyn,j)
-                   is = MAX(nxl,i)
-                   ie = MIN(nxr,i)
+                   ks = MAX(nzb+1,k-1)
+                   ke = MIN(nzt,k+1)
+                   js = MAX(nys,j-1)
+                   je = MIN(nyn,j+1)
+                   is = MAX(nxl,i-1)
+                   ie = MIN(nxr,i+1)
                    grid_particles(ks:ke,js:je,is:ie)%time_loop_done = .FALSE.
+                ELSE
+                   grid_particles(k,j,i)%time_loop_done = .TRUE.
                 ENDIF
 
@@ -392,5 +409,4 @@
 !--    determine new number of particles
        CALL lpm_pack_all_arrays
-
 !
 !--    Initialize variables for the next (sub-) timestep, i.e., for marking 
Index: /palm/trunk/SOURCE/lpm_advec.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_advec.f90	(revision 2416)
+++ /palm/trunk/SOURCE/lpm_advec.f90	(revision 2417)
@@ -25,4 +25,10 @@
 ! -----------------
 ! $Id$
+! Particle loops adapted for sub-box structure, i.e. for each sub-box the 
+! particle loop runs from start_index up to end_index instead from 1 to 
+! number_of_particles. This way, it is possible to skip unnecessary 
+! computations for particles that already completed the LES timestep. 
+! 
+! 2318 2017-07-20 17:27:44Z suehring
 ! Get topography top index via Function call 
 ! 
@@ -267,5 +273,5 @@
 
     DO  nb = 0, 7
-
+       
        i = ip
        j = jp + block_offset(nb)%j_off
@@ -287,5 +293,5 @@
 !
 !--       Determine vertical index of topography top
-          k_wall = get_topography_top_index( jlog,ilog, 's' )
+          k_wall = get_topography_top_index( jlog, ilog, 's' )
 
           IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
@@ -648,309 +654,110 @@
        ELSE ! non-flat topography, e.g., buildings 
 
-          DO  n = 1, number_of_particles
-             i = particles(n)%x * ddx
-             j = particles(n)%y * ddy
-             k = ( zv(n) + 0.5_wp * dz * atmos_ocean_sign ) / dz  &
-                 + offset_ocean_nzt                      ! only exact if eq.dist
-!
-!--          In case that there are buildings it has to be determined 
-!--          how many of the gridpoints defining the particle box are 
-!--          situated within a building
-!--          gp_outside_of_building(1): i,j,k
-!--          gp_outside_of_building(2): i,j+1,k
-!--          gp_outside_of_building(3): i,j,k+1
-!--          gp_outside_of_building(4): i,j+1,k+1
-!--          gp_outside_of_building(5): i+1,j,k
-!--          gp_outside_of_building(6): i+1,j+1,k
-!--          gp_outside_of_building(7): i+1,j,k+1
-!--          gp_outside_of_building(8): i+1,j+1,k+1
-
-             gp_outside_of_building = 0
-             location = 0.0_wp
-             num_gp = 0
-
-!
-!--          Determine vertical index of topography top at (j,i)
-             k_wall = get_topography_top_index( j, i, 's' )
-!
-!--          To do: Reconsider order of computations in order to avoid 
-!--          unnecessary index calculations.
-             IF ( k > k_wall  .OR.  k_wall == 0 )  THEN 
-                num_gp = num_gp + 1
-                gp_outside_of_building(1) = 1
-                location(num_gp,1) = i * dx
-                location(num_gp,2) = j * dy
-                location(num_gp,3) = k * dz - 0.5_wp * dz
-                ei(num_gp)     = e(k,j,i)
-                dissi(num_gp)  = diss(k,j,i)
-                de_dxi(num_gp) = de_dx(k,j,i)
-                de_dyi(num_gp) = de_dy(k,j,i)
-                de_dzi(num_gp) = de_dz(k,j,i)
-             ENDIF
-
-!
-!--          Determine vertical index of topography top at (j+1,i)
-             k_wall = get_topography_top_index( j+1, i, 's' )
-
-             IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
-                num_gp = num_gp + 1
-                gp_outside_of_building(2) = 1
-                location(num_gp,1) = i * dx
-                location(num_gp,2) = (j+1) * dy
-                location(num_gp,3) = k * dz - 0.5_wp * dz
-                ei(num_gp)     = e(k,j+1,i)
-                dissi(num_gp)  = diss(k,j+1,i)
-                de_dxi(num_gp) = de_dx(k,j+1,i)
-                de_dyi(num_gp) = de_dy(k,j+1,i)
-                de_dzi(num_gp) = de_dz(k,j+1,i)
-             ENDIF
-
-!
-!--          Determine vertical index of topography top at (j,i)
-             k_wall = get_topography_top_index( j, i, 's' )
-
-             IF ( k+1 > k_wall  .OR.  k_wall == 0 )  THEN
-                num_gp = num_gp + 1
-                gp_outside_of_building(3) = 1
-                location(num_gp,1) = i * dx
-                location(num_gp,2) = j * dy
-                location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
-                ei(num_gp)     = e(k+1,j,i)
-                dissi(num_gp)  = diss(k+1,j,i)
-                de_dxi(num_gp) = de_dx(k+1,j,i)
-                de_dyi(num_gp) = de_dy(k+1,j,i)
-                de_dzi(num_gp) = de_dz(k+1,j,i)
-             ENDIF
-
-!
-!--          Determine vertical index of topography top at (j+1,i)
-             k_wall = get_topography_top_index( j+1, i, 's' )
-             IF ( k+1 > k_wall  .OR.  k_wall == 0 )  THEN
-                num_gp = num_gp + 1
-                gp_outside_of_building(4) = 1
-                location(num_gp,1) = i * dx
-                location(num_gp,2) = (j+1) * dy
-                location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
-                ei(num_gp)     = e(k+1,j+1,i)
-                dissi(num_gp)  = diss(k+1,j+1,i)
-                de_dxi(num_gp) = de_dx(k+1,j+1,i)
-                de_dyi(num_gp) = de_dy(k+1,j+1,i)
-                de_dzi(num_gp) = de_dz(k+1,j+1,i)
-             ENDIF
-
-!
-!--          Determine vertical index of topography top at (j,i+1)
-             k_wall = get_topography_top_index( j, i+1, 's' )
-             IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
-                num_gp = num_gp + 1
-                gp_outside_of_building(5) = 1
-                location(num_gp,1) = (i+1) * dx
-                location(num_gp,2) = j * dy
-                location(num_gp,3) = k * dz - 0.5_wp * dz
-                ei(num_gp)     = e(k,j,i+1)
-                dissi(num_gp)  = diss(k,j,i+1)
-                de_dxi(num_gp) = de_dx(k,j,i+1)
-                de_dyi(num_gp) = de_dy(k,j,i+1)
-                de_dzi(num_gp) = de_dz(k,j,i+1)
-             ENDIF
-
-!
-!--          Determine vertical index of topography top at (j+1,i+1)
-             k_wall = get_topography_top_index( j+1, i+1, 's' )
-
-             IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
-                num_gp = num_gp + 1
-                gp_outside_of_building(6) = 1
-                location(num_gp,1) = (i+1) * dx
-                location(num_gp,2) = (j+1) * dy
-                location(num_gp,3) = k * dz - 0.5_wp * dz
-                ei(num_gp)     = e(k,j+1,i+1)
-                dissi(num_gp)  = diss(k,j+1,i+1)
-                de_dxi(num_gp) = de_dx(k,j+1,i+1)
-                de_dyi(num_gp) = de_dy(k,j+1,i+1)
-                de_dzi(num_gp) = de_dz(k,j+1,i+1)
-             ENDIF
-
-!
-!--          Determine vertical index of topography top at (j,i+1)
-             k_wall = get_topography_top_index( j, i+1, 's' )
-
-             IF ( k+1 > k_wall  .OR.  k_wall == 0 )  THEN
-                num_gp = num_gp + 1
-                gp_outside_of_building(7) = 1
-                location(num_gp,1) = (i+1) * dx
-                location(num_gp,2) = j * dy
-                location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
-                ei(num_gp)     = e(k+1,j,i+1)
-                dissi(num_gp)  = diss(k+1,j,i+1)
-                de_dxi(num_gp) = de_dx(k+1,j,i+1)
-                de_dyi(num_gp) = de_dy(k+1,j,i+1)
-                de_dzi(num_gp) = de_dz(k+1,j,i+1)
-             ENDIF
-
-!
-!--          Determine vertical index of topography top at (j+1,i+1)
-             k_wall = get_topography_top_index( j+1, i+1, 's' )
-
-             IF ( k+1 > k_wall  .OR.  k_wall == 0)  THEN
-                num_gp = num_gp + 1
-                gp_outside_of_building(8) = 1
-                location(num_gp,1) = (i+1) * dx
-                location(num_gp,2) = (j+1) * dy
-                location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
-                ei(num_gp)     = e(k+1,j+1,i+1)
-                dissi(num_gp)  = diss(k+1,j+1,i+1)
-                de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
-                de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
-                de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
-             ENDIF
-!
-!--          If all gridpoints are situated outside of a building, then the
-!--          ordinary interpolation scheme can be used.
-             IF ( num_gp == 8 )  THEN
-
-                x  = particles(n)%x - i * dx
-                y  = particles(n)%y - j * dy
-                aa = x**2          + y**2
-                bb = ( dx - x )**2 + y**2
-                cc = x**2          + ( dy - y )**2
-                dd = ( dx - x )**2 + ( dy - y )**2
-                gg = aa + bb + cc + dd
-  
-                e_int_l = ( ( gg - aa ) * e(k,j,i)   + ( gg - bb ) * e(k,j,i+1)   &
-                          + ( gg - cc ) * e(k,j+1,i) + ( gg - dd ) * e(k,j+1,i+1) &
-                          ) / ( 3.0_wp * gg )
-   
-                IF ( k == nzt )  THEN
-                   e_int(n) = e_int_l
-                ELSE
-                   e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
-                               ( gg - bb ) * e(k+1,j,i+1) + &
-                               ( gg - cc ) * e(k+1,j+1,i) + &
-                               ( gg - dd ) * e(k+1,j+1,i+1) &
-                             ) / ( 3.0_wp * gg )
-                   e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
-                                       ( e_int_u - e_int_l )
-                ENDIF
-! 
-!--             Needed to avoid NaN particle velocities (this might not be
-!--             required any more)
-                IF ( e_int(n) <= 0.0_wp )  THEN
-                   e_int(n) = 1.0E-20_wp
-                ENDIF
-!
-!--             Interpolate the TKE gradient along x (adopt incides i,j,k
-!--             and all position variables from above (TKE))
-                de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
-                                ( gg - bb ) * de_dx(k,j,i+1) + &
-                                ( gg - cc ) * de_dx(k,j+1,i) + &
-                                ( gg - dd ) * de_dx(k,j+1,i+1) &
-                              ) / ( 3.0_wp * gg )
-
-                IF ( ( k == nzt )  .OR.  ( k == nzb ) )  THEN
-                   de_dx_int(n) = de_dx_int_l
-                ELSE
-                   de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
-                                   ( gg - bb ) * de_dx(k+1,j,i+1) + &
-                                   ( gg - cc ) * de_dx(k+1,j+1,i) + &
-                                   ( gg - dd ) * de_dx(k+1,j+1,i+1) &
-                                 ) / ( 3.0_wp * gg )
-                   de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / &
-                                           dz * ( de_dx_int_u - de_dx_int_l )
-                ENDIF
-
-!
-!--             Interpolate the TKE gradient along y
-                de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
-                                ( gg - bb ) * de_dy(k,j,i+1) + &
-                                ( gg - cc ) * de_dy(k,j+1,i) + &
-                                ( gg - dd ) * de_dy(k,j+1,i+1) &
-                              ) / ( 3.0_wp * gg )
-                IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
-                   de_dy_int(n) = de_dy_int_l
-                ELSE
-                   de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
-                                   ( gg - bb ) * de_dy(k+1,j,i+1) + &
-                                   ( gg - cc ) * de_dy(k+1,j+1,i) + &
-                                   ( gg - dd ) * de_dy(k+1,j+1,i+1) &
-                                 ) / ( 3.0_wp * gg )
-                   de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / &
-                                           dz * ( de_dy_int_u - de_dy_int_l )
-                ENDIF
-
-!
-!--             Interpolate the TKE gradient along z
-                IF ( zv(n) < 0.5_wp * dz )  THEN
-                   de_dz_int(n) = 0.0_wp
-                ELSE
-                   de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
-                                   ( gg - bb ) * de_dz(k,j,i+1) + &
-                                   ( gg - cc ) * de_dz(k,j+1,i) + &
-                                   ( gg - dd ) * de_dz(k,j+1,i+1) &
-                                 ) / ( 3.0_wp * gg )
- 
-                   IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
-                      de_dz_int(n) = de_dz_int_l
-                   ELSE
-                      de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
-                                      ( gg - bb ) * de_dz(k+1,j,i+1) + &
-                                      ( gg - cc ) * de_dz(k+1,j+1,i) + &
-                                      ( gg - dd ) * de_dz(k+1,j+1,i+1) &
-                                    ) / ( 3.0_wp * gg )
-                      de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) /&
-                                           dz * ( de_dz_int_u - de_dz_int_l )
-                   ENDIF
-                ENDIF
-
-!
-!--             Interpolate the dissipation of TKE
-                diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
-                               ( gg - bb ) * diss(k,j,i+1) + &
-                               ( gg - cc ) * diss(k,j+1,i) + &
-                               ( gg - dd ) * diss(k,j+1,i+1) &
-                             ) / ( 3.0_wp * gg )
- 
-                IF ( k == nzt )  THEN
-                   diss_int(n) = diss_int_l
-                ELSE
-                   diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
-                                  ( gg - bb ) * diss(k+1,j,i+1) + &
-                                  ( gg - cc ) * diss(k+1,j+1,i) + &
-                                  ( gg - dd ) * diss(k+1,j+1,i+1) &
-                                ) / ( 3.0_wp * gg )
-                   diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz *&
-                                           ( diss_int_u - diss_int_l )
-                ENDIF
-!
-!--             Set flag for stochastic equation.
-                term_1_2(n) = 1.0_wp
- 
-             ELSE
- 
-!
-!--             If wall between gridpoint 1 and gridpoint 5, then
-!--             Neumann boundary condition has to be applied
-                IF ( gp_outside_of_building(1) == 1  .AND. &
-                     gp_outside_of_building(5) == 0 )  THEN
+          DO  nb = 0, 7
+
+             DO  n = start_index(nb), end_index(nb)
+
+                i = particles(n)%x * ddx
+                j = particles(n)%y * ddy
+                k = ( zv(n) + 0.5_wp * dz * atmos_ocean_sign ) / dz  &
+                    + offset_ocean_nzt                      ! only exact if eq.dist
+!
+!--             In case that there are buildings it has to be determined 
+!--             how many of the gridpoints defining the particle box are 
+!--             situated within a building
+!--             gp_outside_of_building(1): i,j,k
+!--             gp_outside_of_building(2): i,j+1,k
+!--             gp_outside_of_building(3): i,j,k+1
+!--             gp_outside_of_building(4): i,j+1,k+1
+!--             gp_outside_of_building(5): i+1,j,k
+!--             gp_outside_of_building(6): i+1,j+1,k
+!--             gp_outside_of_building(7): i+1,j,k+1
+!--             gp_outside_of_building(8): i+1,j+1,k+1
+
+                gp_outside_of_building = 0
+                location = 0.0_wp
+                num_gp = 0
+
+!
+!--             Determine vertical index of topography top at (j,i)
+                k_wall = get_topography_top_index( j, i, 's' )
+!
+!--             To do: Reconsider order of computations in order to avoid 
+!--             unnecessary index calculations.
+                IF ( k > k_wall  .OR.  k_wall == 0 )  THEN 
                    num_gp = num_gp + 1
-                   location(num_gp,1) = i * dx + 0.5_wp * dx
+                   gp_outside_of_building(1) = 1
+                   location(num_gp,1) = i * dx
                    location(num_gp,2) = j * dy
                    location(num_gp,3) = k * dz - 0.5_wp * dz
                    ei(num_gp)     = e(k,j,i)
                    dissi(num_gp)  = diss(k,j,i)
-                   de_dxi(num_gp) = 0.0_wp
+                   de_dxi(num_gp) = de_dx(k,j,i)
                    de_dyi(num_gp) = de_dy(k,j,i)
                    de_dzi(num_gp) = de_dz(k,j,i)
                 ENDIF
 
-                IF ( gp_outside_of_building(5) == 1  .AND. &
-                     gp_outside_of_building(1) == 0 )  THEN
+!
+!--             Determine vertical index of topography top at (j+1,i)
+                k_wall = get_topography_top_index( j+1, i, 's' )
+
+                IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
                    num_gp = num_gp + 1
-                   location(num_gp,1) = i * dx + 0.5_wp * dx
+                   gp_outside_of_building(2) = 1
+                   location(num_gp,1) = i * dx
+                   location(num_gp,2) = (j+1) * dy
+                   location(num_gp,3) = k * dz - 0.5_wp * dz
+                   ei(num_gp)     = e(k,j+1,i)
+                   dissi(num_gp)  = diss(k,j+1,i)
+                   de_dxi(num_gp) = de_dx(k,j+1,i)
+                   de_dyi(num_gp) = de_dy(k,j+1,i)
+                   de_dzi(num_gp) = de_dz(k,j+1,i)
+                ENDIF
+
+!
+!--             Determine vertical index of topography top at (j,i)
+                k_wall = get_topography_top_index( j, i, 's' )
+
+                IF ( k+1 > k_wall  .OR.  k_wall == 0 )  THEN
+                   num_gp = num_gp + 1
+                   gp_outside_of_building(3) = 1
+                   location(num_gp,1) = i * dx
+                   location(num_gp,2) = j * dy
+                   location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
+                   ei(num_gp)     = e(k+1,j,i)
+                   dissi(num_gp)  = diss(k+1,j,i)
+                   de_dxi(num_gp) = de_dx(k+1,j,i)
+                   de_dyi(num_gp) = de_dy(k+1,j,i)
+                   de_dzi(num_gp) = de_dz(k+1,j,i)
+                ENDIF
+
+!
+!--             Determine vertical index of topography top at (j+1,i)
+                k_wall = get_topography_top_index( j+1, i, 's' )
+                IF ( k+1 > k_wall  .OR.  k_wall == 0 )  THEN
+                   num_gp = num_gp + 1
+                   gp_outside_of_building(4) = 1
+                   location(num_gp,1) = i * dx
+                   location(num_gp,2) = (j+1) * dy
+                   location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
+                   ei(num_gp)     = e(k+1,j+1,i)
+                   dissi(num_gp)  = diss(k+1,j+1,i)
+                   de_dxi(num_gp) = de_dx(k+1,j+1,i)
+                   de_dyi(num_gp) = de_dy(k+1,j+1,i)
+                   de_dzi(num_gp) = de_dz(k+1,j+1,i)
+                ENDIF
+
+!
+!--             Determine vertical index of topography top at (j,i+1)
+                k_wall = get_topography_top_index( j, i+1, 's' )
+                IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
+                   num_gp = num_gp + 1
+                   gp_outside_of_building(5) = 1
+                   location(num_gp,1) = (i+1) * dx
                    location(num_gp,2) = j * dy
                    location(num_gp,3) = k * dz - 0.5_wp * dz
                    ei(num_gp)     = e(k,j,i+1)
                    dissi(num_gp)  = diss(k,j,i+1)
-                   de_dxi(num_gp) = 0.0_wp
+                   de_dxi(num_gp) = de_dx(k,j,i+1)
                    de_dyi(num_gp) = de_dy(k,j,i+1)
                    de_dzi(num_gp) = de_dz(k,j,i+1)
@@ -958,57 +765,16 @@
 
 !
-!--             If wall between gridpoint 5 and gridpoint 6, then
-!--             then Neumann boundary condition has to be applied
-                IF ( gp_outside_of_building(5) == 1  .AND. &
-                     gp_outside_of_building(6) == 0 )  THEN
+!--             Determine vertical index of topography top at (j+1,i+1)
+                k_wall = get_topography_top_index( j+1, i+1, 's' )
+
+                IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
                    num_gp = num_gp + 1
+                   gp_outside_of_building(6) = 1
                    location(num_gp,1) = (i+1) * dx
-                   location(num_gp,2) = j * dy + 0.5_wp * dy
-                   location(num_gp,3) = k * dz - 0.5_wp * dz
-                   ei(num_gp)     = e(k,j,i+1)
-                   dissi(num_gp)  = diss(k,j,i+1)
-                   de_dxi(num_gp) = de_dx(k,j,i+1)
-                   de_dyi(num_gp) = 0.0_wp
-                   de_dzi(num_gp) = de_dz(k,j,i+1)
-                ENDIF
-
-                IF ( gp_outside_of_building(6) == 1  .AND. &
-                     gp_outside_of_building(5) == 0 )  THEN
-                   num_gp = num_gp + 1
-                   location(num_gp,1) = (i+1) * dx
-                   location(num_gp,2) = j * dy + 0.5_wp * dy
+                   location(num_gp,2) = (j+1) * dy
                    location(num_gp,3) = k * dz - 0.5_wp * dz
                    ei(num_gp)     = e(k,j+1,i+1)
                    dissi(num_gp)  = diss(k,j+1,i+1)
                    de_dxi(num_gp) = de_dx(k,j+1,i+1)
-                   de_dyi(num_gp) = 0.0_wp
-                   de_dzi(num_gp) = de_dz(k,j+1,i+1)
-                ENDIF
-
-!
-!--             If wall between gridpoint 2 and gridpoint 6, then
-!--             Neumann boundary condition has to be applied
-                IF ( gp_outside_of_building(2) == 1  .AND. &
-                     gp_outside_of_building(6) == 0 )  THEN
-                   num_gp = num_gp + 1
-                   location(num_gp,1) = i * dx + 0.5_wp * dx
-                   location(num_gp,2) = (j+1) * dy
-                   location(num_gp,3) = k * dz - 0.5_wp * dz
-                   ei(num_gp)     = e(k,j+1,i)
-                   dissi(num_gp)  = diss(k,j+1,i)
-                   de_dxi(num_gp) = 0.0_wp
-                   de_dyi(num_gp) = de_dy(k,j+1,i)
-                   de_dzi(num_gp) = de_dz(k,j+1,i)
-                ENDIF
-
-                IF ( gp_outside_of_building(6) == 1  .AND. &
-                     gp_outside_of_building(2) == 0 )  THEN
-                   num_gp = num_gp + 1
-                   location(num_gp,1) = i * dx + 0.5_wp * dx
-                   location(num_gp,2) = (j+1) * dy
-                   location(num_gp,3) = k * dz - 0.5_wp * dz
-                   ei(num_gp)     = e(k,j+1,i+1)
-                   dissi(num_gp)  = diss(k,j+1,i+1)
-                   de_dxi(num_gp) = 0.0_wp
                    de_dyi(num_gp) = de_dy(k,j+1,i+1)
                    de_dzi(num_gp) = de_dz(k,j+1,i+1)
@@ -1016,277 +782,521 @@
 
 !
-!--             If wall between gridpoint 1 and gridpoint 2, then
-!--             Neumann boundary condition has to be applied
-                IF ( gp_outside_of_building(1) == 1  .AND. &
-                     gp_outside_of_building(2) == 0 )  THEN
+!--             Determine vertical index of topography top at (j,i+1)
+                k_wall = get_topography_top_index( j, i+1, 's' )
+
+                IF ( k+1 > k_wall  .OR.  k_wall == 0 )  THEN
                    num_gp = num_gp + 1
-                   location(num_gp,1) = i * dx
-                   location(num_gp,2) = j * dy + 0.5_wp * dy
-                   location(num_gp,3) = k * dz - 0.5_wp * dz
-                   ei(num_gp)     = e(k,j,i)
-                   dissi(num_gp)  = diss(k,j,i)
-                   de_dxi(num_gp) = de_dx(k,j,i)
-                   de_dyi(num_gp) = 0.0_wp
-                   de_dzi(num_gp) = de_dz(k,j,i)  
-                ENDIF
-
-                IF ( gp_outside_of_building(2) == 1  .AND. &
-                     gp_outside_of_building(1) == 0 )  THEN
-                   num_gp = num_gp + 1
-                   location(num_gp,1) = i * dx
-                   location(num_gp,2) = j * dy + 0.5_wp * dy
-                   location(num_gp,3) = k * dz - 0.5_wp * dz
-                   ei(num_gp)     = e(k,j+1,i)
-                   dissi(num_gp)  = diss(k,j+1,i)
-                   de_dxi(num_gp) = de_dx(k,j+1,i)
-                   de_dyi(num_gp) = 0.0_wp
-                   de_dzi(num_gp) = de_dz(k,j+1,i)
-                ENDIF
-
-!
-!--             If wall between gridpoint 3 and gridpoint 7, then
-!--             Neumann boundary condition has to be applied
-                IF ( gp_outside_of_building(3) == 1  .AND. &
-                     gp_outside_of_building(7) == 0 )  THEN
-                   num_gp = num_gp + 1
-                   location(num_gp,1) = i * dx + 0.5_wp * dx
-                   location(num_gp,2) = j * dy
-                   location(num_gp,3) = k * dz + 0.5_wp * dz
-                   ei(num_gp)     = e(k+1,j,i)
-                   dissi(num_gp)  = diss(k+1,j,i)
-                   de_dxi(num_gp) = 0.0_wp
-                   de_dyi(num_gp) = de_dy(k+1,j,i)
-                   de_dzi(num_gp) = de_dz(k+1,j,i) 
-                ENDIF
-
-                IF ( gp_outside_of_building(7) == 1  .AND. &
-                     gp_outside_of_building(3) == 0 )  THEN
-                   num_gp = num_gp + 1
-                   location(num_gp,1) = i * dx + 0.5_wp * dx
-                   location(num_gp,2) = j * dy
-                   location(num_gp,3) = k * dz + 0.5_wp * dz
-                   ei(num_gp)     = e(k+1,j,i+1)
-                   dissi(num_gp)  = diss(k+1,j,i+1)
-                   de_dxi(num_gp) = 0.0_wp
-                   de_dyi(num_gp) = de_dy(k+1,j,i+1)
-                   de_dzi(num_gp) = de_dz(k+1,j,i+1)
-                ENDIF
-
-!
-!--             If wall between gridpoint 7 and gridpoint 8, then
-!--             Neumann boundary condition has to be applied
-                IF ( gp_outside_of_building(7) == 1  .AND. &
-                     gp_outside_of_building(8) == 0 )  THEN
-                   num_gp = num_gp + 1
-                   location(num_gp,1) = (i+1) * dx
-                   location(num_gp,2) = j * dy + 0.5_wp * dy
-                   location(num_gp,3) = k * dz + 0.5_wp * dz
-                   ei(num_gp)     = e(k+1,j,i+1)
-                   dissi(num_gp)  = diss(k+1,j,i+1)
-                   de_dxi(num_gp) = de_dx(k+1,j,i+1)
-                   de_dyi(num_gp) = 0.0_wp
-                   de_dzi(num_gp) = de_dz(k+1,j,i+1)
-                ENDIF
-
-                IF ( gp_outside_of_building(8) == 1  .AND. &
-                     gp_outside_of_building(7) == 0 )  THEN
-                   num_gp = num_gp + 1
-                   location(num_gp,1) = (i+1) * dx
-                   location(num_gp,2) = j * dy + 0.5_wp * dy
-                   location(num_gp,3) = k * dz + 0.5_wp * dz
-                   ei(num_gp)     = e(k+1,j+1,i+1)
-                   dissi(num_gp)  = diss(k+1,j+1,i+1)
-                   de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
-                   de_dyi(num_gp) = 0.0_wp
-                   de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
-                ENDIF
-
-!
-!--             If wall between gridpoint 4 and gridpoint 8, then
-!--             Neumann boundary condition has to be applied
-                IF ( gp_outside_of_building(4) == 1  .AND. &
-                     gp_outside_of_building(8) == 0 )  THEN
-                   num_gp = num_gp + 1
-                   location(num_gp,1) = i * dx + 0.5_wp * dx
-                   location(num_gp,2) = (j+1) * dy
-                   location(num_gp,3) = k * dz + 0.5_wp * dz
-                   ei(num_gp)     = e(k+1,j+1,i)
-                   dissi(num_gp)  = diss(k+1,j+1,i)
-                   de_dxi(num_gp) = 0.0_wp
-                   de_dyi(num_gp) = de_dy(k+1,j+1,i)
-                   de_dzi(num_gp) = de_dz(k+1,j+1,i)
-                ENDIF
-
-                IF ( gp_outside_of_building(8) == 1  .AND. &
-                     gp_outside_of_building(4) == 0 )  THEN
-                   num_gp = num_gp + 1
-                   location(num_gp,1) = i * dx + 0.5_wp * dx
-                   location(num_gp,2) = (j+1) * dy
-                   location(num_gp,3) = k * dz + 0.5_wp * dz
-                   ei(num_gp)     = e(k+1,j+1,i+1)
-                   dissi(num_gp)  = diss(k+1,j+1,i+1)
-                   de_dxi(num_gp) = 0.0_wp
-                   de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
-                   de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
-                ENDIF
-
-!
-!--             If wall between gridpoint 3 and gridpoint 4, then
-!--             Neumann boundary condition has to be applied
-                IF ( gp_outside_of_building(3) == 1  .AND. &
-                     gp_outside_of_building(4) == 0 )  THEN
-                   num_gp = num_gp + 1
-                   location(num_gp,1) = i * dx
-                   location(num_gp,2) = j * dy + 0.5_wp * dy
-                   location(num_gp,3) = k * dz + 0.5_wp * dz
-                   ei(num_gp)     = e(k+1,j,i)
-                   dissi(num_gp)  = diss(k+1,j,i)
-                   de_dxi(num_gp) = de_dx(k+1,j,i)
-                   de_dyi(num_gp) = 0.0_wp
-                   de_dzi(num_gp) = de_dz(k+1,j,i)
-                ENDIF
-
-                IF ( gp_outside_of_building(4) == 1  .AND. &
-                     gp_outside_of_building(3) == 0 )  THEN
-                   num_gp = num_gp + 1
-                   location(num_gp,1) = i * dx
-                   location(num_gp,2) = j * dy + 0.5_wp * dy
-                   location(num_gp,3) = k * dz + 0.5_wp * dz
-                   ei(num_gp)     = e(k+1,j+1,i)
-                   dissi(num_gp)  = diss(k+1,j+1,i)
-                   de_dxi(num_gp) = de_dx(k+1,j+1,i)
-                   de_dyi(num_gp) = 0.0_wp
-                   de_dzi(num_gp) = de_dz(k+1,j+1,i)
-                ENDIF
-
-!
-!--             If wall between gridpoint 1 and gridpoint 3, then
-!--             Neumann boundary condition has to be applied
-!--             (only one case as only building beneath is possible)
-                IF ( gp_outside_of_building(1) == 0  .AND. &
-                     gp_outside_of_building(3) == 1 )  THEN
-                   num_gp = num_gp + 1
-                   location(num_gp,1) = i * dx
-                   location(num_gp,2) = j * dy
-                   location(num_gp,3) = k * dz
-                   ei(num_gp)     = e(k+1,j,i)
-                   dissi(num_gp)  = diss(k+1,j,i)
-                   de_dxi(num_gp) = de_dx(k+1,j,i)
-                   de_dyi(num_gp) = de_dy(k+1,j,i)
-                   de_dzi(num_gp) = 0.0_wp
-                ENDIF
-
-!
-!--             If wall between gridpoint 5 and gridpoint 7, then
-!--             Neumann boundary condition has to be applied
-!--             (only one case as only building beneath is possible)
-                IF ( gp_outside_of_building(5) == 0  .AND. &
-                     gp_outside_of_building(7) == 1 )  THEN
-                   num_gp = num_gp + 1
+                   gp_outside_of_building(7) = 1
                    location(num_gp,1) = (i+1) * dx
                    location(num_gp,2) = j * dy
-                   location(num_gp,3) = k * dz
+                   location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
                    ei(num_gp)     = e(k+1,j,i+1)
                    dissi(num_gp)  = diss(k+1,j,i+1)
                    de_dxi(num_gp) = de_dx(k+1,j,i+1)
                    de_dyi(num_gp) = de_dy(k+1,j,i+1)
-                   de_dzi(num_gp) = 0.0_wp
-                ENDIF
-
-!
-!--             If wall between gridpoint 2 and gridpoint 4, then
-!--             Neumann boundary condition has to be applied
-!--             (only one case as only building beneath is possible)
-                IF ( gp_outside_of_building(2) == 0  .AND. &
-                     gp_outside_of_building(4) == 1 )  THEN
+                   de_dzi(num_gp) = de_dz(k+1,j,i+1)
+                ENDIF
+
+!
+!--             Determine vertical index of topography top at (j+1,i+1)
+                k_wall = get_topography_top_index( j+1, i+1, 's' )
+
+                IF ( k+1 > k_wall  .OR.  k_wall == 0)  THEN
                    num_gp = num_gp + 1
-                   location(num_gp,1) = i * dx
-                   location(num_gp,2) = (j+1) * dy
-                   location(num_gp,3) = k * dz
-                   ei(num_gp)     = e(k+1,j+1,i)
-                   dissi(num_gp)  = diss(k+1,j+1,i)
-                   de_dxi(num_gp) = de_dx(k+1,j+1,i)
-                   de_dyi(num_gp) = de_dy(k+1,j+1,i)
-                   de_dzi(num_gp) = 0.0_wp
-                ENDIF
-
-! 
-!--             If wall between gridpoint 6 and gridpoint 8, then
-!--             Neumann boundary condition has to be applied
-!--             (only one case as only building beneath is possible)
-                IF ( gp_outside_of_building(6) == 0  .AND. &
-                     gp_outside_of_building(8) == 1 )  THEN
-                   num_gp = num_gp + 1
+                   gp_outside_of_building(8) = 1
                    location(num_gp,1) = (i+1) * dx
                    location(num_gp,2) = (j+1) * dy
-                   location(num_gp,3) = k * dz
+                   location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
                    ei(num_gp)     = e(k+1,j+1,i+1)
                    dissi(num_gp)  = diss(k+1,j+1,i+1)
                    de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
                    de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
-                   de_dzi(num_gp) = 0.0_wp
-                ENDIF
-  
-!
-!--             Carry out the interpolation
-                IF ( num_gp == 1 )  THEN
+                   de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
+                ENDIF
+!
+!--             If all gridpoints are situated outside of a building, then the
+!--             ordinary interpolation scheme can be used.
+                IF ( num_gp == 8 )  THEN
+
+                   x  = particles(n)%x - i * dx
+                   y  = particles(n)%y - j * dy
+                   aa = x**2          + y**2
+                   bb = ( dx - x )**2 + y**2
+                   cc = x**2          + ( dy - y )**2
+                   dd = ( dx - x )**2 + ( dy - y )**2
+                   gg = aa + bb + cc + dd
+     
+                   e_int_l = ( ( gg - aa ) * e(k,j,i)   + ( gg - bb ) * e(k,j,i+1)   &
+                             + ( gg - cc ) * e(k,j+1,i) + ( gg - dd ) * e(k,j+1,i+1) &
+                             ) / ( 3.0_wp * gg )
+      
+                   IF ( k == nzt )  THEN
+                      e_int(n) = e_int_l
+                   ELSE
+                      e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
+                                  ( gg - bb ) * e(k+1,j,i+1) + &
+                                  ( gg - cc ) * e(k+1,j+1,i) + &
+                                  ( gg - dd ) * e(k+1,j+1,i+1) &
+                                ) / ( 3.0_wp * gg )
+                      e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
+                                          ( e_int_u - e_int_l )
+                   ENDIF
 ! 
-!--                If only one of the gridpoints is situated outside of the
-!--                building, it follows that the values at the particle
-!--                location are the same as the gridpoint values
-                   e_int(n)     = ei(num_gp)    
-                   diss_int(n)  = dissi(num_gp) 
-                   de_dx_int(n) = de_dxi(num_gp)
-                   de_dy_int(n) = de_dyi(num_gp)
-                   de_dz_int(n) = de_dzi(num_gp)
+!--                Needed to avoid NaN particle velocities (this might not be
+!--                required any more)
+                   IF ( e_int(n) <= 0.0_wp )  THEN
+                      e_int(n) = 1.0E-20_wp
+                   ENDIF
+!
+!--                Interpolate the TKE gradient along x (adopt incides i,j,k
+!--                and all position variables from above (TKE))
+                   de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
+                                   ( gg - bb ) * de_dx(k,j,i+1) + &
+                                   ( gg - cc ) * de_dx(k,j+1,i) + &
+                                   ( gg - dd ) * de_dx(k,j+1,i+1) &
+                                 ) / ( 3.0_wp * gg )
+
+                   IF ( ( k == nzt )  .OR.  ( k == nzb ) )  THEN
+                      de_dx_int(n) = de_dx_int_l
+                   ELSE
+                      de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
+                                      ( gg - bb ) * de_dx(k+1,j,i+1) + &
+                                      ( gg - cc ) * de_dx(k+1,j+1,i) + &
+                                      ( gg - dd ) * de_dx(k+1,j+1,i+1) &
+                                    ) / ( 3.0_wp * gg )
+                      de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / &
+                                              dz * ( de_dx_int_u - de_dx_int_l )
+                   ENDIF
+
+!
+!--                Interpolate the TKE gradient along y
+                   de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
+                                   ( gg - bb ) * de_dy(k,j,i+1) + &
+                                   ( gg - cc ) * de_dy(k,j+1,i) + &
+                                   ( gg - dd ) * de_dy(k,j+1,i+1) &
+                                 ) / ( 3.0_wp * gg )
+                   IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
+                      de_dy_int(n) = de_dy_int_l
+                   ELSE
+                      de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
+                                      ( gg - bb ) * de_dy(k+1,j,i+1) + &
+                                      ( gg - cc ) * de_dy(k+1,j+1,i) + &
+                                      ( gg - dd ) * de_dy(k+1,j+1,i+1) &
+                                    ) / ( 3.0_wp * gg )
+                      de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / &
+                                              dz * ( de_dy_int_u - de_dy_int_l )
+                   ENDIF
+
+!
+!--                Interpolate the TKE gradient along z
+                   IF ( zv(n) < 0.5_wp * dz )  THEN
+                      de_dz_int(n) = 0.0_wp
+                   ELSE
+                      de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
+                                      ( gg - bb ) * de_dz(k,j,i+1) + &
+                                      ( gg - cc ) * de_dz(k,j+1,i) + &
+                                      ( gg - dd ) * de_dz(k,j+1,i+1) &
+                                    ) / ( 3.0_wp * gg )
+    
+                      IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
+                         de_dz_int(n) = de_dz_int_l
+                      ELSE
+                         de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
+                                         ( gg - bb ) * de_dz(k+1,j,i+1) + &
+                                         ( gg - cc ) * de_dz(k+1,j+1,i) + &
+                                         ( gg - dd ) * de_dz(k+1,j+1,i+1) &
+                                       ) / ( 3.0_wp * gg )
+                         de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) /&
+                                              dz * ( de_dz_int_u - de_dz_int_l )
+                      ENDIF
+                   ENDIF
+
+!
+!--                Interpolate the dissipation of TKE
+                   diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
+                                  ( gg - bb ) * diss(k,j,i+1) + &
+                                  ( gg - cc ) * diss(k,j+1,i) + &
+                                  ( gg - dd ) * diss(k,j+1,i+1) &
+                                ) / ( 3.0_wp * gg )
+    
+                   IF ( k == nzt )  THEN
+                      diss_int(n) = diss_int_l
+                   ELSE
+                      diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
+                                     ( gg - bb ) * diss(k+1,j,i+1) + &
+                                     ( gg - cc ) * diss(k+1,j+1,i) + &
+                                     ( gg - dd ) * diss(k+1,j+1,i+1) &
+                                   ) / ( 3.0_wp * gg )
+                      diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz *&
+                                              ( diss_int_u - diss_int_l )
+                   ENDIF
+!
+!--                Set flag for stochastic equation.
+                   term_1_2(n) = 1.0_wp
+    
+                ELSE
+    
+!
+!--                If wall between gridpoint 1 and gridpoint 5, then
+!--                Neumann boundary condition has to be applied
+                   IF ( gp_outside_of_building(1) == 1  .AND. &
+                        gp_outside_of_building(5) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx + 0.5_wp * dx
+                      location(num_gp,2) = j * dy
+                      location(num_gp,3) = k * dz - 0.5_wp * dz
+                      ei(num_gp)     = e(k,j,i)
+                      dissi(num_gp)  = diss(k,j,i)
+                      de_dxi(num_gp) = 0.0_wp
+                      de_dyi(num_gp) = de_dy(k,j,i)
+                      de_dzi(num_gp) = de_dz(k,j,i)
+                   ENDIF
+
+                   IF ( gp_outside_of_building(5) == 1  .AND. &
+                        gp_outside_of_building(1) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx + 0.5_wp * dx
+                      location(num_gp,2) = j * dy
+                      location(num_gp,3) = k * dz - 0.5_wp * dz
+                      ei(num_gp)     = e(k,j,i+1)
+                      dissi(num_gp)  = diss(k,j,i+1)
+                      de_dxi(num_gp) = 0.0_wp
+                      de_dyi(num_gp) = de_dy(k,j,i+1)
+                      de_dzi(num_gp) = de_dz(k,j,i+1)
+                   ENDIF
+
+!
+!--                If wall between gridpoint 5 and gridpoint 6, then
+!--                then Neumann boundary condition has to be applied
+                   IF ( gp_outside_of_building(5) == 1  .AND. &
+                        gp_outside_of_building(6) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = (i+1) * dx
+                      location(num_gp,2) = j * dy + 0.5_wp * dy
+                      location(num_gp,3) = k * dz - 0.5_wp * dz
+                      ei(num_gp)     = e(k,j,i+1)
+                      dissi(num_gp)  = diss(k,j,i+1)
+                      de_dxi(num_gp) = de_dx(k,j,i+1)
+                      de_dyi(num_gp) = 0.0_wp
+                      de_dzi(num_gp) = de_dz(k,j,i+1)
+                   ENDIF
+
+                   IF ( gp_outside_of_building(6) == 1  .AND. &
+                        gp_outside_of_building(5) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = (i+1) * dx
+                      location(num_gp,2) = j * dy + 0.5_wp * dy
+                      location(num_gp,3) = k * dz - 0.5_wp * dz
+                      ei(num_gp)     = e(k,j+1,i+1)
+                      dissi(num_gp)  = diss(k,j+1,i+1)
+                      de_dxi(num_gp) = de_dx(k,j+1,i+1)
+                      de_dyi(num_gp) = 0.0_wp
+                      de_dzi(num_gp) = de_dz(k,j+1,i+1)
+                   ENDIF
+
+!
+!--                If wall between gridpoint 2 and gridpoint 6, then
+!--                Neumann boundary condition has to be applied
+                   IF ( gp_outside_of_building(2) == 1  .AND. &
+                        gp_outside_of_building(6) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx + 0.5_wp * dx
+                      location(num_gp,2) = (j+1) * dy
+                      location(num_gp,3) = k * dz - 0.5_wp * dz
+                      ei(num_gp)     = e(k,j+1,i)
+                      dissi(num_gp)  = diss(k,j+1,i)
+                      de_dxi(num_gp) = 0.0_wp
+                      de_dyi(num_gp) = de_dy(k,j+1,i)
+                      de_dzi(num_gp) = de_dz(k,j+1,i)
+                   ENDIF
+
+                   IF ( gp_outside_of_building(6) == 1  .AND. &
+                        gp_outside_of_building(2) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx + 0.5_wp * dx
+                      location(num_gp,2) = (j+1) * dy
+                      location(num_gp,3) = k * dz - 0.5_wp * dz
+                      ei(num_gp)     = e(k,j+1,i+1)
+                      dissi(num_gp)  = diss(k,j+1,i+1)
+                      de_dxi(num_gp) = 0.0_wp
+                      de_dyi(num_gp) = de_dy(k,j+1,i+1)
+                      de_dzi(num_gp) = de_dz(k,j+1,i+1)
+                   ENDIF
+
+!
+!--                If wall between gridpoint 1 and gridpoint 2, then
+!--                Neumann boundary condition has to be applied
+                   IF ( gp_outside_of_building(1) == 1  .AND. &
+                        gp_outside_of_building(2) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx
+                      location(num_gp,2) = j * dy + 0.5_wp * dy
+                      location(num_gp,3) = k * dz - 0.5_wp * dz
+                      ei(num_gp)     = e(k,j,i)
+                      dissi(num_gp)  = diss(k,j,i)
+                      de_dxi(num_gp) = de_dx(k,j,i)
+                      de_dyi(num_gp) = 0.0_wp
+                      de_dzi(num_gp) = de_dz(k,j,i)  
+                   ENDIF
+
+                   IF ( gp_outside_of_building(2) == 1  .AND. &
+                        gp_outside_of_building(1) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx
+                      location(num_gp,2) = j * dy + 0.5_wp * dy
+                      location(num_gp,3) = k * dz - 0.5_wp * dz
+                      ei(num_gp)     = e(k,j+1,i)
+                      dissi(num_gp)  = diss(k,j+1,i)
+                      de_dxi(num_gp) = de_dx(k,j+1,i)
+                      de_dyi(num_gp) = 0.0_wp
+                      de_dzi(num_gp) = de_dz(k,j+1,i)
+                   ENDIF
+
+!
+!--                If wall between gridpoint 3 and gridpoint 7, then
+!--                Neumann boundary condition has to be applied
+                   IF ( gp_outside_of_building(3) == 1  .AND. &
+                        gp_outside_of_building(7) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx + 0.5_wp * dx
+                      location(num_gp,2) = j * dy
+                      location(num_gp,3) = k * dz + 0.5_wp * dz
+                      ei(num_gp)     = e(k+1,j,i)
+                      dissi(num_gp)  = diss(k+1,j,i)
+                      de_dxi(num_gp) = 0.0_wp
+                      de_dyi(num_gp) = de_dy(k+1,j,i)
+                      de_dzi(num_gp) = de_dz(k+1,j,i) 
+                   ENDIF
+
+                   IF ( gp_outside_of_building(7) == 1  .AND. &
+                        gp_outside_of_building(3) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx + 0.5_wp * dx
+                      location(num_gp,2) = j * dy
+                      location(num_gp,3) = k * dz + 0.5_wp * dz
+                      ei(num_gp)     = e(k+1,j,i+1)
+                      dissi(num_gp)  = diss(k+1,j,i+1)
+                      de_dxi(num_gp) = 0.0_wp
+                      de_dyi(num_gp) = de_dy(k+1,j,i+1)
+                      de_dzi(num_gp) = de_dz(k+1,j,i+1)
+                   ENDIF
+
+!
+!--                If wall between gridpoint 7 and gridpoint 8, then
+!--                Neumann boundary condition has to be applied
+                   IF ( gp_outside_of_building(7) == 1  .AND. &
+                        gp_outside_of_building(8) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = (i+1) * dx
+                      location(num_gp,2) = j * dy + 0.5_wp * dy
+                      location(num_gp,3) = k * dz + 0.5_wp * dz
+                      ei(num_gp)     = e(k+1,j,i+1)
+                      dissi(num_gp)  = diss(k+1,j,i+1)
+                      de_dxi(num_gp) = de_dx(k+1,j,i+1)
+                      de_dyi(num_gp) = 0.0_wp
+                      de_dzi(num_gp) = de_dz(k+1,j,i+1)
+                   ENDIF
+
+                   IF ( gp_outside_of_building(8) == 1  .AND. &
+                        gp_outside_of_building(7) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = (i+1) * dx
+                      location(num_gp,2) = j * dy + 0.5_wp * dy
+                      location(num_gp,3) = k * dz + 0.5_wp * dz
+                      ei(num_gp)     = e(k+1,j+1,i+1)
+                      dissi(num_gp)  = diss(k+1,j+1,i+1)
+                      de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
+                      de_dyi(num_gp) = 0.0_wp
+                      de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
+                   ENDIF
+
+!
+!--                If wall between gridpoint 4 and gridpoint 8, then
+!--                Neumann boundary condition has to be applied
+                   IF ( gp_outside_of_building(4) == 1  .AND. &
+                        gp_outside_of_building(8) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx + 0.5_wp * dx
+                      location(num_gp,2) = (j+1) * dy
+                      location(num_gp,3) = k * dz + 0.5_wp * dz
+                      ei(num_gp)     = e(k+1,j+1,i)
+                      dissi(num_gp)  = diss(k+1,j+1,i)
+                      de_dxi(num_gp) = 0.0_wp
+                      de_dyi(num_gp) = de_dy(k+1,j+1,i)
+                      de_dzi(num_gp) = de_dz(k+1,j+1,i)
+                   ENDIF
+
+                   IF ( gp_outside_of_building(8) == 1  .AND. &
+                        gp_outside_of_building(4) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx + 0.5_wp * dx
+                      location(num_gp,2) = (j+1) * dy
+                      location(num_gp,3) = k * dz + 0.5_wp * dz
+                      ei(num_gp)     = e(k+1,j+1,i+1)
+                      dissi(num_gp)  = diss(k+1,j+1,i+1)
+                      de_dxi(num_gp) = 0.0_wp
+                      de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
+                      de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
+                   ENDIF
+
+!
+!--                If wall between gridpoint 3 and gridpoint 4, then
+!--                Neumann boundary condition has to be applied
+                   IF ( gp_outside_of_building(3) == 1  .AND. &
+                        gp_outside_of_building(4) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx
+                      location(num_gp,2) = j * dy + 0.5_wp * dy
+                      location(num_gp,3) = k * dz + 0.5_wp * dz
+                      ei(num_gp)     = e(k+1,j,i)
+                      dissi(num_gp)  = diss(k+1,j,i)
+                      de_dxi(num_gp) = de_dx(k+1,j,i)
+                      de_dyi(num_gp) = 0.0_wp
+                      de_dzi(num_gp) = de_dz(k+1,j,i)
+                   ENDIF
+
+                   IF ( gp_outside_of_building(4) == 1  .AND. &
+                        gp_outside_of_building(3) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx
+                      location(num_gp,2) = j * dy + 0.5_wp * dy
+                      location(num_gp,3) = k * dz + 0.5_wp * dz
+                      ei(num_gp)     = e(k+1,j+1,i)
+                      dissi(num_gp)  = diss(k+1,j+1,i)
+                      de_dxi(num_gp) = de_dx(k+1,j+1,i)
+                      de_dyi(num_gp) = 0.0_wp
+                      de_dzi(num_gp) = de_dz(k+1,j+1,i)
+                   ENDIF
+
+!
+!--                If wall between gridpoint 1 and gridpoint 3, then
+!--                Neumann boundary condition has to be applied
+!--                (only one case as only building beneath is possible)
+                   IF ( gp_outside_of_building(1) == 0  .AND. &
+                        gp_outside_of_building(3) == 1 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx
+                      location(num_gp,2) = j * dy
+                      location(num_gp,3) = k * dz
+                      ei(num_gp)     = e(k+1,j,i)
+                      dissi(num_gp)  = diss(k+1,j,i)
+                      de_dxi(num_gp) = de_dx(k+1,j,i)
+                      de_dyi(num_gp) = de_dy(k+1,j,i)
+                      de_dzi(num_gp) = 0.0_wp
+                   ENDIF
+
+!
+!--                If wall between gridpoint 5 and gridpoint 7, then
+!--                Neumann boundary condition has to be applied
+!--                (only one case as only building beneath is possible)
+                   IF ( gp_outside_of_building(5) == 0  .AND. &
+                        gp_outside_of_building(7) == 1 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = (i+1) * dx
+                      location(num_gp,2) = j * dy
+                      location(num_gp,3) = k * dz
+                      ei(num_gp)     = e(k+1,j,i+1)
+                      dissi(num_gp)  = diss(k+1,j,i+1)
+                      de_dxi(num_gp) = de_dx(k+1,j,i+1)
+                      de_dyi(num_gp) = de_dy(k+1,j,i+1)
+                      de_dzi(num_gp) = 0.0_wp
+                   ENDIF
+
+!
+!--                If wall between gridpoint 2 and gridpoint 4, then
+!--                Neumann boundary condition has to be applied
+!--                (only one case as only building beneath is possible)
+                   IF ( gp_outside_of_building(2) == 0  .AND. &
+                        gp_outside_of_building(4) == 1 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx
+                      location(num_gp,2) = (j+1) * dy
+                      location(num_gp,3) = k * dz
+                      ei(num_gp)     = e(k+1,j+1,i)
+                      dissi(num_gp)  = diss(k+1,j+1,i)
+                      de_dxi(num_gp) = de_dx(k+1,j+1,i)
+                      de_dyi(num_gp) = de_dy(k+1,j+1,i)
+                      de_dzi(num_gp) = 0.0_wp
+                   ENDIF
+
+! 
+!--                If wall between gridpoint 6 and gridpoint 8, then
+!--                Neumann boundary condition has to be applied
+!--                (only one case as only building beneath is possible)
+                   IF ( gp_outside_of_building(6) == 0  .AND. &
+                        gp_outside_of_building(8) == 1 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = (i+1) * dx
+                      location(num_gp,2) = (j+1) * dy
+                      location(num_gp,3) = k * dz
+                      ei(num_gp)     = e(k+1,j+1,i+1)
+                      dissi(num_gp)  = diss(k+1,j+1,i+1)
+                      de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
+                      de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
+                      de_dzi(num_gp) = 0.0_wp
+                   ENDIF
+     
+!
+!--                Carry out the interpolation
+                   IF ( num_gp == 1 )  THEN
+! 
+!--                   If only one of the gridpoints is situated outside of the
+!--                   building, it follows that the values at the particle
+!--                   location are the same as the gridpoint values
+                      e_int(n)     = ei(num_gp)    
+                      diss_int(n)  = dissi(num_gp) 
+                      de_dx_int(n) = de_dxi(num_gp)
+                      de_dy_int(n) = de_dyi(num_gp)
+                      de_dz_int(n) = de_dzi(num_gp)
+!
+!--                   Set flag for stochastic equation which disables calculation 
+!--                   of drift and memory term near topography.
+                      term_1_2(n) = 0.0_wp
+                   ELSE IF ( num_gp > 1 )  THEN
+    
+                      d_sum = 0.0_wp
+! 
+!--                   Evaluation of the distances between the gridpoints
+!--                   contributing to the interpolated values, and the particle
+!--                   location
+                      DO  agp = 1, num_gp
+                         d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2  &
+                                      + ( particles(n)%y-location(agp,2) )**2  &
+                                      + ( zv(n)-location(agp,3) )**2
+                         d_sum        = d_sum + d_gp_pl(agp)
+                      ENDDO
+    
+!
+!--                  Finally the interpolation can be carried out
+                      e_int(n)     = 0.0_wp 
+                      diss_int(n)  = 0.0_wp 
+                      de_dx_int(n) = 0.0_wp  
+                      de_dy_int(n) = 0.0_wp  
+                      de_dz_int(n) = 0.0_wp  
+                      DO  agp = 1, num_gp
+                         e_int(n)     = e_int(n)     + ( d_sum - d_gp_pl(agp) ) * &
+                                                ei(agp) / ( (num_gp-1) * d_sum )
+                         diss_int(n)  = diss_int(n)  + ( d_sum - d_gp_pl(agp) ) * &
+                                             dissi(agp) / ( (num_gp-1) * d_sum )
+                         de_dx_int(n) = de_dx_int(n) + ( d_sum - d_gp_pl(agp) ) * &
+                                            de_dxi(agp) / ( (num_gp-1) * d_sum )
+                         de_dy_int(n) = de_dy_int(n) + ( d_sum - d_gp_pl(agp) ) * &
+                                            de_dyi(agp) / ( (num_gp-1) * d_sum )
+                         de_dz_int(n) = de_dz_int(n) + ( d_sum - d_gp_pl(agp) ) * &
+                                            de_dzi(agp) / ( (num_gp-1) * d_sum )
+                      ENDDO
+    
+                   ENDIF
+                   e_int(n)     = MAX( 1E-20_wp, e_int(n)     )
+                   diss_int(n)  = MAX( 1E-20_wp, diss_int(n)  )
+                   de_dx_int(n) = MAX( 1E-20_wp, de_dx_int(n) )
+                   de_dy_int(n) = MAX( 1E-20_wp, de_dy_int(n) ) 
+                   de_dz_int(n) = MAX( 1E-20_wp, de_dz_int(n) ) 
 !
 !--                Set flag for stochastic equation which disables calculation 
 !--                of drift and memory term near topography.
                    term_1_2(n) = 0.0_wp
-                ELSE IF ( num_gp > 1 )  THEN
- 
-                   d_sum = 0.0_wp
-! 
-!--                Evaluation of the distances between the gridpoints
-!--                contributing to the interpolated values, and the particle
-!--                location
-                   DO  agp = 1, num_gp
-                      d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2  &
-                                   + ( particles(n)%y-location(agp,2) )**2  &
-                                   + ( zv(n)-location(agp,3) )**2
-                      d_sum        = d_sum + d_gp_pl(agp)
-                   ENDDO
- 
-!
-!--                Finally the interpolation can be carried out
-                   e_int(n)     = 0.0_wp 
-                   diss_int(n)  = 0.0_wp 
-                   de_dx_int(n) = 0.0_wp  
-                   de_dy_int(n) = 0.0_wp  
-                   de_dz_int(n) = 0.0_wp  
-                   DO  agp = 1, num_gp
-                      e_int(n)     = e_int(n)     + ( d_sum - d_gp_pl(agp) ) * &
-                                             ei(agp) / ( (num_gp-1) * d_sum )
-                      diss_int(n)  = diss_int(n)  + ( d_sum - d_gp_pl(agp) ) * &
-                                          dissi(agp) / ( (num_gp-1) * d_sum )
-                      de_dx_int(n) = de_dx_int(n) + ( d_sum - d_gp_pl(agp) ) * &
-                                         de_dxi(agp) / ( (num_gp-1) * d_sum )
-                      de_dy_int(n) = de_dy_int(n) + ( d_sum - d_gp_pl(agp) ) * &
-                                         de_dyi(agp) / ( (num_gp-1) * d_sum )
-                      de_dz_int(n) = de_dz_int(n) + ( d_sum - d_gp_pl(agp) ) * &
-                                         de_dzi(agp) / ( (num_gp-1) * d_sum )
-                   ENDDO
- 
-                ENDIF
-                e_int(n)     = MAX( 1E-20_wp, e_int(n)     )
-                diss_int(n)  = MAX( 1E-20_wp, diss_int(n)  )
-                de_dx_int(n) = MAX( 1E-20_wp, de_dx_int(n) )
-                de_dy_int(n) = MAX( 1E-20_wp, de_dy_int(n) ) 
-                de_dz_int(n) = MAX( 1E-20_wp, de_dz_int(n) ) 
-!
-!--             Set flag for stochastic equation which disables calculation 
-!--             of drift and memory term near topography.
-                term_1_2(n) = 0.0_wp
-             ENDIF
+                ENDIF
+             ENDDO
           ENDDO
        ENDIF
@@ -1345,94 +1355,103 @@
        ENDDO
 
-       DO  n = 1, number_of_particles
-
-         rg(n,1) = random_gauss( iran_part, 5.0_wp )
-         rg(n,2) = random_gauss( iran_part, 5.0_wp )
-         rg(n,3) = random_gauss( iran_part, 5.0_wp )
-
+       DO  nb = 0, 7
+          DO  n = start_index(nb), end_index(nb)
+             rg(n,1) = random_gauss( iran_part, 5.0_wp )
+             rg(n,2) = random_gauss( iran_part, 5.0_wp )
+             rg(n,3) = random_gauss( iran_part, 5.0_wp )
+          ENDDO
        ENDDO
 
-       DO  n = 1, number_of_particles
-!
-!--       Calculate the Lagrangian timescale according to Weil et al. (2004).
-          lagr_timescale = ( 4.0_wp * e_int(n) + 1E-20_wp ) / &
-                           ( 3.0_wp * fs_int(n) * c_0 * diss_int(n) + 1E-20_wp )
-
-!
-!--       Calculate the next particle timestep. dt_gap is the time needed to
-!--       complete the current LES timestep.
-          dt_gap = dt_3d - particles(n)%dt_sum
-          dt_particle(n) = MIN( dt_3d, 0.025_wp * lagr_timescale, dt_gap )
-
-!
-!--       The particle timestep should not be too small in order to prevent
-!--       the number of particle timesteps of getting too large
-          IF ( dt_particle(n) < dt_min_part  .AND.  dt_min_part < dt_gap )  THEN
-             dt_particle(n) = dt_min_part
-          ENDIF
-
-!
-!--       Calculate the SGS velocity components
-          IF ( particles(n)%age == 0.0_wp )  THEN
-!
-!--          For new particles the SGS components are derived from the SGS
-!--          TKE. Limit the Gaussian random number to the interval
-!--          [-5.0*sigma, 5.0*sigma] in order to prevent the SGS velocities
-!--          from becoming unrealistically large.
-             particles(n)%rvar1 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) *   &
-                                  ( rg(n,1) - 1.0_wp )
-             particles(n)%rvar2 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) *   &
-                                  ( rg(n,2) - 1.0_wp )
-             particles(n)%rvar3 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) *   &
-                                  ( rg(n,3) - 1.0_wp )
-
-          ELSE
-!
-!--          Restriction of the size of the new timestep: compared to the 
-!--          previous timestep the increase must not exceed 200%
-
-             dt_particle_m = particles(n)%age - particles(n)%age_m
-             IF ( dt_particle(n) > 2.0_wp * dt_particle_m )  THEN
-                dt_particle(n) = 2.0_wp * dt_particle_m
+       DO  nb = 0, 7
+          DO  n = start_index(nb), end_index(nb)
+
+!
+!--          Calculate the Lagrangian timescale according to Weil et al. (2004).
+             lagr_timescale = ( 4.0_wp * e_int(n) + 1E-20_wp ) / &
+                              ( 3.0_wp * fs_int(n) * c_0 * diss_int(n) + 1E-20_wp )
+
+!
+!--          Calculate the next particle timestep. dt_gap is the time needed to
+!--          complete the current LES timestep.
+             dt_gap = dt_3d - particles(n)%dt_sum
+             dt_particle(n) = MIN( dt_3d, 0.025_wp * lagr_timescale, dt_gap )
+
+!
+!--          The particle timestep should not be too small in order to prevent
+!--          the number of particle timesteps of getting too large
+             IF ( dt_particle(n) < dt_min_part  .AND.  dt_min_part < dt_gap )  THEN
+                dt_particle(n) = dt_min_part
              ENDIF
 
 !
-!--          For old particles the SGS components are correlated with the
-!--          values from the previous timestep. Random numbers have also to
-!--          be limited (see above).
-!--          As negative values for the subgrid TKE are not allowed, the
-!--          change of the subgrid TKE with time cannot be smaller than
-!--          -e_int(n)/dt_particle. This value is used as a lower boundary
-!--          value for the change of TKE 
-
-             de_dt_min = - e_int(n) / dt_particle(n)
-
-             de_dt = ( e_int(n) - particles(n)%e_m ) / dt_particle_m
-
-             IF ( de_dt < de_dt_min )  THEN
-                de_dt = de_dt_min
+!--          Calculate the SGS velocity components
+             IF ( particles(n)%age == 0.0_wp )  THEN
+!
+!--             For new particles the SGS components are derived from the SGS
+!--             TKE. Limit the Gaussian random number to the interval
+!--             [-5.0*sigma, 5.0*sigma] in order to prevent the SGS velocities
+!--             from becoming unrealistically large.
+                particles(n)%rvar1 = SQRT( 2.0_wp * sgs_wf_part * e_int(n)     &
+                                          + 1E-20_wp ) *                       &
+                                     ( rg(n,1) - 1.0_wp )
+                particles(n)%rvar2 = SQRT( 2.0_wp * sgs_wf_part * e_int(n)     &
+                                          + 1E-20_wp ) *                       &
+                                     ( rg(n,2) - 1.0_wp )
+                particles(n)%rvar3 = SQRT( 2.0_wp * sgs_wf_part * e_int(n)     &
+                                          + 1E-20_wp ) *                       &
+                                     ( rg(n,3) - 1.0_wp )
+
+             ELSE
+!
+!--             Restriction of the size of the new timestep: compared to the 
+!--             previous timestep the increase must not exceed 200%. First,
+!--             check if age > age_m, in order to prevent that particles get zero
+!--             timestep.
+                dt_particle_m = MERGE( dt_particle(n),                         &
+                                       particles(n)%age - particles(n)%age_m,  &
+                                       particles(n)%age - particles(n)%age_m < &
+                                       1E-8_wp )
+                IF ( dt_particle(n) > 2.0_wp * dt_particle_m )  THEN
+                   dt_particle(n) = 2.0_wp * dt_particle_m
+                ENDIF
+
+!
+!--             For old particles the SGS components are correlated with the
+!--             values from the previous timestep. Random numbers have also to
+!--             be limited (see above).
+!--             As negative values for the subgrid TKE are not allowed, the
+!--             change of the subgrid TKE with time cannot be smaller than
+!--             -e_int(n)/dt_particle. This value is used as a lower boundary
+!--             value for the change of TKE 
+                de_dt_min = - e_int(n) / dt_particle(n)
+
+                de_dt = ( e_int(n) - particles(n)%e_m ) / dt_particle_m
+
+                IF ( de_dt < de_dt_min )  THEN
+                   de_dt = de_dt_min
+                ENDIF
+
+                CALL weil_stochastic_eq(particles(n)%rvar1, fs_int(n), e_int(n),& 
+                                        de_dx_int(n), de_dt, diss_int(n),       &
+                                        dt_particle(n), rg(n,1), term_1_2(n) )
+
+                CALL weil_stochastic_eq(particles(n)%rvar2, fs_int(n), e_int(n),& 
+                                        de_dy_int(n), de_dt, diss_int(n),       &
+                                        dt_particle(n), rg(n,2), term_1_2(n) )
+
+                CALL weil_stochastic_eq(particles(n)%rvar3, fs_int(n), e_int(n),& 
+                                        de_dz_int(n), de_dt, diss_int(n),       &
+                                        dt_particle(n), rg(n,3), term_1_2(n) )
+
              ENDIF
 
-             CALL weil_stochastic_eq(particles(n)%rvar1, fs_int(n), e_int(n),  & 
-                                     de_dx_int(n), de_dt, diss_int(n),         &
-                                     dt_particle(n), rg(n,1), term_1_2(n) )
-
-             CALL weil_stochastic_eq(particles(n)%rvar2, fs_int(n), e_int(n),  & 
-                                     de_dy_int(n), de_dt, diss_int(n),         &
-                                     dt_particle(n), rg(n,2), term_1_2(n) )
-
-             CALL weil_stochastic_eq(particles(n)%rvar3, fs_int(n), e_int(n),  & 
-                                     de_dz_int(n), de_dt, diss_int(n),         &
-                                     dt_particle(n), rg(n,3), term_1_2(n) )
-
-          ENDIF
-
-          u_int(n) = u_int(n) + particles(n)%rvar1
-          v_int(n) = v_int(n) + particles(n)%rvar2
-          w_int(n) = w_int(n) + particles(n)%rvar3
-!
-!--       Store the SGS TKE of the current timelevel which is needed for
-!--       for calculating the SGS particle velocities at the next timestep
-          particles(n)%e_m = e_int(n)
+             u_int(n) = u_int(n) + particles(n)%rvar1
+             v_int(n) = v_int(n) + particles(n)%rvar2
+             w_int(n) = w_int(n) + particles(n)%rvar3
+!
+!--          Store the SGS TKE of the current timelevel which is needed for
+!--          for calculating the SGS particle velocities at the next timestep
+             particles(n)%e_m = e_int(n)
+          ENDDO
        ENDDO
 
@@ -1444,44 +1463,112 @@
 
     ENDIF
-!
-!-- Store the old age of the particle ( needed to prevent that a
-!-- particle crosses several PEs during one timestep, and for the
-!-- evaluation of the subgrid particle velocity fluctuations )
-    particles(1:number_of_particles)%age_m = particles(1:number_of_particles)%age
 
     dens_ratio = particle_groups(particles(1:number_of_particles)%group)%density_ratio
 
     IF ( ANY( dens_ratio == 0.0_wp ) )  THEN
-       DO  n = 1, number_of_particles
-
-!
-!--       Particle advection
-          IF ( dens_ratio(n) == 0.0_wp )  THEN
-!
-!--          Pure passive transport (without particle inertia)
-             particles(n)%x = xv(n) + u_int(n) * dt_particle(n)
-             particles(n)%y = yv(n) + v_int(n) * dt_particle(n)
-             particles(n)%z = zv(n) + w_int(n) * dt_particle(n)
-
-             particles(n)%speed_x = u_int(n)
-             particles(n)%speed_y = v_int(n)
-             particles(n)%speed_z = w_int(n)
-
-          ELSE
+       DO  nb = 0, 7
+          DO  n = start_index(nb), end_index(nb)
+
+!
+!--          Particle advection
+             IF ( dens_ratio(n) == 0.0_wp )  THEN
+!
+!--             Pure passive transport (without particle inertia)
+                particles(n)%x = xv(n) + u_int(n) * dt_particle(n)
+                particles(n)%y = yv(n) + v_int(n) * dt_particle(n)
+                particles(n)%z = zv(n) + w_int(n) * dt_particle(n)
+
+                particles(n)%speed_x = u_int(n)
+                particles(n)%speed_y = v_int(n)
+                particles(n)%speed_z = w_int(n)
+
+             ELSE
+!
+!--             Transport of particles with inertia
+                particles(n)%x = particles(n)%x + particles(n)%speed_x * &
+                                                  dt_particle(n)
+                particles(n)%y = particles(n)%y + particles(n)%speed_y * &
+                                                  dt_particle(n)
+                particles(n)%z = particles(n)%z + particles(n)%speed_z * &
+                                                  dt_particle(n)
+
+!
+!--             Update of the particle velocity
+                IF ( cloud_droplets )  THEN
+!
+!--                Terminal velocity is computed for vertical direction (Rogers et 
+!--                al., 1993, J. Appl. Meteorol.)
+                   diameter = particles(n)%radius * 2000.0_wp !diameter in mm
+                   IF ( diameter <= d0_rog )  THEN
+                      w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) )
+                   ELSE
+                      w_s = a_rog - b_rog * EXP( -c_rog * diameter )
+                   ENDIF
+
+!
+!--                If selected, add random velocities following Soelch and Kaercher
+!--                (2010, Q. J. R. Meteorol. Soc.)
+                   IF ( use_sgs_for_particles )  THEN
+                      lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
+                      RL             = EXP( -1.0_wp * dt_3d / lagr_timescale )
+                      sigma          = SQRT( e(kp,jp,ip) )
+
+                      rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
+                      rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
+                      rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
+
+                      particles(n)%rvar1 = RL * particles(n)%rvar1 +              &
+                                           SQRT( 1.0_wp - RL**2 ) * sigma * rg1
+                      particles(n)%rvar2 = RL * particles(n)%rvar2 +              &
+                                           SQRT( 1.0_wp - RL**2 ) * sigma * rg2
+                      particles(n)%rvar3 = RL * particles(n)%rvar3 +              &
+                                           SQRT( 1.0_wp - RL**2 ) * sigma * rg3
+
+                      particles(n)%speed_x = u_int(n) + particles(n)%rvar1
+                      particles(n)%speed_y = v_int(n) + particles(n)%rvar2
+                      particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s
+                   ELSE
+                      particles(n)%speed_x = u_int(n)
+                      particles(n)%speed_y = v_int(n)
+                      particles(n)%speed_z = w_int(n) - w_s
+                   ENDIF
+
+                ELSE
+
+                   IF ( use_sgs_for_particles )  THEN
+                      exp_arg  = particle_groups(particles(n)%group)%exp_arg
+                      exp_term = EXP( -exp_arg * dt_particle(n) )
+                   ELSE
+                      exp_arg  = particle_groups(particles(n)%group)%exp_arg
+                      exp_term = particle_groups(particles(n)%group)%exp_term
+                   ENDIF
+                   particles(n)%speed_x = particles(n)%speed_x * exp_term +         &
+                                          u_int(n) * ( 1.0_wp - exp_term )
+                   particles(n)%speed_y = particles(n)%speed_y * exp_term +         &
+                                          v_int(n) * ( 1.0_wp - exp_term )
+                   particles(n)%speed_z = particles(n)%speed_z * exp_term +         &
+                                          ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * &
+                                          g / exp_arg ) * ( 1.0_wp - exp_term )
+                ENDIF
+
+             ENDIF
+          ENDDO
+       ENDDO
+    
+    ELSE
+
+       DO  nb = 0, 7
+          DO  n = start_index(nb), end_index(nb)
 !
 !--          Transport of particles with inertia
-             particles(n)%x = particles(n)%x + particles(n)%speed_x * &
-                                               dt_particle(n)
-             particles(n)%y = particles(n)%y + particles(n)%speed_y * &
-                                               dt_particle(n)
-             particles(n)%z = particles(n)%z + particles(n)%speed_z * &
-                                               dt_particle(n)
-
+             particles(n)%x = xv(n) + particles(n)%speed_x * dt_particle(n)
+             particles(n)%y = yv(n) + particles(n)%speed_y * dt_particle(n)
+             particles(n)%z = zv(n) + particles(n)%speed_z * dt_particle(n)
 !
 !--          Update of the particle velocity
              IF ( cloud_droplets )  THEN
 !
-!--             Terminal velocity is computed for vertical direction (Rogers et 
-!--             al., 1993, J. Appl. Meteorol.)
+!--             Terminal velocity is computed for vertical direction (Rogers et al., 
+!--             1993, J. Appl. Meteorol.)
                 diameter = particles(n)%radius * 2000.0_wp !diameter in mm
                 IF ( diameter <= d0_rog )  THEN
@@ -1495,26 +1582,26 @@
 !--             (2010, Q. J. R. Meteorol. Soc.)
                 IF ( use_sgs_for_particles )  THEN
-                   lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
-                   RL             = EXP( -1.0_wp * dt_3d / lagr_timescale )
-                   sigma          = SQRT( e(kp,jp,ip) )
-
-                   rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
-                   rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
-                   rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
-
-                   particles(n)%rvar1 = RL * particles(n)%rvar1 +              &
-                                        SQRT( 1.0_wp - RL**2 ) * sigma * rg1
-                   particles(n)%rvar2 = RL * particles(n)%rvar2 +              &
-                                        SQRT( 1.0_wp - RL**2 ) * sigma * rg2
-                   particles(n)%rvar3 = RL * particles(n)%rvar3 +              &
-                                        SQRT( 1.0_wp - RL**2 ) * sigma * rg3
-
-                   particles(n)%speed_x = u_int(n) + particles(n)%rvar1
-                   particles(n)%speed_y = v_int(n) + particles(n)%rvar2
-                   particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s
+                    lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
+                    RL             = EXP( -1.0_wp * dt_3d / lagr_timescale )
+                    sigma          = SQRT( e(kp,jp,ip) )
+
+                    rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
+                    rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
+                    rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
+
+                    particles(n)%rvar1 = RL * particles(n)%rvar1 +                &
+                                         SQRT( 1.0_wp - RL**2 ) * sigma * rg1
+                    particles(n)%rvar2 = RL * particles(n)%rvar2 +                &
+                                         SQRT( 1.0_wp - RL**2 ) * sigma * rg2
+                    particles(n)%rvar3 = RL * particles(n)%rvar3 +                &
+                                         SQRT( 1.0_wp - RL**2 ) * sigma * rg3
+
+                    particles(n)%speed_x = u_int(n) + particles(n)%rvar1
+                    particles(n)%speed_y = v_int(n) + particles(n)%rvar2
+                    particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s
                 ELSE
-                   particles(n)%speed_x = u_int(n)
-                   particles(n)%speed_y = v_int(n)
-                   particles(n)%speed_z = w_int(n) - w_s
+                    particles(n)%speed_x = u_int(n)
+                    particles(n)%speed_y = v_int(n)
+                    particles(n)%speed_z = w_int(n) - w_s
                 ENDIF
 
@@ -1528,102 +1615,39 @@
                    exp_term = particle_groups(particles(n)%group)%exp_term
                 ENDIF
-                particles(n)%speed_x = particles(n)%speed_x * exp_term +         &
+                particles(n)%speed_x = particles(n)%speed_x * exp_term +             &
                                        u_int(n) * ( 1.0_wp - exp_term )
-                particles(n)%speed_y = particles(n)%speed_y * exp_term +         &
+                particles(n)%speed_y = particles(n)%speed_y * exp_term +             &
                                        v_int(n) * ( 1.0_wp - exp_term )
-                particles(n)%speed_z = particles(n)%speed_z * exp_term +         &
-                                       ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * &
-                                       g / exp_arg ) * ( 1.0_wp - exp_term )
+                particles(n)%speed_z = particles(n)%speed_z * exp_term +             &
+                                       ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * g / &
+                                       exp_arg ) * ( 1.0_wp - exp_term )
              ENDIF
-
+          ENDDO
+       ENDDO
+
+    ENDIF
+
+!
+!-- Store the old age of the particle ( needed to prevent that a
+!-- particle crosses several PEs during one timestep, and for the
+!-- evaluation of the subgrid particle velocity fluctuations )
+    particles(1:number_of_particles)%age_m = particles(1:number_of_particles)%age
+
+    DO  nb = 0, 7
+       DO  n = start_index(nb), end_index(nb)
+!
+!--       Increment the particle age and the total time that the particle
+!--       has advanced within the particle timestep procedure
+          particles(n)%age    = particles(n)%age    + dt_particle(n)
+          particles(n)%dt_sum = particles(n)%dt_sum + dt_particle(n)
+
+!
+!--       Check whether there is still a particle that has not yet completed 
+!--       the total LES timestep
+          IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8_wp )  THEN
+             dt_3d_reached_l = .FALSE.
           ENDIF
 
        ENDDO
-    
-    ELSE
-
-       DO  n = 1, number_of_particles
-
-!--       Transport of particles with inertia
-          particles(n)%x = xv(n) + particles(n)%speed_x * dt_particle(n)
-          particles(n)%y = yv(n) + particles(n)%speed_y * dt_particle(n)
-          particles(n)%z = zv(n) + particles(n)%speed_z * dt_particle(n)
-!
-!--       Update of the particle velocity
-          IF ( cloud_droplets )  THEN
-!
-!--          Terminal velocity is computed for vertical direction (Rogers et al., 
-!--          1993, J. Appl. Meteorol.)
-             diameter = particles(n)%radius * 2000.0_wp !diameter in mm
-             IF ( diameter <= d0_rog )  THEN
-                w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) )
-             ELSE
-                w_s = a_rog - b_rog * EXP( -c_rog * diameter )
-             ENDIF
-
-!
-!--          If selected, add random velocities following Soelch and Kaercher
-!--          (2010, Q. J. R. Meteorol. Soc.)
-             IF ( use_sgs_for_particles )  THEN
-                 lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
-                 RL             = EXP( -1.0_wp * dt_3d / lagr_timescale )
-                 sigma          = SQRT( e(kp,jp,ip) )
-
-                 rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
-                 rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
-                 rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
-
-                 particles(n)%rvar1 = RL * particles(n)%rvar1 +                &
-                                      SQRT( 1.0_wp - RL**2 ) * sigma * rg1
-                 particles(n)%rvar2 = RL * particles(n)%rvar2 +                &
-                                      SQRT( 1.0_wp - RL**2 ) * sigma * rg2
-                 particles(n)%rvar3 = RL * particles(n)%rvar3 +                &
-                                      SQRT( 1.0_wp - RL**2 ) * sigma * rg3
-
-                 particles(n)%speed_x = u_int(n) + particles(n)%rvar1
-                 particles(n)%speed_y = v_int(n) + particles(n)%rvar2
-                 particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s
-             ELSE
-                 particles(n)%speed_x = u_int(n)
-                 particles(n)%speed_y = v_int(n)
-                 particles(n)%speed_z = w_int(n) - w_s
-             ENDIF
-
-          ELSE
-
-             IF ( use_sgs_for_particles )  THEN
-                exp_arg  = particle_groups(particles(n)%group)%exp_arg
-                exp_term = EXP( -exp_arg * dt_particle(n) )
-             ELSE
-                exp_arg  = particle_groups(particles(n)%group)%exp_arg
-                exp_term = particle_groups(particles(n)%group)%exp_term
-             ENDIF
-             particles(n)%speed_x = particles(n)%speed_x * exp_term +             &
-                                    u_int(n) * ( 1.0_wp - exp_term )
-             particles(n)%speed_y = particles(n)%speed_y * exp_term +             &
-                                    v_int(n) * ( 1.0_wp - exp_term )
-             particles(n)%speed_z = particles(n)%speed_z * exp_term +             &
-                                    ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * g / &
-                                    exp_arg ) * ( 1.0_wp - exp_term )
-          ENDIF
-
-       ENDDO
-
-    ENDIF
-
-    DO  n = 1, number_of_particles
-!
-!--    Increment the particle age and the total time that the particle
-!--    has advanced within the particle timestep procedure
-       particles(n)%age    = particles(n)%age    + dt_particle(n)
-       particles(n)%dt_sum = particles(n)%dt_sum + dt_particle(n)
-
-!
-!--    Check whether there is still a particle that has not yet completed 
-!--    the total LES timestep
-       IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8_wp )  THEN
-          dt_3d_reached_l = .FALSE.
-       ENDIF
-
     ENDDO
 
@@ -1693,5 +1717,4 @@
 !-- subgrid-scale velocity component is set to zero, in order to prevent a 
 !-- velocity build-up. 
-
 !-- This case, set also previous subgrid-scale component to zero. 
     v_sgs = v_sgs * fac + term1 + term2 + term3
Index: /palm/trunk/SOURCE/lpm_pack_arrays.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_pack_arrays.f90	(revision 2416)
+++ /palm/trunk/SOURCE/lpm_pack_arrays.f90	(revision 2417)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! New routine which sorts particles into particles that completed and not 
+! completed the LES timestep. 
+! 
+! 2101 2017-01-05 16:42:31Z suehring
 !
 ! 2000 2016-08-20 18:09:15Z knoop
@@ -61,4 +65,6 @@
 !> deletion and move data with higher index values to these free indices.
 !> Determine the new number of particles.
+!> Moreover, particles are also sorted into groups finished and not finished
+!> its timestep. 
 !------------------------------------------------------------------------------!
  MODULE lpm_pack_arrays_mod
@@ -70,5 +76,5 @@
 
     PRIVATE
-    PUBLIC lpm_pack_all_arrays, lpm_pack_arrays
+    PUBLIC lpm_pack_all_arrays, lpm_pack_arrays, lpm_sort
 
     INTERFACE lpm_pack_all_arrays
@@ -80,5 +86,10 @@
     END INTERFACE lpm_pack_arrays
 
-CONTAINS
+    INTERFACE lpm_sort
+       MODULE PROCEDURE lpm_sort
+    END INTERFACE lpm_sort
+
+
+ CONTAINS
 
 !------------------------------------------------------------------------------!
@@ -236,4 +247,101 @@
     END SUBROUTINE lpm_pack_and_sort
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Sort particles in each sub-grid box into two groups: particles that already
+!> completed the LES timestep, and particles that need further timestepping to
+!> complete the LES timestep. 
+!------------------------------------------------------------------------------!
+    SUBROUTINE lpm_sort
+
+       USE control_parameters,                                                 &
+           ONLY:  dt_3d
+
+       USE indices,                                                            &
+           ONLY: nxl, nxr, nys, nyn, nzb, nzt
+
+       USE kinds
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) :: end_index     !< particle end index for each sub-box
+       INTEGER(iwp) :: i             !< index of particle grid box in x-direction
+       INTEGER(iwp) :: j             !< index of particle grid box in y-direction
+       INTEGER(iwp) :: k             !< index of particle grid box in z-direction
+       INTEGER(iwp) :: n             !< running index for number of particles
+       INTEGER(iwp) :: nb            !< index of subgrid boux
+       INTEGER(iwp) :: nf            !< indices for particles in each sub-box that already finalized their substeps
+       INTEGER(iwp) :: nnf           !< indices for particles in each sub-box that need further treatment
+       INTEGER(iwp) :: num_finalized !< number of particles in each sub-box that already finalized their substeps
+       INTEGER(iwp) :: start_index   !< particle start index for each sub-box
+
+       TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: sort_particles  !< temporary particle array 
+
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb+1, nzt
+
+                number_of_particles = prt_count(k,j,i)
+                IF ( number_of_particles <= 0 )  CYCLE
+
+                particles => grid_particles(k,j,i)%particles(1:number_of_particles)
+
+                DO  nb = 0, 7
+!
+!--                Obtain start and end index for each subgrid box
+                   start_index = grid_particles(k,j,i)%start_index(nb)
+                   end_index   = grid_particles(k,j,i)%end_index(nb)
+!
+!--                Allocate temporary array used for sorting. 
+                   ALLOCATE( sort_particles(start_index:end_index) )
+!
+!--                Determine number of particles already completed the LES 
+!--                timestep, and write them into a temporary array. 
+                   nf = start_index
+                   num_finalized = 0
+                   DO  n = start_index, end_index
+                      IF ( dt_3d - particles(n)%dt_sum < 1E-8_wp )  THEN
+                         sort_particles(nf) = particles(n)
+                         nf                 = nf + 1
+                         num_finalized      = num_finalized + 1
+                      ENDIF
+                   ENDDO
+!
+!--                Determine number of particles that not completed the LES 
+!--                timestep, and write them into a temporary array. 
+                   nnf = nf
+                   DO  n = start_index, end_index
+                      IF ( dt_3d - particles(n)%dt_sum > 1E-8_wp )  THEN
+                         sort_particles(nnf) = particles(n)
+                         nnf                 = nnf + 1
+                      ENDIF
+                   ENDDO
+!
+!--                Write back sorted particles
+                   particles(start_index:end_index) =                          &
+                                           sort_particles(start_index:end_index)
+!
+!--                Determine updated start_index, used to masked already 
+!--                completed particles. 
+                   grid_particles(k,j,i)%start_index(nb) =                     &
+                                      grid_particles(k,j,i)%start_index(nb)    &
+                                    + num_finalized
+!
+!--                Deallocate dummy array
+                   DEALLOCATE ( sort_particles )
+!
+!--                Finally, if number of non-completed particles is non zero 
+!--                in any of the sub-boxes, set control flag appropriately. 
+                   IF ( nnf > nf )                                             &
+                      grid_particles(k,j,i)%time_loop_done = .FALSE.
+
+                ENDDO
+             ENDDO
+          ENDDO
+       ENDDO
+
+    END SUBROUTINE lpm_sort
+
 
  END module lpm_pack_arrays_mod
Index: /palm/trunk/SOURCE/user_lpm_advec.f90
===================================================================
--- /palm/trunk/SOURCE/user_lpm_advec.f90	(revision 2416)
+++ /palm/trunk/SOURCE/user_lpm_advec.f90	(revision 2417)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! Particle loops adapted for sub-box structure, i.e. for each sub-box the 
+! particle loop runs from start_index up to end_index instead from 1 to 
+! number_of_particles.
+! 
+! 2101 2017-01-05 16:42:31Z suehring
 !
 ! 2000 2016-08-20 18:09:15Z knoop
@@ -70,8 +75,12 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  ip   !<
-    INTEGER(iwp) ::  jp   !< 
-    INTEGER(iwp) ::  kp   !< 
-    INTEGER(iwp) ::  n    !< 
+    INTEGER(iwp) ::  ip   !< index of particle grid box, x-direction
+    INTEGER(iwp) ::  jp   !< index of particle grid box, y-direction
+    INTEGER(iwp) ::  kp   !< index of particle grid box, z-direction
+    INTEGER(iwp) ::  n    !< particle index
+    INTEGER(iwp) ::  nb   !< index of sub-box particles are sorted in
+
+    INTEGER(iwp), DIMENSION(0:7)  ::  start_index !< start particle index for current sub-box
+    INTEGER(iwp), DIMENSION(0:7)  ::  end_index   !< start particle index for current sub-box
 
 !
@@ -79,7 +88,14 @@
 !   number_of_particles = prt_count(kp,jp,ip)
 !   particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
+!
+!   start_index = grid_particles(kp,jp,ip)%start_index
+!   end_index   = grid_particles(kp,jp,ip)%end_index
+!
 !   IF ( number_of_particles <= 0 )  CYCLE
 !   DO  n = 1, number_of_particles
-!
+!   DO  nb = 0, 7
+!      DO  n = start_index(nb), end_index(nb)
+!         particles(n)%xxx = 
+!      ENDDO
 !   ENDDO
 
