Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 1928)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 1929)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix in check for use_upstream_for_tke
 !
 ! Former revisions:
@@ -894,5 +894,5 @@
     IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets  .AND.               &
          .NOT. use_upstream_for_tke  .AND.                                       &
-         ( scalar_advec /= 'ws-scheme'  .OR.  scalar_advec /= 'ws-scheme-mono' ) &
+         ( scalar_advec /= 'ws-scheme'  .AND.  scalar_advec /= 'ws-scheme-mono' )&
        )  THEN
        use_upstream_for_tke = .TRUE.
Index: /palm/trunk/SOURCE/lpm.f90
===================================================================
--- /palm/trunk/SOURCE/lpm.f90	(revision 1928)
+++ /palm/trunk/SOURCE/lpm.f90	(revision 1929)
@@ -19,5 +19,6 @@
 ! Current revisions:
 ! ------------------
-!
+! Call wall boundary conditions only if particles are in the vertical range of 
+! topography.
 ! 
 ! Former revisions:
@@ -111,5 +112,5 @@
 
     USE indices,                                                               &
-        ONLY: nxl, nxr, nys, nyn, nzb, nzt
+        ONLY: nxl, nxr, nys, nyn, nzb, nzb_max, nzt, nzb_w_inner
 
     USE kinds
@@ -300,6 +301,8 @@
                 CALL lpm_advec(i,j,k)
 !
-!--             Particle reflection from walls
-                IF ( topography /= 'flat' )  THEN
+!--             Particle reflection from walls. Only applied if the particles 
+!--             are in the vertical range of the topography. (Here, some
+!--             optimization is still possible.)
+                IF ( topography /= 'flat' .AND. k < nzb_max + 2 )  THEN
                    CALL lpm_boundary_conds( 'walls' )
                 ENDIF
@@ -334,5 +337,4 @@
        steps = steps + 1
        dt_3d_reached_l = ALL(grid_particles(:,:,:)%time_loop_done)
-
 !
 !--    Find out, if all particles on every PE have completed the LES timestep
@@ -359,5 +361,4 @@
 !--    Horizontal boundary conditions including exchange between subdmains
        CALL lpm_exchange_horiz
-
 !
 !--    Pack particles (eliminate those marked for deletion),
Index: /palm/trunk/SOURCE/lpm_advec.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_advec.f90	(revision 1928)
+++ /palm/trunk/SOURCE/lpm_advec.f90	(revision 1929)
@@ -19,5 +19,11 @@
 ! Current revisions:
 ! ------------------
-!
+! Put stochastic equation in an extra subroutine.
+! Set flag for stochastic equation to communicate whether a particle is near 
+! topography. This case, memory and drift term are disabled in the Weil equation. 
+!
+! Enable vertical logarithmic interpolation also above topography. This case, 
+! set a lower limit for the friction velocity, as it can become very small 
+! in narrow street canyons, leading to too large particle speeds. 
 ! 
 ! Former revisions:
@@ -110,7 +116,6 @@
         ONLY:  block_offset, c_0, dt_min_part, grid_particles,                 &
                iran_part, log_z_z0, number_of_particles, number_of_sublayers,  &
-               particles, particle_groups, offset_ocean_nzt, sgs_wfu_part,     &
-               sgs_wfv_part, sgs_wfw_part, use_sgs_for_particles,              &
-               vertical_particle_advection, z0_av_global
+               particles, particle_groups, offset_ocean_nzt, sgs_wf_part,      &
+               use_sgs_for_particles, vertical_particle_advection, z0_av_global
         
     USE statistics,                                                            &
@@ -119,48 +124,51 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  agp                         !<
-    INTEGER(iwp) ::  gp_outside_of_building(1:8) !<
-    INTEGER(iwp) ::  i                           !<
-    INTEGER(iwp) ::  ip                          !<
-    INTEGER(iwp) ::  j                           !<
-    INTEGER(iwp) ::  jp                          !<
-    INTEGER(iwp) ::  k                           !<
-    INTEGER(iwp) ::  kp                          !<
-    INTEGER(iwp) ::  kw                          !<
-    INTEGER(iwp) ::  n                           !<
-    INTEGER(iwp) ::  nb                          !<
-    INTEGER(iwp) ::  num_gp                      !<
-
-    INTEGER(iwp), DIMENSION(0:7) ::  start_index !<
-    INTEGER(iwp), DIMENSION(0:7) ::  end_index   !<
-
-    REAL(wp) ::  aa                 !<
-    REAL(wp) ::  bb                 !<
-    REAL(wp) ::  cc                 !<
-    REAL(wp) ::  d_sum              !<
-    REAL(wp) ::  d_z_p_z0           !<
-    REAL(wp) ::  dd                 !<
-    REAL(wp) ::  de_dx_int_l        !<
-    REAL(wp) ::  de_dx_int_u        !<
-    REAL(wp) ::  de_dy_int_l        !<
-    REAL(wp) ::  de_dy_int_u        !<
-    REAL(wp) ::  de_dt              !<
-    REAL(wp) ::  de_dt_min          !<
-    REAL(wp) ::  de_dz_int_l        !<
-    REAL(wp) ::  de_dz_int_u        !<
+    INTEGER(iwp) ::  agp                         !< loop variable
+    INTEGER(iwp) ::  gp_outside_of_building(1:8) !< number of grid points used for particle interpolation in case of topography 
+    INTEGER(iwp) ::  i                           !< index variable along x
+    INTEGER(iwp) ::  ip                          !< index variable along x 
+    INTEGER(iwp) ::  ilog                        !< index variable along x 
+    INTEGER(iwp) ::  j                           !< index variable along y 
+    INTEGER(iwp) ::  jp                          !< index variable along y 
+    INTEGER(iwp) ::  jlog                        !< index variable along y 
+    INTEGER(iwp) ::  k                           !< index variable along z 
+    INTEGER(iwp) ::  kp                          !< index variable along z 
+    INTEGER(iwp) ::  kw                          !< index variable along z
+    INTEGER(iwp) ::  n                           !< loop variable over all particles in a grid box
+    INTEGER(iwp) ::  nb                          !< block number particles are sorted in
+    INTEGER(iwp) ::  num_gp                      !< number of adjacent grid points inside topography 
+
+    INTEGER(iwp), DIMENSION(0:7) ::  start_index !< start particle index for current block
+    INTEGER(iwp), DIMENSION(0:7) ::  end_index   !< start particle index for current block
+
+    REAL(wp) ::  aa                 !< dummy argument for horizontal particle interpolation 
+    REAL(wp) ::  bb                 !< dummy argument for horizontal particle interpolation 
+    REAL(wp) ::  cc                 !< dummy argument for horizontal particle interpolation 
+    REAL(wp) ::  d_sum              !< dummy argument for horizontal particle interpolation in case of topography 
+    REAL(wp) ::  d_z_p_z0           !< inverse of interpolation length for logarithmic interpolation 
+    REAL(wp) ::  dd                 !< dummy argument for horizontal particle interpolation 
+    REAL(wp) ::  de_dx_int_l        !< x/y-interpolated TKE gradient (x) at particle position at lower vertical level
+    REAL(wp) ::  de_dx_int_u        !< x/y-interpolated TKE gradient (x) at particle position at upper vertical level
+    REAL(wp) ::  de_dy_int_l        !< x/y-interpolated TKE gradient (y) at particle position at lower vertical level
+    REAL(wp) ::  de_dy_int_u        !< x/y-interpolated TKE gradient (y) at particle position at upper vertical level
+    REAL(wp) ::  de_dt              !< temporal derivative of TKE experienced by the particle
+    REAL(wp) ::  de_dt_min          !< lower level for temporal TKE derivative
+    REAL(wp) ::  de_dz_int_l        !< x/y-interpolated TKE gradient (z) at particle position at lower vertical level
+    REAL(wp) ::  de_dz_int_u        !< x/y-interpolated TKE gradient (z) at particle position at upper vertical level
     REAL(wp) ::  diameter           !< diamter of droplet
-    REAL(wp) ::  diss_int_l         !<
-    REAL(wp) ::  diss_int_u         !<
-    REAL(wp) ::  dt_gap             !<
-    REAL(wp) ::  dt_particle_m      !<
-    REAL(wp) ::  e_int_l            !<
-    REAL(wp) ::  e_int_u            !<
-    REAL(wp) ::  e_mean_int         !<
+    REAL(wp) ::  diss_int_l         !< x/y-interpolated dissipation at particle position at lower vertical level
+    REAL(wp) ::  diss_int_u         !< x/y-interpolated dissipation at particle position at upper vertical level
+    REAL(wp) ::  dt_gap             !< remaining time until particle time integration reaches LES time
+    REAL(wp) ::  dt_particle_m      !< previous particle time step
+    REAL(wp) ::  e_int_l            !< x/y-interpolated TKE at particle position at lower vertical level
+    REAL(wp) ::  e_int_u            !< x/y-interpolated TKE at particle position at upper vertical level
+    REAL(wp) ::  e_mean_int         !< horizontal mean TKE at particle height
     REAL(wp) ::  exp_arg            !<
     REAL(wp) ::  exp_term           !<
-    REAL(wp) ::  gg                 !<
-    REAL(wp) ::  height_p           !<
+    REAL(wp) ::  gg                 !< dummy argument for horizontal particle interpolation 
+    REAL(wp) ::  height_p           !< dummy argument for logarithmic interpolation
     REAL(wp) ::  lagr_timescale     !< Lagrangian timescale
-    REAL(wp) ::  location(1:30,1:3) !<
+    REAL(wp) ::  location(1:30,1:3) !< wall locations
+    REAL(wp) ::  log_z_z0_int       !< logarithmus used for surface_layer interpolation
     REAL(wp) ::  random_gauss       !<
     REAL(wp) ::  RL                 !< Lagrangian autocorrelation coefficient
@@ -169,16 +177,16 @@
     REAL(wp) ::  rg3                !< Gaussian distributed random number
     REAL(wp) ::  sigma              !< velocity standard deviation
-    REAL(wp) ::  u_int_l            !<
-    REAL(wp) ::  u_int_u            !<
-    REAL(wp) ::  us_int             !<
-    REAL(wp) ::  v_int_l            !<
-    REAL(wp) ::  v_int_u            !<
+    REAL(wp) ::  u_int_l            !< x/y-interpolated u-component at particle position at lower vertical level
+    REAL(wp) ::  u_int_u            !< x/y-interpolated u-component at particle position at upper vertical level
+    REAL(wp) ::  us_int             !< friction velocity at particle grid box
+    REAL(wp) ::  v_int_l            !< x/y-interpolated v-component at particle position at lower vertical level
+    REAL(wp) ::  v_int_u            !< x/y-interpolated v-component at particle position at upper vertical level
     REAL(wp) ::  vv_int             !<
-    REAL(wp) ::  w_int_l            !<
-    REAL(wp) ::  w_int_u            !<
+    REAL(wp) ::  w_int_l            !< x/y-interpolated w-component at particle position at lower vertical level
+    REAL(wp) ::  w_int_u            !< x/y-interpolated w-component at particle position at upper vertical level
     REAL(wp) ::  w_s                !< terminal velocity of droplets
-    REAL(wp) ::  x                  !<
-    REAL(wp) ::  y                  !<
-    REAL(wp) ::  z_p                !<
+    REAL(wp) ::  x                  !< dummy argument for horizontal particle interpolation 
+    REAL(wp) ::  y                  !< dummy argument for horizontal particle interpolation 
+    REAL(wp) ::  z_p                !< surface layer height (0.5 dz)
 
     REAL(wp), PARAMETER ::  a_rog = 9.65_wp      !< parameter for fall velocity
@@ -189,28 +197,28 @@
     REAL(wp), PARAMETER ::  d0_rog = 0.745_wp    !< separation diameter
 
-    REAL(wp), DIMENSION(1:30) ::  d_gp_pl !<
-    REAL(wp), DIMENSION(1:30) ::  de_dxi  !<
-    REAL(wp), DIMENSION(1:30) ::  de_dyi  !<
-    REAL(wp), DIMENSION(1:30) ::  de_dzi  !<
-    REAL(wp), DIMENSION(1:30) ::  dissi   !<
-    REAL(wp), DIMENSION(1:30) ::  ei      !<
-
+    REAL(wp), DIMENSION(1:30) ::  d_gp_pl !< dummy argument for particle interpolation scheme in case of topography
+    REAL(wp), DIMENSION(1:30) ::  de_dxi  !< horizontal TKE gradient along x at adjacent wall 
+    REAL(wp), DIMENSION(1:30) ::  de_dyi  !< horizontal TKE gradient along y at adjacent wall 
+    REAL(wp), DIMENSION(1:30) ::  de_dzi  !< horizontal TKE gradient along z at adjacent wall 
+    REAL(wp), DIMENSION(1:30) ::  dissi   !< dissipation at adjacent wall
+    REAL(wp), DIMENSION(1:30) ::  ei      !< TKE at adjacent wall
+
+    REAL(wp), DIMENSION(number_of_particles) ::  term_1_2     !< flag to communicate whether a particle is near topography or not
     REAL(wp), DIMENSION(number_of_particles) ::  dens_ratio   !<
-    REAL(wp), DIMENSION(number_of_particles) ::  de_dx_int    !<
-    REAL(wp), DIMENSION(number_of_particles) ::  de_dy_int    !<
-    REAL(wp), DIMENSION(number_of_particles) ::  de_dz_int    !<
-    REAL(wp), DIMENSION(number_of_particles) ::  diss_int     !<
-    REAL(wp), DIMENSION(number_of_particles) ::  dt_particle  !<
-    REAL(wp), DIMENSION(number_of_particles) ::  e_int        !<
-    REAL(wp), DIMENSION(number_of_particles) ::  fs_int       !<
-    REAL(wp), DIMENSION(number_of_particles) ::  log_z_z0_int !<
-    REAL(wp), DIMENSION(number_of_particles) ::  u_int        !<
-    REAL(wp), DIMENSION(number_of_particles) ::  v_int        !<
-    REAL(wp), DIMENSION(number_of_particles) ::  w_int        !<
-    REAL(wp), DIMENSION(number_of_particles) ::  xv           !<
-    REAL(wp), DIMENSION(number_of_particles) ::  yv           !<
-    REAL(wp), DIMENSION(number_of_particles) ::  zv           !<
-
-    REAL(wp), DIMENSION(number_of_particles, 3) ::  rg !<
+    REAL(wp), DIMENSION(number_of_particles) ::  de_dx_int    !< horizontal TKE gradient along x at particle position
+    REAL(wp), DIMENSION(number_of_particles) ::  de_dy_int    !< horizontal TKE gradient along y at particle position
+    REAL(wp), DIMENSION(number_of_particles) ::  de_dz_int    !< horizontal TKE gradient along z at particle position
+    REAL(wp), DIMENSION(number_of_particles) ::  diss_int     !< dissipation at particle position
+    REAL(wp), DIMENSION(number_of_particles) ::  dt_particle  !< particle time step
+    REAL(wp), DIMENSION(number_of_particles) ::  e_int        !< TKE at particle position
+    REAL(wp), DIMENSION(number_of_particles) ::  fs_int       !< weighting factor for subgrid-scale particle speed
+    REAL(wp), DIMENSION(number_of_particles) ::  u_int        !< u-component of particle speed
+    REAL(wp), DIMENSION(number_of_particles) ::  v_int        !< v-component of particle speed 
+    REAL(wp), DIMENSION(number_of_particles) ::  w_int        !< w-component of particle speed
+    REAL(wp), DIMENSION(number_of_particles) ::  xv           !< x-position
+    REAL(wp), DIMENSION(number_of_particles) ::  yv           !< y-position
+    REAL(wp), DIMENSION(number_of_particles) ::  zv           !< z-position
+
+    REAL(wp), DIMENSION(number_of_particles, 3) ::  rg !< vector of Gaussian distributed random numbers 
 
     CALL cpu_log( log_point_s(44), 'lpm_advec', 'continue' )
@@ -236,4 +244,5 @@
        j = jp + block_offset(nb)%j_off
        k = kp + block_offset(nb)%k_off
+
 
 !
@@ -249,13 +258,16 @@
 !--       First, check if particle is located below first vertical grid level 
 !--       (Prandtl-layer height)
-          IF ( constant_flux_layer  .AND.  particles(n)%z < z_p )  THEN
+          ilog = ( particles(n)%x + 0.5_wp * dx ) * ddx
+          jlog = ( particles(n)%y + 0.5_wp * dy ) * ddy
+
+          IF ( constant_flux_layer .AND. zv(n) - zw(nzb_s_inner(jlog,ilog)) < z_p )  THEN
 !
 !--          Resolved-scale horizontal particle velocity is zero below z0. 
-             IF ( particles(n)%z < z0_av_global )  THEN
+             IF ( zv(n) - zw(nzb_s_inner(jlog,ilog)) < z0_av_global )  THEN
                 u_int(n) = 0.0_wp
              ELSE
 !
 !--             Determine the sublayer. Further used as index.
-                height_p     = ( particles(n)%z - z0_av_global )            &
+                height_p     = ( zv(n) - zw(nzb_s_inner(jlog,ilog)) - z0_av_global )            &
                                      * REAL( number_of_sublayers, KIND=wp ) &
                                      * d_z_p_z0 
@@ -263,9 +275,15 @@
 !--             Calculate LOG(z/z0) for exact particle height. Therefore,   
 !--             interpolate linearly between precalculated logarithm. 
-                log_z_z0_int(n) = log_z_z0(INT(height_p))                      &
+                log_z_z0_int = log_z_z0(INT(height_p))                         &
                                  + ( height_p - INT(height_p) )                &
                                  * ( log_z_z0(INT(height_p)+1)                 &
                                       - log_z_z0(INT(height_p))                &
                                    ) 
+!
+!--             Limit friction velocity. In narrow canyons or holes the 
+!--             friction velocity can become very small, resulting in a too
+!--             large particle speed.
+                us_int   = MAX( 0.5_wp * ( us(jlog,ilog) + us(jlog,ilog-1) ),  &
+                                0.01_wp )  
 !
 !--             Neutral solution is applied for all situations, e.g. also for 
@@ -275,10 +293,7 @@
 !--             as sensitivity studies revealed no significant effect of 
 !--             using the neutral solution also for un/stable situations.
-!--             Calculated left and bottom index on u grid.
-                us_int   = 0.5_wp * ( us(j,i) + us(j,i-1) )  
-
-                u_int(n) = -usws(j,i) / ( us_int * kappa + 1E-10_wp )          & 
-                         * log_z_z0_int(n) - u_gtrans
-
+                u_int(n) = -usws(jlog,ilog) / ( us_int * kappa + 1E-10_wp )          & 
+                            * log_z_z0_int - u_gtrans
+                
              ENDIF
 !
@@ -308,4 +323,5 @@
                            ( u_int_u - u_int_l )
              ENDIF
+
           ENDIF
 
@@ -319,11 +335,35 @@
        DO  n = start_index(nb), end_index(nb)
 
-          IF ( constant_flux_layer  .AND.  particles(n)%z < z_p )  THEN
-
-             IF ( particles(n)%z < z0_av_global )  THEN
+          ilog = ( particles(n)%x + 0.5_wp * dx ) * ddx
+          jlog = ( particles(n)%y + 0.5_wp * dy ) * ddy
+          IF ( constant_flux_layer .AND. zv(n) - zw(nzb_s_inner(jlog,ilog)) < z_p )  THEN
+
+             IF ( zv(n) - zw(nzb_s_inner(jlog,ilog)) < z0_av_global )  THEN
 !
 !--             Resolved-scale horizontal particle velocity is zero below z0. 
                 v_int(n) = 0.0_wp
              ELSE       
+!
+!--             Determine the sublayer. Further used as index. Please note, 
+!--             logarithmus can not be reused from above, as in in case of
+!--             topography particle on u-grid can be above surface-layer height,
+!--             whereas it can be below on v-grid. 
+                height_p     = ( zv(n) - zw(nzb_s_inner(jlog,ilog)) - z0_av_global ) &
+                                  * REAL( number_of_sublayers, KIND=wp )    &
+                                  * d_z_p_z0 
+!
+!--             Calculate LOG(z/z0) for exact particle height. Therefore,   
+!--             interpolate linearly between precalculated logarithm. 
+                log_z_z0_int = log_z_z0(INT(height_p))                         &
+                                 + ( height_p - INT(height_p) )                &
+                                 * ( log_z_z0(INT(height_p)+1)                 &
+                                      - log_z_z0(INT(height_p))                &
+                                   ) 
+!
+!--             Limit friction velocity. In narrow canyons or holes the 
+!--             friction velocity can become very small, resulting in a too
+!--             large particle speed.
+                us_int   = MAX( 0.5_wp * ( us(jlog,ilog) + us(jlog-1,ilog) ), &
+                                0.01_wp )    
 !
 !--             Neutral solution is applied for all situations, e.g. also for 
@@ -333,10 +373,9 @@
 !--             as sensitivity studies revealed no significant effect of 
 !--             using the neutral solution also for un/stable situations.
-!--             Calculated left and bottom index on v grid.
-                us_int   = 0.5_wp * ( us(j,i) + us(j-1,i) )    
-
-                v_int(n) = -vsws(j,i) / ( us_int * kappa + 1E-10_wp )          &
-                         * log_z_z0_int(n) - v_gtrans
-             ENDIF
+                v_int(n) = -vsws(jlog,ilog) / ( us_int * kappa + 1E-10_wp )          &
+                         * log_z_z0_int - v_gtrans
+
+             ENDIF
+
           ELSE
              x  = xv(n) - i * dx
@@ -361,4 +400,5 @@
                                   ( v_int_u - v_int_l )
              ENDIF
+
           ENDIF
 
@@ -367,5 +407,5 @@
        i = ip + block_offset(nb)%i_off
        j = jp + block_offset(nb)%j_off
-       k = kp-1
+       k = kp - 1
 !
 !--    Same procedure for interpolation of the w velocity-component
@@ -535,4 +575,8 @@
                 ENDIF
 
+!
+!--             Set flag for stochastic equation.
+                term_1_2(n) = 1.0_wp
+
              ENDDO
           ENDDO
@@ -541,5 +585,4 @@
 
           DO  n = 1, number_of_particles
-
              i = particles(n)%x * ddx
              j = particles(n)%y * ddy
@@ -575,7 +618,5 @@
                 de_dzi(num_gp) = de_dz(k,j,i)
              ENDIF
-
-             IF ( k > nzb_s_inner(j+1,i)  .OR.  nzb_s_inner(j+1,i) == 0 ) &
-             THEN
+             IF ( k > nzb_s_inner(j+1,i)  .OR.  nzb_s_inner(j+1,i) == 0 )  THEN
                 num_gp = num_gp + 1
                 gp_outside_of_building(2) = 1
@@ -603,6 +644,5 @@
              ENDIF
 
-             IF ( k+1 > nzb_s_inner(j+1,i)  .OR.  nzb_s_inner(j+1,i) == 0 ) &
-             THEN
+             IF ( k+1 > nzb_s_inner(j+1,i)  .OR.  nzb_s_inner(j+1,i) == 0 )  THEN
                 num_gp = num_gp + 1
                 gp_outside_of_building(4) = 1
@@ -617,6 +657,5 @@
              ENDIF
 
-             IF ( k > nzb_s_inner(j,i+1)  .OR.  nzb_s_inner(j,i+1) == 0 ) &
-             THEN
+             IF ( k > nzb_s_inner(j,i+1)  .OR.  nzb_s_inner(j,i+1) == 0 )  THEN
                 num_gp = num_gp + 1
                 gp_outside_of_building(5) = 1
@@ -631,6 +670,5 @@
              ENDIF
 
-             IF ( k > nzb_s_inner(j+1,i+1)  .OR.  nzb_s_inner(j+1,i+1) == 0 ) &
-             THEN
+             IF ( k > nzb_s_inner(j+1,i+1)  .OR.  nzb_s_inner(j+1,i+1) == 0 )  THEN
                 num_gp = num_gp + 1
                 gp_outside_of_building(6) = 1
@@ -645,6 +683,5 @@
              ENDIF
 
-             IF ( k+1 > nzb_s_inner(j,i+1)  .OR.  nzb_s_inner(j,i+1) == 0 ) &
-             THEN
+             IF ( k+1 > nzb_s_inner(j,i+1)  .OR.  nzb_s_inner(j,i+1) == 0 )  THEN
                 num_gp = num_gp + 1
                 gp_outside_of_building(7) = 1
@@ -659,6 +696,5 @@
              ENDIF
 
-             IF ( k+1 > nzb_s_inner(j+1,i+1)  .OR.  nzb_s_inner(j+1,i+1) == 0)&
-             THEN
+             IF ( k+1 > nzb_s_inner(j+1,i+1)  .OR.  nzb_s_inner(j+1,i+1) == 0)  THEN
                 num_gp = num_gp + 1
                 gp_outside_of_building(8) = 1
@@ -672,5 +708,4 @@
                 de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
              ENDIF
-
 !
 !--          If all gridpoints are situated outside of a building, then the
@@ -685,9 +720,9 @@
                 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
@@ -699,7 +734,7 @@
                              ) / ( 3.0_wp * gg )
                    e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
-                                     ( e_int_u - e_int_l )
-                ENDIF
-!
+                                       ( e_int_u - e_int_l )
+                ENDIF
+! 
 !--             Needed to avoid NaN particle velocities (this might not be
 !--             required any more)
@@ -757,5 +792,5 @@
                                    ( 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
@@ -778,5 +813,5 @@
                                ( gg - dd ) * diss(k,j+1,i+1) &
                              ) / ( 3.0_wp * gg )
-
+ 
                 IF ( k == nzt )  THEN
                    diss_int(n) = diss_int_l
@@ -790,7 +825,10 @@
                                            ( 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
@@ -810,5 +848,5 @@
 
                 IF ( gp_outside_of_building(5) == 1  .AND. &
-                   gp_outside_of_building(1) == 0 )  THEN
+                     gp_outside_of_building(1) == 0 )  THEN
                    num_gp = num_gp + 1
                    location(num_gp,1) = i * dx + 0.5_wp * dx
@@ -893,5 +931,5 @@
                    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)
+                   de_dzi(num_gp) = de_dz(k,j,i)  
                 ENDIF
 
@@ -1076,5 +1114,5 @@
                 ENDIF
 
-!
+! 
 !--             If wall between gridpoint 6 and gridpoint 8, then
 !--             Neumann boundary condition has to be applied
@@ -1092,21 +1130,25 @@
                    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)
+                   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
@@ -1118,12 +1160,12 @@
                       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 
+                   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) ) * &
@@ -1138,7 +1180,15 @@
                                          de_dzi(agp) / ( (num_gp-1) * d_sum )
                    ENDDO
-
-                ENDIF
-
+ 
+                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
           ENDDO
@@ -1209,6 +1259,6 @@
 !
 !--       Calculate the Lagrangian timescale according to Weil et al. (2004).
-          lagr_timescale = ( 4.0_wp * e_int(n) ) / &
-                           ( 3.0_wp * fs_int(n) * c_0 * diss_int(n) )
+          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 )
 
 !
@@ -1233,9 +1283,9 @@
 !--          [-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_wfu_part * e_int(n) ) *   &
+             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_wfv_part * e_int(n) ) *   &
+             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_wfw_part * e_int(n) ) *   &
+             particles(n)%rvar3 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) *   &
                                   ( rg(n,3) - 1.0_wp )
 
@@ -1267,45 +1317,21 @@
              ENDIF
 
-             particles(n)%rvar1 = particles(n)%rvar1 - fs_int(n) * c_0 *       &
-                           diss_int(n) * particles(n)%rvar1 * dt_particle(n) / &
-                           ( 4.0_wp * sgs_wfu_part * e_int(n) ) +              &
-                           ( 2.0_wp * sgs_wfu_part * de_dt *                   &
-                             particles(n)%rvar1 /                              &
-                             ( 2.0_wp * sgs_wfu_part * e_int(n) ) +            &
-                             de_dx_int(n)                                      &
-                           ) * dt_particle(n) / 2.0_wp          +              &
-                           SQRT( fs_int(n) * c_0 * diss_int(n) ) *             &
-                           ( rg(n,1) - 1.0_wp ) *                              &
-                           SQRT( dt_particle(n) )
-
-             particles(n)%rvar2 = particles(n)%rvar2 - fs_int(n) * c_0 *       &
-                           diss_int(n) * particles(n)%rvar2 * dt_particle(n) / &
-                           ( 4.0_wp * sgs_wfv_part * e_int(n) ) +              &
-                           ( 2.0_wp * sgs_wfv_part * de_dt *                   &
-                             particles(n)%rvar2 /                              &
-                             ( 2.0_wp * sgs_wfv_part * e_int(n) ) +            &
-                             de_dy_int(n)                                      &
-                           ) * dt_particle(n) / 2.0_wp          +              &
-                           SQRT( fs_int(n) * c_0 * diss_int(n) ) *             &
-                           ( rg(n,2) - 1.0_wp ) *                              &
-                           SQRT( dt_particle(n) )
-
-             particles(n)%rvar3 = particles(n)%rvar3 - fs_int(n) * c_0 *       &
-                           diss_int(n) * particles(n)%rvar3 * dt_particle(n) / &
-                           ( 4.0_wp * sgs_wfw_part * e_int(n) ) +              &
-                           ( 2.0_wp * sgs_wfw_part * de_dt *                   &
-                             particles(n)%rvar3 /                              &
-                             ( 2.0_wp * sgs_wfw_part * e_int(n) ) +            &
-                             de_dz_int(n)                                      &
-                           ) * dt_particle(n) / 2.0_wp          +              &
-                           SQRT( fs_int(n) * c_0 * diss_int(n) ) *             &
-                           ( rg(n,3) - 1.0_wp ) *                              &
-                           SQRT( dt_particle(n) )
+             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
@@ -1507,3 +1533,66 @@
     CALL cpu_log( log_point_s(44), 'lpm_advec', 'pause' )
 
+
  END SUBROUTINE lpm_advec
+
+! Description:
+! ------------
+!> Calculation of subgrid-scale particle speed using the stochastic model 
+!> of Weil et al. (2004, JAS, 61, 2877-2887).
+!------------------------------------------------------------------------------! 
+ SUBROUTINE weil_stochastic_eq( v_sgs, fs_n, e_n, dedxi_n, dedt_n, diss_n,     &
+                                dt_n, rg_n, fac )
+
+    USE kinds
+
+    USE particle_attributes,                                                   &
+        ONLY:  c_0, sgs_wf_part
+
+    IMPLICIT NONE
+
+    REAL(wp) ::  a1      !< dummy argument
+    REAL(wp) ::  dedt_n  !< time derivative of TKE at particle position 
+    REAL(wp) ::  dedxi_n !< horizontal derivative of TKE at particle position
+    REAL(wp) ::  diss_n  !< dissipation at particle position 
+    REAL(wp) ::  dt_n    !< particle timestep
+    REAL(wp) ::  e_n     !< TKE at particle position
+    REAL(wp) ::  fac     !< flag to identify adjacent topography
+    REAL(wp) ::  fs_n    !< weighting factor to prevent that subgrid-scale particle speed becomes too large
+    REAL(wp) ::  sgs_w   !< constant (1/3)
+    REAL(wp) ::  rg_n    !< random number
+    REAL(wp) ::  term1   !< memory term
+    REAL(wp) ::  term2   !< drift correction term
+    REAL(wp) ::  term3   !< random term
+    REAL(wp) ::  v_sgs   !< subgrid-scale velocity component 
+
+!
+!-- Please note, terms 1 and 2 (drift and memory term, respectively) are 
+!-- multiplied by a flag to switch of both terms near topography. 
+!-- This is necessary, as both terms may cause a subgrid-scale velocity build up 
+!-- if particles are trapped in regions with very small TKE, e.g. in narrow street 
+!-- canyons resolved by only a few grid points. Hence, term 1 and term 2 are 
+!-- disabled if one of the adjacent grid points belongs to topography. 
+!-- Moreover, in this case, the  previous subgrid-scale component is also set
+!-- to zero.
+
+    a1 = fs_n * c_0 * diss_n
+!
+!-- Memory term
+    term1 = - a1 * v_sgs * dt_n / ( 4.0_wp * sgs_wf_part * e_n + 1E-20_wp )    &
+                 * fac
+!
+!-- Drift correction term
+    term2 = ( ( dedt_n * v_sgs / e_n ) + dedxi_n ) * 0.5_wp * dt_n              &
+                 * fac
+!
+!-- Random term
+    term3 = SQRT( MAX( a1, 1E-20 ) ) * ( rg_n - 1.0_wp ) * SQRT( dt_n )
+!
+!-- In cese one of the adjacent grid-boxes belongs to topograhy, the previous
+!-- 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
+
+ END SUBROUTINE weil_stochastic_eq
Index: /palm/trunk/SOURCE/lpm_boundary_conds.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_boundary_conds.f90	(revision 1928)
+++ /palm/trunk/SOURCE/lpm_boundary_conds.f90	(revision 1929)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-!
+! Rewritten wall reflection 
 ! 
 ! Former revisions:
@@ -102,45 +102,61 @@
     CHARACTER (LEN=*) ::  range     !<
     
-    INTEGER(iwp) ::  i              !<
-    INTEGER(iwp) ::  inc            !<
-    INTEGER(iwp) ::  ir             !<
-    INTEGER(iwp) ::  i1             !<
-    INTEGER(iwp) ::  i2             !<
-    INTEGER(iwp) ::  i3             !<
-    INTEGER(iwp) ::  i5             !<
-    INTEGER(iwp) ::  j              !<
-    INTEGER(iwp) ::  jr             !<
-    INTEGER(iwp) ::  j1             !<
-    INTEGER(iwp) ::  j2             !<
-    INTEGER(iwp) ::  j3             !<
-    INTEGER(iwp) ::  j5             !<
-    INTEGER(iwp) ::  k              !<
-    INTEGER(iwp) ::  k1             !<
-    INTEGER(iwp) ::  k2             !<
-    INTEGER(iwp) ::  k3             !<
-    INTEGER(iwp) ::  k5             !<
-    INTEGER(iwp) ::  n              !<
-    INTEGER(iwp) ::  t_index        !<
-    INTEGER(iwp) ::  t_index_number !<
+    INTEGER(iwp) ::  inc            !< dummy for sorting algorithmus
+    INTEGER(iwp) ::  ir             !< dummy for sorting algorithmus
+    INTEGER(iwp) ::  i1             !< grid index (x) of old particle position
+    INTEGER(iwp) ::  i2             !< grid index (x) of current particle position
+    INTEGER(iwp) ::  i3             !< grid index (x) of intermediate particle position
+    INTEGER(iwp) ::  jr             !< dummy for sorting algorithmus
+    INTEGER(iwp) ::  j1             !< grid index (y) of old particle position
+    INTEGER(iwp) ::  j2             !< grid index (x) of current particle position
+    INTEGER(iwp) ::  j3             !< grid index (x) of intermediate particle position
+    INTEGER(iwp) ::  n              !< particle number
+    INTEGER(iwp) ::  t_index        !< running index for intermediate particle timesteps in reflection algorithmus
+    INTEGER(iwp) ::  t_index_number !< number of intermediate particle timesteps in reflection algorithmus
+    INTEGER(iwp) ::  tmp_x          !< dummy for sorting algorithmus
+    INTEGER(iwp) ::  tmp_y          !< dummy for sorting algorithmus
+
+    INTEGER(iwp), DIMENSION(0:10) :: x_ind(0:10) = 0 !< index array (x) of intermediate particle positions 
+    INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (x) of intermediate particle positions 
     
-    LOGICAL  ::  reflect_x   !<
-    LOGICAL  ::  reflect_y   !<
-    LOGICAL  ::  reflect_z   !<
-
-    REAL(wp) ::  dt_particle !<
-    REAL(wp) ::  pos_x       !<
-    REAL(wp) ::  pos_x_old   !<
-    REAL(wp) ::  pos_y       !<
-    REAL(wp) ::  pos_y_old   !<
-    REAL(wp) ::  pos_z       !<
-    REAL(wp) ::  pos_z_old   !<
-    REAL(wp) ::  prt_x       !<
-    REAL(wp) ::  prt_y       !<
-    REAL(wp) ::  prt_z       !<
-    REAL(wp) ::  t(1:200)    !<
-    REAL(wp) ::  tmp_t       !<
-    REAL(wp) ::  xline       !<
-    REAL(wp) ::  yline       !<
-    REAL(wp) ::  zline       !<
+    LOGICAL  ::  cross_wall_x    !< flag to check if particle reflection along x is necessary
+    LOGICAL  ::  cross_wall_y    !< flag to check if particle reflection along y is necessary
+    LOGICAL  ::  downwards       !< flag to check if particle reflection along z is necessary (only if particle move downwards)
+    LOGICAL  ::  reflect_x       !< flag to check if particle is already reflected along x
+    LOGICAL  ::  reflect_y       !< flag to check if particle is already reflected along y
+    LOGICAL  ::  reflect_z       !< flag to check if particle is already reflected along z
+    LOGICAL  ::  tmp_reach_x     !< dummy for sorting algorithmus
+    LOGICAL  ::  tmp_reach_y     !< dummy for sorting algorithmus
+    LOGICAL  ::  tmp_reach_z     !< dummy for sorting algorithmus
+    LOGICAL  ::  x_wall_reached  !< flag to check if particle has already reached wall
+    LOGICAL  ::  y_wall_reached  !< flag to check if particle has already reached wall
+
+    LOGICAL, DIMENSION(0:10) ::  reach_x  !< flag to check if particle is at a yz-wall
+    LOGICAL, DIMENSION(0:10) ::  reach_y  !< flag to check if particle is at a xz-wall
+    LOGICAL, DIMENSION(0:10) ::  reach_z  !< flag to check if particle is at a xy-wall
+
+    REAL(wp) ::  dt_particle    !< particle timestep
+    REAL(wp) ::  dum            !< dummy argument
+    REAL(wp) ::  eps = 1E-10_wp !< security number to check if particle has reached a wall
+    REAL(wp) ::  pos_x          !< intermediate particle position (x)
+    REAL(wp) ::  pos_x_old      !< particle position (x) at previous particle timestep
+    REAL(wp) ::  pos_y          !< intermediate particle position (y)
+    REAL(wp) ::  pos_y_old      !< particle position (y) at previous particle timestep
+    REAL(wp) ::  pos_z          !< intermediate particle position (z)
+    REAL(wp) ::  pos_z_old      !< particle position (z) at previous particle timestep
+    REAL(wp) ::  prt_x          !< current particle position (x)
+    REAL(wp) ::  prt_y          !< current particle position (y)
+    REAL(wp) ::  prt_z          !< current particle position (z)
+    REAL(wp) ::  t_old          !< previous reflection time
+    REAL(wp) ::  tmp_t          !< dummy for sorting algorithmus
+    REAL(wp) ::  xwall          !< location of wall in x
+    REAL(wp) ::  ywall          !< location of wall in y
+    REAL(wp) ::  zwall1         !< location of wall in z (old grid box)
+    REAL(wp) ::  zwall2         !< location of wall in z (current grid box)
+    REAL(wp) ::  zwall3         !< location of wall in z (old y, current x)
+    REAL(wp) ::  zwall4         !< location of wall in z (current y, old x)
+
+    REAL(wp), DIMENSION(0:10) ::  t  !< reflection time
+
 
     IF ( range == 'bottom/top' )  THEN
@@ -211,532 +227,332 @@
     ELSEIF ( range == 'walls' )  THEN
 
+
        CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'start' )
 
-       reflect_x = .FALSE.
-       reflect_y = .FALSE.
-       reflect_z = .FALSE.
-
        DO  n = 1, number_of_particles
-
+!
+!--       Recalculate particle timestep
           dt_particle = particles(n)%age - particles(n)%age_m
-
+!
+!--       Obtain x/y indices for current particle position
           i2 = ( particles(n)%x + 0.5_wp * dx ) * ddx
           j2 = ( particles(n)%y + 0.5_wp * dy ) * ddy
-          k2 = particles(n)%z / dz + 1 + offset_ocean_nzt_m1
-
+!
+!--       Save current particle positions
           prt_x = particles(n)%x
           prt_y = particles(n)%y
           prt_z = particles(n)%z
-
-!
-!--       If the particle position is below the surface, it has to be reflected
-          IF ( k2 <= nzb_s_inner(j2,i2)  .AND.  nzb_s_inner(j2,i2) /=0 )  THEN
-
-             pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle
-             pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle
-             pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle
-             i1 = ( pos_x_old + 0.5_wp * dx ) * ddx
-             j1 = ( pos_y_old + 0.5_wp * dy ) * ddy
-             k1 = pos_z_old / dz + offset_ocean_nzt_m1
-
-!
-!--          Case 1
-             IF ( particles(n)%x > pos_x_old  .AND.  particles(n)%y > pos_y_old )&
-             THEN
-                t_index = 1
-
-                DO  i = i1, i2
-                   xline      = i * dx + 0.5_wp * dx
-                   t(t_index) = ( xline - pos_x_old ) / &
-                                ( particles(n)%x - pos_x_old )
-                   t_index    = t_index + 1
+!
+!--       Recalculate old particle positions
+          pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle
+          pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle
+          pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle
+!
+!--       Obtain x/y indices for old particle positions
+          i1 = ( pos_x_old + 0.5_wp * dx ) * ddx
+          j1 = ( pos_y_old + 0.5_wp * dy ) * ddy
+!
+!--       Determine horizontal as well as vertical walls at which particle can 
+!--       be potentially reflected. 
+!--       Start with walls aligned in yz layer.
+!--       Wall to the right 
+          IF ( prt_x > pos_x_old )  THEN
+             xwall = ( i1 + 0.5_wp ) * dx
+!
+!--       Wall to the left
+          ELSE
+             xwall = ( i1 - 0.5_wp ) * dx
+          ENDIF
+!
+!--       Walls aligned in xz layer
+!--       Wall to the north
+          IF ( prt_y > pos_y_old )  THEN
+             ywall = ( j1 + 0.5_wp ) * dy
+!--       Wall to the south
+          ELSE
+             ywall = ( j1 - 0.5_wp ) * dy
+          ENDIF
+!
+!--       Walls aligned in xy layer at which particle can be possiblly reflected
+          zwall1 = zw(nzb_s_inner(j2,i2))
+          zwall2 = zw(nzb_s_inner(j1,i1))
+          zwall3 = zw(nzb_s_inner(j1,i2))
+          zwall4 = zw(nzb_s_inner(j2,i1))
+!
+!--       Initialize flags to check if particle reflection is necessary
+          downwards    = .FALSE.
+          cross_wall_x = .FALSE.
+          cross_wall_y = .FALSE.
+!
+!--       Initialize flags to check if a wall is reached
+          reach_x      = .FALSE.
+          reach_y      = .FALSE.
+          reach_z      = .FALSE.
+!
+!--       Initialize flags to check if a particle was already reflected
+          reflect_x = .FALSE.
+          reflect_y = .FALSE.
+          reflect_z = .FALSE.
+!
+!--       Initialize flags to check if a vertical wall is already crossed.
+!--       ( Required to obtain correct indices. )
+          x_wall_reached = .FALSE.
+          y_wall_reached = .FALSE.
+!
+!--       Initialize time array 
+          t     = 0.0_wp
+!
+!--       Check if particle can reach any wall. This case, calculate the 
+!--       fractional time needed to reach this wall. Store this fractional
+!--       timestep in array t. Moreover, store indices for these grid
+!--       boxes where the respective wall belongs to.  
+!--       Start with x-direction.
+          t_index    = 1
+          t(t_index) = ( xwall - pos_x_old )                                   &
+                     / MERGE( MAX( prt_x - pos_x_old,  1E-30_wp ),             &
+                              MIN( prt_x - pos_x_old, -1E-30_wp ),             &
+                              prt_x > pos_x_old )
+          x_ind(t_index)   = i2
+          y_ind(t_index)   = j1
+          reach_x(t_index) = .TRUE.
+          reach_y(t_index) = .FALSE.
+          reach_z(t_index) = .FALSE.
+!
+!--       Store these values only if particle really reaches any wall. t must 
+!--       be in a interval between [0:1].
+          IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )  THEN
+             t_index      = t_index + 1
+             cross_wall_x = .TRUE.
+          ENDIF
+!
+!--       y-direction
+          t(t_index) = ( ywall - pos_y_old )                                   &
+                     / MERGE( MAX( prt_y - pos_y_old,  1E-30_wp ),             &
+                              MIN( prt_y - pos_y_old, -1E-30_wp ),             &
+                              prt_y > pos_y_old )
+          x_ind(t_index)   = i1
+          y_ind(t_index)   = j2
+          reach_x(t_index) = .FALSE.
+          reach_y(t_index) = .TRUE.
+          reach_z(t_index) = .FALSE.
+          IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )  THEN
+             t_index      = t_index + 1
+             cross_wall_y = .TRUE.
+          ENDIF
+!
+!--       z-direction
+!--       At first, check if particle moves downwards. Only in this case a 
+!--       particle can be reflected vertically.
+          IF ( prt_z < pos_z_old )  THEN
+
+             downwards = .TRUE.
+             dum       =  1.0_wp / MERGE( MAX( prt_z - pos_z_old,  1E-30_wp ), &
+                                          MIN( prt_z - pos_z_old, -1E-30_wp ), &
+                                          prt_z > pos_z_old )
+
+             t(t_index)       = ( zwall1 - pos_z_old ) * dum
+             x_ind(t_index)   = i2
+             y_ind(t_index)   = j2
+             reach_x(t_index) = .FALSE.
+             reach_y(t_index) = .FALSE.
+             reach_z(t_index) = .TRUE.
+             IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )            &
+                t_index = t_index + 1
+
+             reach_x(t_index) = .FALSE.
+             reach_y(t_index) = .FALSE.
+             reach_z(t_index) = .TRUE.
+             t(t_index)       = ( zwall2 - pos_z_old ) * dum
+             x_ind(t_index)   = i1
+             y_ind(t_index)   = j1
+             IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )            &
+                t_index = t_index + 1
+
+             reach_x(t_index) = .FALSE.
+             reach_y(t_index) = .FALSE.
+             reach_z(t_index) = .TRUE.
+             t(t_index)       = ( zwall3 - pos_z_old ) * dum
+             x_ind(t_index)   = i2
+             y_ind(t_index)   = j1
+             IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )            &
+                t_index = t_index + 1
+
+             reach_x(t_index) = .FALSE.
+             reach_y(t_index) = .FALSE.
+             reach_z(t_index) = .TRUE.
+             t(t_index)       = ( zwall4 - pos_z_old ) * dum
+             x_ind(t_index)   = i1
+             y_ind(t_index)   = j2
+             IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )            &
+                t_index = t_index + 1
+
+          END IF
+          t_index_number = t_index - 1
+!
+!--       Carry out reflection only if particle reaches any wall
+          IF ( cross_wall_x .OR. cross_wall_y .OR. downwards )  THEN
+!
+!--          Sort fractional timesteps in ascending order. Also sort the 
+!--          corresponding indices and flag according to the time interval a  
+!--          particle reaches the respective wall.
+             inc = 1
+             jr  = 1
+             DO WHILE ( inc <= t_index_number )
+                inc = 3 * inc + 1
+             ENDDO
+
+             DO WHILE ( inc > 1 )
+                inc = inc / 3
+                DO  ir = inc+1, t_index_number
+                   tmp_t       = t(ir)
+                   tmp_x       = x_ind(ir)
+                   tmp_y       = y_ind(ir)
+                   tmp_reach_x = reach_x(ir)
+                   tmp_reach_y = reach_y(ir)
+                   tmp_reach_z = reach_z(ir)
+                   jr    = ir
+                   DO WHILE ( t(jr-inc) > tmp_t )
+                      t(jr)       = t(jr-inc)
+                      x_ind(jr)   = x_ind(jr-inc)
+                      y_ind(jr)   = y_ind(jr-inc)
+                      reach_x(jr) = reach_x(jr-inc)
+                      reach_y(jr) = reach_y(jr-inc)
+                      reach_z(jr) = reach_z(jr-inc)
+                      jr    = jr - inc
+                      IF ( jr <= inc )  EXIT
+                   ENDDO
+                   t(jr)       = tmp_t
+                   x_ind(jr)   = tmp_x
+                   y_ind(jr)   = tmp_y
+                   reach_x(jr) = tmp_reach_x
+                   reach_y(jr) = tmp_reach_y
+                   reach_z(jr) = tmp_reach_z
                 ENDDO
-
-                DO  j = j1, j2
-                   yline      = j * dy + 0.5_wp * dy
-                   t(t_index) = ( yline - pos_y_old ) / &
-                                ( particles(n)%y - pos_y_old )
-                   t_index    = t_index + 1
-                ENDDO
-
-                IF ( particles(n)%z < pos_z_old )  THEN
-                   DO  k = k1, k2 , -1
-                      zline      = k * dz
-                      t(t_index) = ( pos_z_old - zline ) / &
-                                   ( pos_z_old - particles(n)%z )
-                      t_index    = t_index + 1
-                   ENDDO
+             ENDDO
+!
+!--          Initialize temporary particle positions
+             pos_x = pos_x_old
+             pos_y = pos_y_old
+             pos_z = pos_z_old
+!
+!--          Loop over all times a particle possibly moves into a new grid box
+             t_old = 0.0_wp
+             DO t_index = 1, t_index_number 
+!           
+!--             Calculate intermediate particle position according to the 
+!--             timesteps a particle reaches any wall.
+                pos_x = pos_x + ( t(t_index) - t_old ) * dt_particle           &
+                                                       * particles(n)%speed_x
+                pos_y = pos_y + ( t(t_index) - t_old ) * dt_particle           &
+                                                       * particles(n)%speed_y
+                pos_z = pos_z + ( t(t_index) - t_old ) * dt_particle           &
+                                                       * particles(n)%speed_z
+!
+!--             Obtain x/y grid indices for intermediate particle position from
+!--             sorted index array
+                i3 = x_ind(t_index)
+                j3 = y_ind(t_index)
+!
+!--             Check which wall is already reached
+                IF ( .NOT. x_wall_reached )  x_wall_reached = reach_x(t_index) 
+                IF ( .NOT. y_wall_reached )  y_wall_reached = reach_y(t_index) 
+!
+!--             Check if a particle needs to be reflected at any yz-wall. If 
+!--             necessary, carry out reflection. Please note, a security 
+!--             constant is required, as the particle position do not 
+!--             necessarily exactly match the wall location due to rounding 
+!--             errors.   
+                IF ( ABS( pos_x - xwall ) < eps      .AND.                     &
+                     pos_z <= zw(nzb_s_inner(j3,i3)) .AND.                     &
+                     reach_x(t_index)                .AND.                     &
+                     .NOT. reflect_x )  THEN
+!
+!--                Reflection in x-direction. 
+!--                Ensure correct reflection by MIN/MAX functions, depending on
+!--                direction of particle transport. 
+!--                Due to rounding errors pos_x do not exactly matches the wall
+!--                location, leading to erroneous reflection.             
+                   pos_x = MERGE( MIN( 2.0_wp * xwall - pos_x, xwall ),        &
+                                  MAX( 2.0_wp * xwall - pos_x, xwall ),        &
+                                  particles(n)%x > xwall )
+!
+!--                Change sign of particle speed                     
+                   particles(n)%speed_x = - particles(n)%speed_x
+!
+!--                Change also sign of subgrid-scale particle speed
+                   particles(n)%rvar1 = - particles(n)%rvar1
+!
+!--                Set flag that reflection along x is already done
+                   reflect_x          = .TRUE.
+!
+!--                As particle do not crosses any further yz-wall during 
+!--                this timestep, set further x-indices to the current one.
+                   x_ind(t_index:t_index_number) = i1
+!
+!--             If particle already reached the wall but was not reflected, 
+!--             set further x-indices to the new one.
+                ELSEIF ( x_wall_reached .AND. .NOT. reflect_x )  THEN
+                    x_ind(t_index:t_index_number) = i2
                 ENDIF
-
-                t_index_number = t_index - 1
-
-!
-!--             Sorting t
-                inc = 1
-                jr  = 1
-                DO WHILE ( inc <= t_index_number )
-                   inc = 3 * inc + 1
-                ENDDO
-
-                DO WHILE ( inc > 1 )
-                   inc = inc / 3
-                   DO  ir = inc+1, t_index_number
-                      tmp_t = t(ir)
-                      jr    = ir
-                      DO WHILE ( t(jr-inc) > tmp_t )
-                         t(jr) = t(jr-inc)
-                         jr    = jr - inc
-                         IF ( jr <= inc )  EXIT
-                      ENDDO
-                      t(jr) = tmp_t
-                   ENDDO
-                ENDDO
-
-         case1: DO  t_index = 1, t_index_number
-
-                   pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
-                   pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
-                   pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
-
-                   i3 = ( pos_x + 0.5_wp * dx ) * ddx   
-                   j3 = ( pos_y + 0.5_wp * dy ) * ddy
-                   k3 = pos_z / dz + offset_ocean_nzt_m1
-
-                   i5 = pos_x * ddx
-                   j5 = pos_y * ddy
-                   k5 = pos_z / dz + offset_ocean_nzt_m1
-
-                   IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
-                        nzb_s_inner(j5,i3) /= 0 )  THEN
-
-                      IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
-                           k3 == nzb_s_inner(j5,i3) )  THEN
-                         reflect_z = .TRUE.
-                         EXIT case1
-                      ENDIF
+!
+!--             Check if a particle needs to be reflected at any xz-wall. If 
+!--             necessary, carry out reflection.
+                IF ( ABS( pos_y - ywall ) < eps      .AND.                     &
+                     pos_z <= zw(nzb_s_inner(j3,i3)) .AND.                     &
+                     reach_y(t_index)                .AND.                     &
+                     .NOT. reflect_y ) THEN
+
+                   pos_y = MERGE( MIN( 2.0_wp * ywall - pos_y, ywall ),        &
+                                  MAX( 2.0_wp * ywall - pos_y, ywall ),        &
+                                  particles(n)%y > ywall ) 
+
+                   particles(n)%speed_y = - particles(n)%speed_y      
+                   particles(n)%rvar2   = - particles(n)%rvar2       
+     
+                   reflect_y            = .TRUE.
+                   y_ind(t_index:t_index_number) = j1
+
+                ELSEIF ( y_wall_reached .AND. .NOT. reflect_y )  THEN
+                   y_ind(t_index:t_index_number) = j2
+                ENDIF
+!
+!--             Check if a particle needs to be reflected at any xy-wall. If 
+!--             necessary, carry out reflection.
+                IF ( downwards .AND. reach_z(t_index) .AND.                    &
+                     .NOT. reflect_z )  THEN
+                   IF ( pos_z - zw(nzb_s_inner(j3,i3)) < eps ) THEN
+  
+                      pos_z = MAX( 2.0_wp * zw(nzb_s_inner(j3,i3)) - pos_z,    &
+                                   zw(nzb_s_inner(j3,i3)) )
+
+                      particles(n)%speed_z = - particles(n)%speed_z
+                      particles(n)%rvar3   = - particles(n)%rvar3
+
+                      reflect_z            = .TRUE.
 
                    ENDIF
 
-                   IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
-                        nzb_s_inner(j3,i5) /= 0 )  THEN
-
-                      IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
-                           k3 == nzb_s_inner(j3,i5) )  THEN
-                         reflect_z = .TRUE.
-                         EXIT case1
-                      ENDIF
-
-                   ENDIF
-
-                   IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
-                        nzb_s_inner(j3,i3) /= 0 )  THEN
-
-                      IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
-                           k3 == nzb_s_inner(j3,i3) )  THEN
-                         reflect_z = .TRUE.
-                         EXIT case1
-                      ENDIF
-
-                      IF ( pos_y == ( j3 * dy - 0.5_wp * dy )  .AND. &
-                           pos_z < nzb_s_inner(j3,i3) * dz )  THEN
-                         reflect_y = .TRUE.
-                         EXIT case1
-                      ENDIF
-
-                      IF ( pos_x == ( i3 * dx - 0.5_wp * dx )  .AND. &
-                           pos_z < nzb_s_inner(j3,i3) * dz )  THEN
-                         reflect_x = .TRUE.
-                         EXIT case1
-                      ENDIF
-
-                   ENDIF
-
-                ENDDO case1
-
-!
-!--          Case 2
-             ELSEIF ( particles(n)%x > pos_x_old  .AND. &
-                      particles(n)%y < pos_y_old )  THEN
-
-                t_index = 1
-
-                DO  i = i1, i2
-                   xline      = i * dx + 0.5_wp * dx
-                   t(t_index) = ( xline - pos_x_old ) / &
-                                ( particles(n)%x - pos_x_old )
-                   t_index    = t_index + 1
-                ENDDO
-
-                DO  j = j1, j2, -1
-                   yline      = j * dy - 0.5_wp * dy
-                   t(t_index) = ( pos_y_old - yline ) / &
-                                ( pos_y_old - particles(n)%y )
-                   t_index    = t_index + 1
-                ENDDO
-
-                IF ( particles(n)%z < pos_z_old )  THEN
-                   DO  k = k1, k2 , -1
-                      zline      = k * dz
-                      t(t_index) = ( pos_z_old - zline ) / &
-                                   ( pos_z_old - particles(n)%z )
-                      t_index    = t_index + 1
-                   ENDDO
                 ENDIF
-                t_index_number = t_index-1
-
-!
-!--             Sorting t
-                inc = 1
-                jr  = 1
-                DO WHILE ( inc <= t_index_number )
-                   inc = 3 * inc + 1
-                ENDDO
-
-                DO WHILE ( inc > 1 )
-                   inc = inc / 3
-                   DO  ir = inc+1, t_index_number
-                      tmp_t = t(ir)
-                      jr    = ir
-                      DO WHILE ( t(jr-inc) > tmp_t )
-                         t(jr) = t(jr-inc)
-                         jr    = jr - inc
-                         IF ( jr <= inc )  EXIT
-                      ENDDO
-                      t(jr) = tmp_t
-                   ENDDO
-                ENDDO
-
-         case2: DO  t_index = 1, t_index_number
-
-                   pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
-                   pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
-                   pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
-
-                   i3 = ( pos_x + 0.5_wp * dx ) * ddx
-                   j3 = ( pos_y + 0.5_wp * dy ) * ddy
-                   k3 = pos_z / dz + offset_ocean_nzt_m1
-
-                   i5 = pos_x * ddx
-                   j5 = pos_y * ddy
-                   k5 = pos_z / dz + offset_ocean_nzt_m1
-
-                   IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
-                        nzb_s_inner(j3,i5) /= 0 )  THEN
-
-                      IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
-                           k3 == nzb_s_inner(j3,i5) )  THEN
-                         reflect_z = .TRUE.
-                         EXIT case2
-                      ENDIF
-
-                   ENDIF
-
-                   IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
-                        nzb_s_inner(j3,i3) /= 0 )  THEN
-
-                      IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
-                           k3 == nzb_s_inner(j3,i3) )  THEN
-                         reflect_z = .TRUE.
-                         EXIT case2
-                      ENDIF
-
-                      IF ( pos_x == ( i3 * dx - 0.5_wp * dx )  .AND. &
-                           pos_z < nzb_s_inner(j3,i3) * dz )  THEN
-                         reflect_x = .TRUE.
-                         EXIT case2
-                      ENDIF
-
-                   ENDIF
-
-                   IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
-                        nzb_s_inner(j5,i3) /= 0 )  THEN
-
-                      IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
-                           k3 == nzb_s_inner(j5,i3) )  THEN
-                         reflect_z = .TRUE.
-                         EXIT case2
-                      ENDIF
-
-                      IF ( pos_y == ( j5 * dy + 0.5_wp * dy )  .AND. &
-                           pos_z < nzb_s_inner(j5,i3) * dz )  THEN
-                         reflect_y = .TRUE.
-                         EXIT case2
-                      ENDIF
-
-                   ENDIF
-
-                ENDDO case2
-
-!
-!--          Case 3
-             ELSEIF ( particles(n)%x < pos_x_old  .AND. &
-                      particles(n)%y > pos_y_old )  THEN
-
-                t_index = 1
-
-                DO  i = i1, i2, -1
-                   xline      = i * dx - 0.5_wp * dx
-                   t(t_index) = ( pos_x_old - xline ) / &
-                                ( pos_x_old - particles(n)%x )
-                   t_index    = t_index + 1
-                ENDDO
-
-                DO  j = j1, j2
-                   yline      = j * dy + 0.5_wp * dy
-                   t(t_index) = ( yline - pos_y_old ) / &
-                                ( particles(n)%y - pos_y_old )
-                   t_index    = t_index + 1
-                ENDDO
-
-                IF ( particles(n)%z < pos_z_old )  THEN
-                   DO  k = k1, k2 , -1
-                      zline      = k * dz
-                      t(t_index) = ( pos_z_old - zline ) / &
-                                   ( pos_z_old - particles(n)%z )
-                      t_index    = t_index + 1
-                   ENDDO
-                ENDIF
-                t_index_number = t_index - 1
-
-!
-!--             Sorting t
-                inc = 1
-                jr  = 1
-
-                DO WHILE ( inc <= t_index_number )
-                   inc = 3 * inc + 1
-                ENDDO
-
-                DO WHILE ( inc > 1 )
-                   inc = inc / 3
-                   DO  ir = inc+1, t_index_number
-                      tmp_t = t(ir)
-                      jr    = ir
-                      DO WHILE ( t(jr-inc) > tmp_t )
-                         t(jr) = t(jr-inc)
-                         jr    = jr - inc
-                         IF ( jr <= inc )  EXIT
-                      ENDDO
-                      t(jr) = tmp_t
-                   ENDDO
-                ENDDO
-
-         case3: DO  t_index = 1, t_index_number
-
-                   pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
-                   pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
-                   pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
-
-                   i3 = ( pos_x + 0.5_wp * dx ) * ddx
-                   j3 = ( pos_y + 0.5_wp * dy ) * ddy
-                   k3 = pos_z / dz + offset_ocean_nzt_m1
-
-                   i5 = pos_x * ddx
-                   j5 = pos_y * ddy
-                   k5 = pos_z / dz + offset_ocean_nzt_m1
-
-                   IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
-                        nzb_s_inner(j5,i3) /= 0 )  THEN
-
-                      IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
-                           k3 == nzb_s_inner(j5,i3) )  THEN
-                         reflect_z = .TRUE.
-                         EXIT case3
-                      ENDIF
-
-                   ENDIF
-
-                   IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
-                        nzb_s_inner(j3,i3) /= 0 )  THEN
-
-                      IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
-                           k3 == nzb_s_inner(j3,i3) )  THEN
-                         reflect_z = .TRUE.
-                         EXIT case3
-                      ENDIF
-
-                      IF ( pos_y == ( j3 * dy - 0.5_wp * dy )  .AND. &
-                           pos_z < nzb_s_inner(j3,i3) * dz )  THEN
-                         reflect_y = .TRUE.
-                         EXIT case3
-                      ENDIF
-
-                   ENDIF
-
-                   IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
-                        nzb_s_inner(j3,i5) /= 0 )  THEN
-
-                      IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
-                           k3 == nzb_s_inner(j3,i5) )  THEN
-                         reflect_z = .TRUE.
-                         EXIT case3
-                      ENDIF
-
-                      IF ( pos_x == ( i5 * dx + 0.5_wp * dx )  .AND. &
-                           pos_z < nzb_s_inner(j3,i5) * dz )  THEN
-                         reflect_x = .TRUE.
-                         EXIT case3
-                      ENDIF
-
-                   ENDIF
-
-                ENDDO case3
-
-!
-!--          Case 4
-             ELSEIF ( particles(n)%x < pos_x_old  .AND. &
-                      particles(n)%y < pos_y_old )  THEN
-
-                t_index = 1
-
-                DO  i = i1, i2, -1
-                   xline      = i * dx - 0.5_wp * dx
-                   t(t_index) = ( pos_x_old - xline ) / &
-                                ( pos_x_old - particles(n)%x )
-                   t_index    = t_index + 1
-                ENDDO
-
-                DO  j = j1, j2, -1
-                   yline      = j * dy - 0.5_wp * dy
-                   t(t_index) = ( pos_y_old - yline ) / &
-                                ( pos_y_old - particles(n)%y )
-                   t_index    = t_index + 1
-                ENDDO
-
-                IF ( particles(n)%z < pos_z_old )  THEN
-                   DO  k = k1, k2 , -1
-                      zline      = k * dz
-                      t(t_index) = ( pos_z_old - zline ) / &
-                                   ( pos_z_old-particles(n)%z )
-                      t_index    = t_index + 1
-                   ENDDO
-                ENDIF
-                t_index_number = t_index-1
-
-!
-!--             Sorting t
-                inc = 1
-                jr  = 1
-
-                DO WHILE ( inc <= t_index_number )
-                   inc = 3 * inc + 1
-                ENDDO
-
-                DO WHILE ( inc > 1 )
-                   inc = inc / 3
-                   DO  ir = inc+1, t_index_number
-                      tmp_t = t(ir)
-                      jr = ir
-                      DO WHILE ( t(jr-inc) > tmp_t )
-                         t(jr) = t(jr-inc)
-                         jr    = jr - inc
-                         IF ( jr <= inc )  EXIT
-                      ENDDO
-                      t(jr) = tmp_t
-                   ENDDO
-                ENDDO
-
-         case4: DO  t_index = 1, t_index_number
-
-                   pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
-                   pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
-                   pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
-
-                   i3 = ( pos_x + 0.5_wp * dx ) * ddx   
-                   j3 = ( pos_y + 0.5_wp * dy ) * ddy
-                   k3 = pos_z / dz + offset_ocean_nzt_m1
-
-                   i5 = pos_x * ddx
-                   j5 = pos_y * ddy
-                   k5 = pos_z / dz + offset_ocean_nzt_m1
-
-                   IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
-                        nzb_s_inner(j3,i3) /= 0 )  THEN
-
-                      IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
-                           k3 == nzb_s_inner(j3,i3) )  THEN
-                         reflect_z = .TRUE.
-                         EXIT case4
-                      ENDIF
-
-                   ENDIF
-
-                   IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
-                        nzb_s_inner(j3,i5) /= 0 )  THEN
-
-                      IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
-                           k3 == nzb_s_inner(j3,i5) )  THEN
-                         reflect_z = .TRUE.
-                         EXIT case4
-                      ENDIF
-
-                      IF ( pos_x == ( i5 * dx + 0.5_wp * dx )  .AND. &
-                           nzb_s_inner(j3,i5) /=0  .AND.          &
-                           pos_z < nzb_s_inner(j3,i5) * dz )  THEN
-                         reflect_x = .TRUE.
-                         EXIT case4
-                      ENDIF
-
-                   ENDIF
-
-                   IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
-                        nzb_s_inner(j5,i3) /= 0 )  THEN
-
-                      IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
-                           k5 == nzb_s_inner(j5,i3) )  THEN
-                         reflect_z = .TRUE.
-                         EXIT case4
-                      ENDIF
-
-                      IF ( pos_y == ( j5 * dy + 0.5_wp * dy )  .AND. &
-                           nzb_s_inner(j5,i3) /= 0  .AND.         &
-                           pos_z < nzb_s_inner(j5,i3) * dz )  THEN
-                         reflect_y = .TRUE.
-                         EXIT case4
-                      ENDIF
-
-                   ENDIF
- 
-                ENDDO case4
+!
+!--             Swap time
+                t_old = t(t_index)
+
+             ENDDO
+!
+!--          If a particle was reflected, calculate final position from last
+!--          intermediate position.
+             IF ( reflect_x .OR. reflect_y .OR. reflect_z )  THEN
+
+                particles(n)%x = pos_x + ( 1.0_wp - t_old ) * dt_particle      &
+                                                         * particles(n)%speed_x
+                particles(n)%y = pos_y + ( 1.0_wp - t_old ) * dt_particle      &
+                                                         * particles(n)%speed_y
+                particles(n)%z = pos_z + ( 1.0_wp - t_old ) * dt_particle      &
+                                                         * particles(n)%speed_z
 
              ENDIF
 
-          ENDIF      ! Check, if particle position is below the surface
-
-!
-!--       Do the respective reflection, in case that one of the above conditions
-!--       is found to be true
-          IF ( reflect_z )  THEN
-
-             particles(n)%z       = 2.0_wp * pos_z - prt_z
-             particles(n)%speed_z = - particles(n)%speed_z
-
-             IF ( use_sgs_for_particles )  THEN
-                particles(n)%rvar3 = - particles(n)%rvar3
-             ENDIF
-             reflect_z = .FALSE.
-
-          ELSEIF ( reflect_y )  THEN
-
-             particles(n)%y       = 2.0_wp * pos_y - prt_y
-             particles(n)%speed_y = - particles(n)%speed_y
-
-             IF ( use_sgs_for_particles )  THEN
-                particles(n)%rvar2 = - particles(n)%rvar2
-             ENDIF
-             reflect_y = .FALSE.
-
-          ELSEIF ( reflect_x )  THEN
-
-             particles(n)%x       = 2.0_wp * pos_x - prt_x
-             particles(n)%speed_x = - particles(n)%speed_x
-
-             IF ( use_sgs_for_particles )  THEN 
-                particles(n)%rvar1 = - particles(n)%rvar1
-             ENDIF
-             reflect_x = .FALSE.
-
-          ENDIF        
-    
+          ENDIF
+
        ENDDO
 
Index: /palm/trunk/SOURCE/lpm_exchange_horiz.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_exchange_horiz.f90	(revision 1928)
+++ /palm/trunk/SOURCE/lpm_exchange_horiz.f90	(revision 1929)
@@ -19,5 +19,13 @@
 ! Current revisions:
 ! ------------------
-! 
+! Bugfixes: 
+! - reallocation of new particles 
+!   ( did not work for small number of min_nr_particle )
+! - dynamical reallocation of north-south exchange arrays ( particles got lost )
+! - north-south exchange ( nr_move_north, nr_move_south were overwritten by zero )
+! - horizontal particle boundary conditions in serial mode
+!
+! Remove unused variables
+! Descriptions in variable declaration blocks added
 ! 
 ! Former revisions:
@@ -121,6 +129,6 @@
     INTEGER(iwp)            ::  nr_move_south               !<
 
-    TYPE(particle_type), DIMENSION(NR_2_direction_move) ::  move_also_north
-    TYPE(particle_type), DIMENSION(NR_2_direction_move) ::  move_also_south
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  move_also_north
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  move_also_south
 
     SAVE
@@ -161,35 +169,27 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  i                                       !<
-    INTEGER(iwp) ::  ip                                      !<
-    INTEGER(iwp) ::  j                                       !<
-    INTEGER(iwp) ::  jp                                      !<
-    INTEGER(iwp) ::  kp                                      !<
-    INTEGER(iwp) ::  n                                       !<
-    INTEGER(iwp) ::  trlp_count                              !<
-    INTEGER(iwp) ::  trlp_count_recv                         !<
-    INTEGER(iwp) ::  trlpt_count                             !<
-    INTEGER(iwp) ::  trlpt_count_recv                        !<
-    INTEGER(iwp) ::  trnp_count                              !<
-    INTEGER(iwp) ::  trnp_count_recv                         !<
-    INTEGER(iwp) ::  trnpt_count                             !<
-    INTEGER(iwp) ::  trnpt_count_recv                        !<
-    INTEGER(iwp) ::  trrp_count                              !<
-    INTEGER(iwp) ::  trrp_count_recv                         !<
-    INTEGER(iwp) ::  trrpt_count                             !<
-    INTEGER(iwp) ::  trrpt_count_recv                        !<
-    INTEGER(iwp) ::  trsp_count                              !<
-    INTEGER(iwp) ::  trsp_count_recv                         !<
-    INTEGER(iwp) ::  trspt_count                             !<
-    INTEGER(iwp) ::  trspt_count_recv                        !<
-
-    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvlp  !<
-    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvnp  !<
-    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvrp  !<
-    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvsp  !<
-    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trlp  !<
-    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trnp  !<
-    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trrp  !<
-    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trsp  !<
+    INTEGER(iwp) ::  i                 !< grid index (x) of particle positition
+    INTEGER(iwp) ::  ip                !< index variable along x
+    INTEGER(iwp) ::  j                 !< grid index (y) of particle positition
+    INTEGER(iwp) ::  jp                !< index variable along y
+    INTEGER(iwp) ::  kp                !< index variable along z
+    INTEGER(iwp) ::  n                 !< particle index variable 
+    INTEGER(iwp) ::  trlp_count        !< number of particles send to left PE
+    INTEGER(iwp) ::  trlp_count_recv   !< number of particles receive from right PE
+    INTEGER(iwp) ::  trnp_count        !< number of particles send to north PE
+    INTEGER(iwp) ::  trnp_count_recv   !< number of particles receive from south PE
+    INTEGER(iwp) ::  trrp_count        !< number of particles send to right PE
+    INTEGER(iwp) ::  trrp_count_recv   !< number of particles receive from left PE
+    INTEGER(iwp) ::  trsp_count        !< number of particles send to south PE
+    INTEGER(iwp) ::  trsp_count_recv   !< number of particles receive from north PE
+
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvlp  !< particles received from right PE
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvnp  !< particles received from south PE
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvrp  !< particles received from left PE
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvsp  !< particles received from north PE
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trlp  !< particles send to left PE
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trnp  !< particles send to north PE
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trrp  !< particles send to right PE
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trsp  !< particles send to south PE
 
     CALL cpu_log( log_point_s(23), 'lpm_exchange_horiz', 'start' )
@@ -209,12 +209,8 @@
 !-- be adjusted.
     trlp_count  = 0
-    trlpt_count = 0
     trrp_count  = 0
-    trrpt_count = 0
 
     trlp_count_recv   = 0
-    trlpt_count_recv  = 0
     trrp_count_recv   = 0
-    trrpt_count_recv  = 0
 
     IF ( pdims(1) /= 1 )  THEN
@@ -232,6 +228,6 @@
                    IF ( particles(n)%particle_mask )  THEN
                       i = ( particles(n)%x + 0.5_wp * dx ) * ddx
-   !
-   !--                Above calculation does not work for indices less than zero
+!
+!--                   Above calculation does not work for indices less than zero
                       IF ( particles(n)%x < -0.5_wp * dx )  i = -1
 
@@ -249,7 +245,5 @@
 
        IF ( trlp_count  == 0 )  trlp_count  = 1
-       IF ( trlpt_count == 0 )  trlpt_count = 1
        IF ( trrp_count  == 0 )  trrp_count  = 1
-       IF ( trrpt_count == 0 )  trrpt_count = 1
 
        ALLOCATE( trlp(trlp_count), trrp(trrp_count) )
@@ -259,12 +253,10 @@
 
        trlp_count  = 0
-       trlpt_count = 0
        trrp_count  = 0
-       trrpt_count = 0
 
     ENDIF
 !
 !-- Compute only first (nxl) and last (nxr) loop iterration
-    DO  ip = nxl, nxr,nxr-nxl
+    DO  ip = nxl, nxr, nxr-nxl
        DO  jp = nys, nyn
           DO  kp = nzb+1, nzt
@@ -279,15 +271,15 @@
 
                    i = ( particles(n)%x + 0.5_wp * dx ) * ddx
-   !
-   !--             Above calculation does not work for indices less than zero
+!
+!--                Above calculation does not work for indices less than zero
                    IF ( particles(n)%x < - 0.5_wp * dx )  i = -1
 
                    IF ( i <  nxl )  THEN
                       IF ( i < 0 )  THEN
-   !
-   !--                   Apply boundary condition along x
+!
+!--                   Apply boundary condition along x
                          IF ( ibc_par_lr == 0 )  THEN
-   !
-   !--                      Cyclic condition
+!
+!--                         Cyclic condition
                             IF ( pdims(1) == 1 )  THEN
                                particles(n)%x        = ( nx + 1 ) * dx + particles(n)%x
@@ -312,12 +304,12 @@
 
                          ELSEIF ( ibc_par_lr == 1 )  THEN
-   !
-   !--                      Particle absorption
+!
+!--                         Particle absorption
                             particles(n)%particle_mask = .FALSE.
                             deleted_particles = deleted_particles + 1
 
                          ELSEIF ( ibc_par_lr == 2 )  THEN
-   !
-   !--                      Particle reflection
+!
+!--                         Particle reflection
                             particles(n)%x       = -particles(n)%x
                             particles(n)%speed_x = -particles(n)%speed_x
@@ -325,7 +317,7 @@
                          ENDIF
                       ELSE
-   !
-   !--                   Store particle data in the transfer array, which will be 
-   !--                   send to the neighbouring PE
+!
+!--                      Store particle data in the transfer array, which will be 
+!--                      send to the neighbouring PE
                          trlp_count = trlp_count + 1
                          trlp(trlp_count) = particles(n)
@@ -337,9 +329,9 @@
                    ELSEIF ( i > nxr )  THEN
                       IF ( i > nx )  THEN
-   !
-   !--                   Apply boundary condition along x
+!
+!--                      Apply boundary condition along x
                          IF ( ibc_par_lr == 0 )  THEN
-   !
-   !--                      Cyclic condition
+!
+!--                         Cyclic condition
                             IF ( pdims(1) == 1 )  THEN
                                particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
@@ -358,12 +350,12 @@
 
                          ELSEIF ( ibc_par_lr == 1 )  THEN
-   !
-   !--                      Particle absorption
+!
+!--                         Particle absorption
                             particles(n)%particle_mask = .FALSE.
                             deleted_particles = deleted_particles + 1
 
                          ELSEIF ( ibc_par_lr == 2 )  THEN
-   !
-   !--                      Particle reflection
+!
+!--                         Particle reflection
                             particles(n)%x       = 2 * ( nx * dx ) - particles(n)%x
                             particles(n)%speed_x = -particles(n)%speed_x
@@ -371,7 +363,7 @@
                          ENDIF
                       ELSE
-   !
-   !--                   Store particle data in the transfer array, which will be send
-   !--                   to the neighbouring PE
+!
+!--                      Store particle data in the transfer array, which will be send
+!--                      to the neighbouring PE
                          trrp_count = trrp_count + 1
                          trrp(trrp_count) = particles(n)
@@ -383,4 +375,5 @@
                    ENDIF
                 ENDIF
+
              ENDDO
           ENDDO
@@ -389,4 +382,12 @@
 
 !
+!-- Allocate arrays required for north-south exchange, as these
+!-- are used directly after particles are exchange along x-direction.
+    ALLOCATE( move_also_north(1:NR_2_direction_move) )
+    ALLOCATE( move_also_south(1:NR_2_direction_move) )
+
+    nr_move_north = 0
+    nr_move_south = 0
+!
 !-- Send left boundary, receive right boundary (but first exchange how many
 !-- and check, if particle storage must be extended)
@@ -427,5 +428,4 @@
 
     ENDIF
-
 
 !
@@ -435,15 +435,11 @@
 !-- Find out first the number of particles to be transferred and allocate
 !-- temporary arrays needed to store them.
-!-- For a one-dimensional decomposition along x, no transfer is necessary,
+!-- For a one-dimensional decomposition along y, no transfer is necessary,
 !-- because the particle remains on the PE.
     trsp_count  = nr_move_south
-    trspt_count = 0
     trnp_count  = nr_move_north
-    trnpt_count = 0
 
     trsp_count_recv   = 0
-    trspt_count_recv  = 0
     trnp_count_recv   = 0
-    trnpt_count_recv  = 0
 
     IF ( pdims(2) /= 1 )  THEN
@@ -452,5 +448,5 @@
 !--    data
        DO  ip = nxl, nxr
-          DO  jp = nys, nyn, nyn-nys                                 !compute only first (nys) and last (nyn) loop iterration
+          DO  jp = nys, nyn, nyn-nys    !compute only first (nys) and last (nyn) loop iterration
              DO  kp = nzb+1, nzt
                 number_of_particles = prt_count(kp,jp,ip)
@@ -476,7 +472,5 @@
 
        IF ( trsp_count  == 0 )  trsp_count  = 1
-       IF ( trspt_count == 0 )  trspt_count = 1
        IF ( trnp_count  == 0 )  trnp_count  = 1
-       IF ( trnpt_count == 0 )  trnpt_count = 1
 
        ALLOCATE( trsp(trsp_count), trnp(trnp_count) )
@@ -486,8 +480,6 @@
 
        trsp_count  = nr_move_south
-       trspt_count = 0
        trnp_count  = nr_move_north
-       trnpt_count = 0
-
+       
        trsp(1:nr_move_south) = move_also_south(1:nr_move_south)
        trnp(1:nr_move_north) = move_also_north(1:nr_move_north)
@@ -506,4 +498,5 @@
 !--             be moved.
                 IF ( particles(n)%particle_mask )  THEN
+
                    j = ( particles(n)%y + 0.5_wp * dy ) * ddy
 !
@@ -521,12 +514,12 @@
                                particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
                                particles(n)%origin_y = ( ny + 1 ) * dy + &
-                               particles(n)%origin_y
+                                                     particles(n)%origin_y
                             ELSE
-                               trsp_count = trsp_count + 1
-                               trsp(trsp_count) = particles(n)
+                               trsp_count         = trsp_count + 1
+                               trsp(trsp_count)   = particles(n)
                                trsp(trsp_count)%y = ( ny + 1 ) * dy + &
-                               trsp(trsp_count)%y
+                                                 trsp(trsp_count)%y
                                trsp(trsp_count)%origin_y = trsp(trsp_count)%origin_y &
-                               + ( ny + 1 ) * dy
+                                                + ( ny + 1 ) * dy
                                particles(n)%particle_mask = .FALSE.
                                deleted_particles = deleted_particles + 1
@@ -536,5 +529,5 @@
                                   !++ why is 1 subtracted in next statement???
                                   trsp(trsp_count)%origin_y =                        &
-                                  trsp(trsp_count)%origin_y - 1
+                                                  trsp(trsp_count)%origin_y - 1
                                ENDIF
 
@@ -545,5 +538,5 @@
 !--                         Particle absorption
                             particles(n)%particle_mask = .FALSE.
-                            deleted_particles = deleted_particles + 1
+                            deleted_particles          = deleted_particles + 1
 
                          ELSEIF ( ibc_par_ns == 2 )  THEN
@@ -568,21 +561,21 @@
                       IF ( j > ny )  THEN
 !
-!--                       Apply boundary condition along x
+!--                       Apply boundary condition along y
                          IF ( ibc_par_ns == 0 )  THEN
 !
 !--                         Cyclic condition
                             IF ( pdims(2) == 1 )  THEN
-                               particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
+                               particles(n)%y        = particles(n)%y - ( ny + 1 ) * dy
                                particles(n)%origin_y =                         &
-                                  particles(n)%origin_y - ( ny + 1 ) * dy
+                                          particles(n)%origin_y - ( ny + 1 ) * dy
                             ELSE
-                               trnp_count = trnp_count + 1
-                               trnp(trnp_count) = particles(n)
+                               trnp_count         = trnp_count + 1
+                               trnp(trnp_count)   = particles(n)
                                trnp(trnp_count)%y =                            &
-                                  trnp(trnp_count)%y - ( ny + 1 ) * dy
+                                          trnp(trnp_count)%y - ( ny + 1 ) * dy
                                trnp(trnp_count)%origin_y =                     &
-                                  trnp(trnp_count)%origin_y - ( ny + 1 ) * dy
+                                         trnp(trnp_count)%origin_y - ( ny + 1 ) * dy
                                particles(n)%particle_mask = .FALSE.
-                               deleted_particles = deleted_particles + 1
+                               deleted_particles          = deleted_particles + 1
                             ENDIF
 
@@ -628,5 +621,5 @@
 
        ALLOCATE(rvnp(MAX(1,trnp_count_recv)))
-
+ 
        CALL MPI_SENDRECV( trsp(1)%radius, trsp_count, mpi_particle_type,      &
                           psouth, 1, rvnp(1)%radius,                             &
@@ -661,95 +654,119 @@
     ENDIF
 
+    DEALLOCATE( move_also_north )
+    DEALLOCATE( move_also_south )
+
 #else
 
-!
-!-- Apply boundary conditions
-    DO  n = 1, number_of_particles
-
-       IF ( particles(n)%x < -0.5_wp * dx )  THEN
-
-          IF ( ibc_par_lr == 0 )  THEN
-!
-!--          Cyclic boundary. Relevant coordinate has to be changed.
-             particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
-
-          ELSEIF ( ibc_par_lr == 1 )  THEN
-!
-!--          Particle absorption
-             particles(n)%particle_mask = .FALSE.
-             deleted_particles = deleted_particles + 1
-
-          ELSEIF ( ibc_par_lr == 2 )  THEN
-!
-!--          Particle reflection
-             particles(n)%x       = -dx - particles(n)%x
-             particles(n)%speed_x = -particles(n)%speed_x
-          ENDIF
-
-       ELSEIF ( particles(n)%x >= ( nx + 0.5_wp ) * dx )  THEN
-
-          IF ( ibc_par_lr == 0 )  THEN
-!
-!--          Cyclic boundary. Relevant coordinate has to be changed.
-             particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
-
-          ELSEIF ( ibc_par_lr == 1 )  THEN
-!
-!--          Particle absorption
-             particles(n)%particle_mask = .FALSE.
-             deleted_particles = deleted_particles + 1
-
-          ELSEIF ( ibc_par_lr == 2 )  THEN
-!
-!--          Particle reflection
-             particles(n)%x       = ( nx + 1 ) * dx - particles(n)%x
-             particles(n)%speed_x = -particles(n)%speed_x
-          ENDIF
-
-       ENDIF
-
-       IF ( particles(n)%y < -0.5_wp * dy )  THEN
-
-          IF ( ibc_par_ns == 0 )  THEN
-!
-!--          Cyclic boundary. Relevant coordinate has to be changed.
-             particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
-
-          ELSEIF ( ibc_par_ns == 1 )  THEN
-!
-!--          Particle absorption
-             particles(n)%particle_mask = .FALSE.
-             deleted_particles = deleted_particles + 1
-
-          ELSEIF ( ibc_par_ns == 2 )  THEN
-!
-!--          Particle reflection
-             particles(n)%y       = -dy - particles(n)%y
-             particles(n)%speed_y = -particles(n)%speed_y
-          ENDIF
-
-       ELSEIF ( particles(n)%y >= ( ny + 0.5_wp ) * dy )  THEN
-
-          IF ( ibc_par_ns == 0 )  THEN
-!
-!--          Cyclic boundary. Relevant coordinate has to be changed.
-             particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
-
-          ELSEIF ( ibc_par_ns == 1 )  THEN
-!
-!--          Particle absorption
-             particles(n)%particle_mask = .FALSE.
-             deleted_particles = deleted_particles + 1
-
-          ELSEIF ( ibc_par_ns == 2 )  THEN
-!
-!--          Particle reflection
-             particles(n)%y       = ( ny + 1 ) * dy - particles(n)%y
-             particles(n)%speed_y = -particles(n)%speed_y
-          ENDIF
-
-       ENDIF
+    DO  ip = nxl, nxr, nxr-nxl
+       DO  jp = nys, nyn
+          DO  kp = nzb+1, nzt
+             number_of_particles = prt_count(kp,jp,ip)
+             IF ( number_of_particles <= 0 )  CYCLE
+             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
+             DO  n = 1, number_of_particles
+!
+!--             Apply boundary conditions
+
+                IF ( particles(n)%x < -0.5_wp * dx )  THEN
+
+                   IF ( ibc_par_lr == 0 )  THEN
+!
+!--                   Cyclic boundary. Relevant coordinate has to be changed.
+                      particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
+
+                   ELSEIF ( ibc_par_lr == 1 )  THEN
+!
+!--                   Particle absorption
+                      particles(n)%particle_mask = .FALSE.
+                      deleted_particles = deleted_particles + 1
+
+                   ELSEIF ( ibc_par_lr == 2 )  THEN
+!
+!--                   Particle reflection
+                      particles(n)%x       = -dx - particles(n)%x
+                      particles(n)%speed_x = -particles(n)%speed_x
+                   ENDIF
+
+                ELSEIF ( particles(n)%x >= ( nx + 0.5_wp ) * dx )  THEN
+
+                   IF ( ibc_par_lr == 0 )  THEN
+!
+!--                   Cyclic boundary. Relevant coordinate has to be changed.
+                      particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
+
+                   ELSEIF ( ibc_par_lr == 1 )  THEN
+!
+!--                   Particle absorption
+                      particles(n)%particle_mask = .FALSE.
+                      deleted_particles = deleted_particles + 1
+
+                   ELSEIF ( ibc_par_lr == 2 )  THEN
+!
+!--                   Particle reflection
+                      particles(n)%x       = ( nx + 1 ) * dx - particles(n)%x
+                      particles(n)%speed_x = -particles(n)%speed_x
+                   ENDIF
+
+                ENDIF
+             ENDDO
+          ENDDO
+       ENDDO
     ENDDO
 
+    DO  ip = nxl, nxr
+       DO  jp = nys, nyn, nyn-nys
+          DO  kp = nzb+1, nzt
+             number_of_particles = prt_count(kp,jp,ip)
+             IF ( number_of_particles <= 0 )  CYCLE
+             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
+             DO  n = 1, number_of_particles
+
+                IF ( particles(n)%y < -0.5_wp * dy )  THEN
+
+                   IF ( ibc_par_ns == 0 )  THEN
+!
+!--                   Cyclic boundary. Relevant coordinate has to be changed.
+                      particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
+
+                   ELSEIF ( ibc_par_ns == 1 )  THEN
+!
+!--                   Particle absorption
+                      particles(n)%particle_mask = .FALSE.
+                      deleted_particles = deleted_particles + 1
+
+                   ELSEIF ( ibc_par_ns == 2 )  THEN
+!
+!--                   Particle reflection
+                      particles(n)%y       = -dy - particles(n)%y
+                      particles(n)%speed_y = -particles(n)%speed_y
+                   ENDIF
+
+                ELSEIF ( particles(n)%y >= ( ny + 0.5_wp ) * dy )  THEN
+
+                   IF ( ibc_par_ns == 0 )  THEN
+!
+!--                   Cyclic boundary. Relevant coordinate has to be changed.
+                      particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
+
+                   ELSEIF ( ibc_par_ns == 1 )  THEN
+!
+!--                   Particle absorption
+                      particles(n)%particle_mask = .FALSE.
+                      deleted_particles = deleted_particles + 1
+
+                   ELSEIF ( ibc_par_ns == 2 )  THEN
+!
+!--                   Particle reflection
+                      particles(n)%y       = ( ny + 1 ) * dy - particles(n)%y
+                      particles(n)%speed_y = -particles(n)%speed_y
+                   ENDIF
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+       ENDDO
+    ENDDO
 #endif
 
@@ -782,23 +799,24 @@
     IMPLICIT NONE
 
-    INTEGER(iwp)        ::  ip        !<
-    INTEGER(iwp)        ::  jp        !<
-    INTEGER(iwp)        ::  kp        !<
-    INTEGER(iwp)        ::  n         !<
-    INTEGER(iwp)        ::  pindex    !<
+    INTEGER(iwp)        ::  ip        !< grid index (x) of particle
+    INTEGER(iwp)        ::  jp        !< grid index (x) of particle
+    INTEGER(iwp)        ::  kp        !< grid index (x) of particle
+    INTEGER(iwp)        ::  n         !< index variable of particle
+    INTEGER(iwp)        ::  pindex    !< dummy argument for new number of particles per grid box
 
     LOGICAL             ::  pack_done !<
 
-    TYPE(particle_type), DIMENSION(:), INTENT(IN)       :: particle_array
+    TYPE(particle_type), DIMENSION(:), INTENT(IN)  ::  particle_array !< new particles in a grid box
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  temp_ns        !< temporary particle array for reallocation
 
     pack_done     = .FALSE.
 
-    nr_move_north = 0
-    nr_move_south = 0
-
     DO n = 1, SIZE(particle_array)
+
+       IF ( .NOT. particle_array(n)%particle_mask )  CYCLE
+
        ip = ( particle_array(n)%x + 0.5_wp * dx ) * ddx
        jp = ( particle_array(n)%y + 0.5_wp * dy ) * ddy
-       kp = particle_array(n)%z / dz + 1 + offset_ocean_nzt
+       kp =   particle_array(n)%z / dz + 1 + offset_ocean_nzt
 
        IF ( ip >= nxl  .AND.  ip <= nxr  .AND.  jp >= nys  .AND.  jp <= nyn    &
@@ -824,7 +842,25 @@
           prt_count(kp,jp,ip) = pindex
        ELSE
-          IF ( jp == nys - 1 )  THEN
+          IF ( jp <= nys - 1 )  THEN
              nr_move_south = nr_move_south+1
+!
+!--          Before particle information is swapped to exchange-array, check 
+!--          if enough memory is allocated. If required, reallocate exchange
+!--          array.
+             IF ( nr_move_south > SIZE(move_also_south) )  THEN
+!
+!--             At first, allocate further temporary array to swap particle 
+!--             information.
+                ALLOCATE( temp_ns(SIZE(move_also_south)+NR_2_direction_move) )
+                temp_ns(1:nr_move_south-1) = move_also_south(1:nr_move_south-1)
+                DEALLOCATE( move_also_south )
+                ALLOCATE( move_also_south(SIZE(temp_ns)) )
+                move_also_south(1:nr_move_south-1) = temp_ns(1:nr_move_south-1)
+                DEALLOCATE( temp_ns )
+
+             ENDIF
+
              move_also_south(nr_move_south) = particle_array(n)
+
              IF ( jp == -1 )  THEN
                 move_also_south(nr_move_south)%y =                             &
@@ -833,6 +869,23 @@
                    move_also_south(nr_move_south)%origin_y + ( ny + 1 ) * dy
              ENDIF
-          ELSEIF ( jp == nyn+1 )  THEN
+          ELSEIF ( jp >= nyn+1 )  THEN
              nr_move_north = nr_move_north+1
+!
+!--          Before particle information is swapped to exchange-array, check 
+!--          if enough memory is allocated. If required, reallocate exchange
+!--          array.
+             IF ( nr_move_north > SIZE(move_also_north) )  THEN
+!
+!--             At first, allocate further temporary array to swap particle 
+!--             information.
+                ALLOCATE( temp_ns(SIZE(move_also_north)+NR_2_direction_move) )
+                temp_ns(1:nr_move_north-1) = move_also_south(1:nr_move_north-1)
+                DEALLOCATE( move_also_north )
+                ALLOCATE( move_also_north(SIZE(temp_ns)) )
+                move_also_north(1:nr_move_north-1) = temp_ns(1:nr_move_north-1)
+                DEALLOCATE( temp_ns )
+
+             ENDIF
+
              move_also_north(nr_move_north) = particle_array(n)
              IF ( jp == ny+1 )  THEN
@@ -867,18 +920,18 @@
     IMPLICIT NONE
 
-    INTEGER(iwp)        ::  i           !<
-    INTEGER(iwp)        ::  ip          !<
-    INTEGER(iwp)        ::  j           !<
-    INTEGER(iwp)        ::  jp          !<
-    INTEGER(iwp)        ::  k           !<
-    INTEGER(iwp)        ::  kp          !<
-    INTEGER(iwp)        ::  n           !<
-    INTEGER(iwp)        ::  np_old_cell !<
-    INTEGER(iwp)        ::  n_start     !<
-    INTEGER(iwp)        ::  pindex      !<
+    INTEGER(iwp)        ::  i           !< grid index (x) of particle position
+    INTEGER(iwp)        ::  ip          !< index variable along x
+    INTEGER(iwp)        ::  j           !< grid index (y) of particle position
+    INTEGER(iwp)        ::  jp          !< index variable along y
+    INTEGER(iwp)        ::  k           !< grid index (z) of particle position
+    INTEGER(iwp)        ::  kp          !< index variable along z
+    INTEGER(iwp)        ::  n           !< index variable for particle array
+    INTEGER(iwp)        ::  np_old_cell !< number of particles per grid box before moving
+    INTEGER(iwp)        ::  n_start     !< start index
+    INTEGER(iwp)        ::  pindex      !< dummy argument for number of new particle per grid box
 
     LOGICAL             ::  pack_done   !<
 
-    TYPE(particle_type), DIMENSION(:), POINTER  ::  particles_old_cell !<
+    TYPE(particle_type), DIMENSION(:), POINTER  ::  particles_old_cell !< particles before moving
 
     CALL cpu_log( log_point_s(41), 'lpm_move_particle', 'start' )
@@ -947,10 +1000,7 @@
                             IF ( pack_done )  THEN
                                CALL realloc_particles_array(i,j,k)
-                               pindex = prt_count(k,j,i)+1
                             ELSE
                                CALL lpm_pack_arrays
                                prt_count(k,j,i) = number_of_particles
-
-                               pindex = prt_count(k,j,i)+1
 !
 !--                            If number of particles in the new grid box 
@@ -959,5 +1009,4 @@
                                IF ( pindex > SIZE(grid_particles(k,j,i)%particles) )  THEN
                                   CALL realloc_particles_array(i,j,k)
-                                  pindex = prt_count(k,j,i)+1
                                ENDIF
 
@@ -987,5 +1036,7 @@
 ! Description:
 ! ------------
-!> @todo Missing subroutine description.
+!> If the allocated memory for the particle array do not suffice to add arriving
+!> particles from neighbour grid cells, this subrouting reallocates the 
+!> particle array to assure enough memory is available. 
 !------------------------------------------------------------------------------!
  SUBROUTINE realloc_particles_array (i,j,k,size_in)
@@ -1012,5 +1063,5 @@
     ENDIF
 
-    new_size = MAX( new_size, min_nr_particle )
+    new_size = MAX( new_size, min_nr_particle, old_size + 1 )
 
     IF ( old_size <= 500 )  THEN
Index: /palm/trunk/SOURCE/lpm_init.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_init.f90	(revision 1928)
+++ /palm/trunk/SOURCE/lpm_init.f90	(revision 1929)
@@ -19,5 +19,10 @@
 ! Current revisions:
 ! -----------------
-!
+! Bugfix in determining initial particle height and grid index in case of 
+! seed_follows_topography.
+! Bugfix concerning random positions, ensure that particles do not move more 
+! than one grid length.
+! Bugfix logarithmic interpolation.
+! Initial setting of sgs_wf_part.
 ! 
 ! Former revisions:
@@ -31,5 +36,6 @@
 !
 ! 1873 2016-04-18 14:50:06Z maronga
-! Module renamed (removed _mod)
+! Module renamed (removed _mod
+
 ! 
 ! 1871 2016-04-15 11:46:09Z hoffmann
@@ -162,5 +168,5 @@
                 prt_count, psb, psl, psn, psr, pss, pst,                       &
                 radius, random_start_position, read_particles_from_restartfile,&
-                seed_follows_topography, sort_count,                           &
+                seed_follows_topography, sgs_wf_part, sort_count,              &
                 total_number_of_particles,                                     &
                 use_sgs_for_particles,                                         &
@@ -190,5 +196,5 @@
     PUBLIC lpm_init, lpm_create_particle
 
-CONTAINS
+ CONTAINS
 
 !------------------------------------------------------------------------------!
@@ -290,9 +296,12 @@
 
 !
-!-- Allocate arrays required for calculating particle SGS velocities
+!-- Allocate arrays required for calculating particle SGS velocities. 
+!-- Initialize prefactor required for stoachastic Weil equation.
     IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
        ALLOCATE( de_dx(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
                  de_dy(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
                  de_dz(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+
+       sgs_wf_part = 1.0_wp / 3.0_wp    
     ENDIF
 
@@ -334,5 +343,5 @@
 !
 !--    Precalculate LOG(z/z0)
-       height_p    = 0.0_wp
+       height_p    = z0_av_global
        DO  k = 1, number_of_sublayers
 
@@ -341,5 +350,4 @@
 
        ENDDO
-
 
     ENDIF
@@ -444,5 +452,5 @@
                                       0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
                                       0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0, 0, 0, &
-                                      0, .FALSE., -1)
+                                      0, .FALSE., -1 )
 
        particle_groups = particle_groups_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
@@ -477,5 +485,4 @@
 
        CALL lpm_create_particle (PHASE_INIT)
-
 !
 !--    User modification of initial particles
@@ -499,5 +506,4 @@
     number_of_particles = prt_count(nzb+1,nys,nxl)
     particles => grid_particles(nzb+1,nys,nxl)%particles(1:number_of_particles)
-
 !
 !-- Formats
@@ -520,30 +526,31 @@
 
     USE particle_attributes,                                                   &
-        ONLY: monodisperse_aerosols
+        ONLY: deleted_particles, monodisperse_aerosols
 
     IMPLICIT  NONE
 
-    INTEGER(iwp)               ::  alloc_size  !<
-    INTEGER(iwp)               ::  i           !<
-    INTEGER(iwp)               ::  ip          !<
-    INTEGER(iwp)               ::  j           !<
-    INTEGER(iwp)               ::  jp          !<
-    INTEGER(iwp)               ::  kp          !<
-    INTEGER(iwp)               ::  loop_stride !<
-    INTEGER(iwp)               ::  n           !<
-    INTEGER(iwp)               ::  new_size    !<
-
-    INTEGER(iwp), INTENT(IN)   ::  phase       !<
-
-    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  local_count !<
-    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  local_start !<
-
-    LOGICAL                    ::  first_stride !<
-
-    REAL(wp)                   ::  pos_x !<
-    REAL(wp)                   ::  pos_y !<
-    REAL(wp)                   ::  pos_z !<
-
-    TYPE(particle_type),TARGET ::  tmp_particle !<
+    INTEGER(iwp)               ::  alloc_size  !< relative increase of allocated memory for particles
+    INTEGER(iwp)               ::  i           !< loop variable ( particle groups )
+    INTEGER(iwp)               ::  ip          !< index variable along x
+    INTEGER(iwp)               ::  j           !< loop variable ( particles per point )
+    INTEGER(iwp)               ::  jp          !< index variable along y
+    INTEGER(iwp)               ::  kp          !< index variable along z
+    INTEGER(iwp)               ::  loop_stride !< loop variable for initialization
+    INTEGER(iwp)               ::  n           !< loop variable ( number of particles )
+    INTEGER(iwp)               ::  new_size    !< new size of allocated memory for particles
+
+    INTEGER(iwp), INTENT(IN)   ::  phase       !< mode of inititialization
+
+    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  local_count !< start address of new particle
+    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  local_start !< start address of new particle
+
+    LOGICAL                    ::  first_stride !< flag for initialization
+
+    REAL(wp)                   ::  pos_x      !< increment for particle position in x      
+    REAL(wp)                   ::  pos_y      !< increment for particle position in y  
+    REAL(wp)                   ::  pos_z      !< increment for particle position in z     
+    REAL(wp)                   ::  rand_contr !< dummy argument for random position
+
+    TYPE(particle_type),TARGET ::  tmp_particle !< temporary particle used for initialization
 
 !
@@ -618,4 +625,5 @@
                             tmp_particle%tail_id       = 0     !unused
 
+
 !
 !--                         Determine the grid indices of the particle position
@@ -628,9 +636,14 @@
 !--                            Particle height is given relative to topography
                                kp = kp + nzb_w_inner(jp,ip)
-                               tmp_particle%z = tmp_particle%z + zw(kp)
+                               tmp_particle%z = tmp_particle%z +               &
+                                                         zw(nzb_w_inner(jp,ip))
                                IF ( kp > nzt )  THEN
                                   pos_x = pos_x + pdx(i)
                                   CYCLE xloop
                                ENDIF
+                            ELSEIF ( .NOT. seed_follows_topography .AND.       &
+                                     tmp_particle%z <= zw(nzb_w_inner(jp,ip)) )  THEN
+                               pos_x = pos_x + pdx(i)
+                               CYCLE xloop                               
                             ENDIF
 
@@ -644,4 +657,5 @@
                                ENDIF
                                grid_particles(kp,jp,ip)%particles(local_count(kp,jp,ip)) = tmp_particle
+
                             ENDIF
                          ENDDO
@@ -691,8 +705,10 @@
                       ENDIF
                    ENDIF
+
                 ENDDO
              ENDDO
           ENDDO
        ENDIF
+
     ENDDO
 
@@ -707,5 +723,5 @@
 
 !
-!-- Add random fluctuation to particle positions
+!-- Add random fluctuation to particle positions.
     IF ( random_start_position )  THEN
        DO  ip = nxl, nxr
@@ -715,25 +731,54 @@
                 IF ( number_of_particles <= 0 )  CYCLE
                 particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
-
-                DO  n = local_start(kp,jp,ip), number_of_particles              !Move only new particles
+!
+!--             Move only new particles. Moreover, limit random fluctuation 
+!--             in order to prevent that particles move more than one grid box, 
+!--             which would lead to problems concerning particle exchange 
+!--             between processors in case pdx/pdy are larger than dx/dy, 
+!--             respectively.  
+                DO  n = local_start(kp,jp,ip), number_of_particles
                    IF ( psl(particles(n)%group) /= psr(particles(n)%group) )  THEN
+                      rand_contr = ( random_function( iran_part ) - 0.5_wp ) * &
+                                     pdx(particles(n)%group)
                       particles(n)%x = particles(n)%x +                        &
-                                   ( random_function( iran_part ) - 0.5_wp ) * &
-                                   pdx(particles(n)%group)
+                              MERGE( rand_contr, SIGN( dx, rand_contr ), &
+                                     ABS( rand_contr ) < dx                    &
+                                   ) 
                    ENDIF
                    IF ( pss(particles(n)%group) /= psn(particles(n)%group) )  THEN
+                      rand_contr = ( random_function( iran_part ) - 0.5_wp ) * &
+                                     pdy(particles(n)%group)
                       particles(n)%y = particles(n)%y +                        &
-                                   ( random_function( iran_part ) - 0.5_wp ) * &
-                                   pdy(particles(n)%group)
+                              MERGE( rand_contr, SIGN( dy, rand_contr ), &
+                                     ABS( rand_contr ) < dy                    &
+                                   ) 
                    ENDIF
                    IF ( psb(particles(n)%group) /= pst(particles(n)%group) )  THEN
+                      rand_contr = ( random_function( iran_part ) - 0.5_wp ) * &
+                                     pdz(particles(n)%group)
                       particles(n)%z = particles(n)%z +                        &
-                                   ( random_function( iran_part ) - 0.5_wp ) * &
-                                   pdz(particles(n)%group)
+                              MERGE( rand_contr, SIGN( dz, rand_contr ), &
+                                     ABS( rand_contr ) < dz                    &
+                                   ) 
                    ENDIF
                 ENDDO
 !
-!--             Identify particles located outside the model domain
+!--             Identify particles located outside the model domain and reflect
+!--             or absorb them if necessary.
                 CALL lpm_boundary_conds( 'bottom/top' )
+!
+!--             Furthermore, remove particles located in topography. Note, as
+!--             the particle speed is still zero at this point, wall 
+!--             reflection boundary conditions will not work in this case.
+                particles =>                                                   &
+                       grid_particles(kp,jp,ip)%particles(1:number_of_particles)
+                DO  n = local_start(kp,jp,ip), number_of_particles
+                   i = ( particles(n)%x + 0.5_wp * dx ) * ddx
+                   j = ( particles(n)%y + 0.5_wp * dy ) * ddy
+                   IF ( particles(n)%z <= zw(nzb_w_inner(j,i)) )  THEN
+                      particles(n)%particle_mask = .FALSE.
+                      deleted_particles = deleted_particles + 1
+                   ENDIF
+                ENDDO
              ENDDO
           ENDDO
Index: /palm/trunk/SOURCE/lpm_init_sgs_tke.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_init_sgs_tke.f90	(revision 1928)
+++ /palm/trunk/SOURCE/lpm_init_sgs_tke.f90	(revision 1929)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! ------------------
-!
+! sgs_wfu_par, sgs_wfv_par and sgs_wfw_par are replaced by sgs_wf_par
 ! 
 ! Former revisions:
@@ -72,5 +72,5 @@
 
     USE particle_attributes,                                                   &
-        ONLY:  sgs_wfu_part, sgs_wfv_part, sgs_wfw_part
+        ONLY:  sgs_wf_part
 
     USE pegrid
@@ -95,10 +95,10 @@
                   k  > nzb_s_inner(j,i+1) )                                    &
              THEN
-                de_dx(k,j,i) = 2.0_wp * sgs_wfu_part *                         &
+                de_dx(k,j,i) = 2.0_wp * sgs_wf_part *                          &
                                ( e(k,j,i+1) - e(k,j,i) ) * ddx
              ELSEIF ( k  > nzb_s_inner(j,i-1)  .AND.  k > nzb_s_inner(j,i)     &
                       .AND.  k <= nzb_s_inner(j,i+1) )                         &
              THEN
-                de_dx(k,j,i) = 2.0_wp * sgs_wfu_part *                         &
+                de_dx(k,j,i) = 2.0_wp * sgs_wf_part *                          &
                                ( e(k,j,i) - e(k,j,i-1) ) * ddx
              ELSEIF ( k < nzb_s_inner(j,i)  .AND.  k < nzb_s_inner(j,i+1) )    &
@@ -109,5 +109,5 @@
                 de_dx(k,j,i) = 0.0_wp
              ELSE
-                de_dx(k,j,i) = sgs_wfu_part * ( e(k,j,i+1) - e(k,j,i-1) ) * ddx
+                de_dx(k,j,i) = sgs_wf_part * ( e(k,j,i+1) - e(k,j,i-1) ) * ddx
              ENDIF
 
@@ -115,10 +115,10 @@
                   k  > nzb_s_inner(j+1,i) )                                    &
              THEN
-                de_dy(k,j,i) = 2.0_wp * sgs_wfv_part *                         &
+                de_dy(k,j,i) = 2.0_wp * sgs_wf_part *                          &
                                ( e(k,j+1,i) - e(k,j,i) ) * ddy
              ELSEIF ( k  > nzb_s_inner(j-1,i)  .AND.  k  > nzb_s_inner(j,i)    &
                       .AND.  k <= nzb_s_inner(j+1,i) )                         &
              THEN
-                de_dy(k,j,i) = 2.0_wp * sgs_wfv_part *                         &
+                de_dy(k,j,i) = 2.0_wp * sgs_wf_part *                          &
                                ( e(k,j,i) - e(k,j-1,i) ) * ddy
              ELSEIF ( k < nzb_s_inner(j,i)  .AND.  k < nzb_s_inner(j+1,i) )    &
@@ -129,5 +129,5 @@
                 de_dy(k,j,i) = 0.0_wp
              ELSE
-                de_dy(k,j,i) = sgs_wfv_part * ( e(k,j+1,i) - e(k,j-1,i) ) * ddy
+                de_dy(k,j,i) = sgs_wf_part * ( e(k,j+1,i) - e(k,j-1,i) ) * ddy
              ENDIF
 
@@ -142,5 +142,5 @@
 
           DO  k = nzb_s_inner(j,i)+2, nzt-1
-             de_dz(k,j,i)  = 2.0_wp * sgs_wfw_part *                           &
+             de_dz(k,j,i)  = 2.0_wp * sgs_wf_part *                            &
                              ( e(k+1,j,i) - e(k-1,j,i) ) / ( zu(k+1)-zu(k-1) )
           ENDDO
@@ -148,5 +148,5 @@
           k = nzb_s_inner(j,i)
           de_dz(nzb:k,j,i) = 0.0_wp
-          de_dz(k+1,j,i)   = 2.0_wp * sgs_wfw_part *                           &
+          de_dz(k+1,j,i)   = 2.0_wp * sgs_wf_part *                            &
                              ( e(k+2,j,i) - e(k+1,j,i) ) / ( zu(k+2) - zu(k+1) )
           de_dz(nzt,j,i)   = 0.0_wp
Index: /palm/trunk/SOURCE/mod_particle_attributes.f90
===================================================================
--- /palm/trunk/SOURCE/mod_particle_attributes.f90	(revision 1928)
+++ /palm/trunk/SOURCE/mod_particle_attributes.f90	(revision 1929)
@@ -19,5 +19,6 @@
 ! Current revisions:
 ! ------------------
-!
+! -sgs_wfu_par, sgs_wfv_par, sgs_wfw_par
+! + sgs_wf_par
 !
 ! Former revisions:
@@ -109,7 +110,5 @@
                  initial_weighting_factor = 1.0_wp,                            &
                  particle_advection_start = 0.0_wp,                            &
-                 sgs_wfu_part = 0.3333333_wp, sgs_wfv_part = 0.3333333_wp,     &
-                 sgs_wfw_part = 0.3333333_wp, time_prel = 0.0_wp,              &
-                 time_sort_particles = 0.0_wp,                                 &
+                 sgs_wf_part, time_prel = 0.0_wp, time_sort_particles = 0.0_wp,&
                  time_write_particle_data = 0.0_wp, z0_av_global
 
Index: /palm/trunk/SOURCE/pres.f90
===================================================================
--- /palm/trunk/SOURCE/pres.f90	(revision 1928)
+++ /palm/trunk/SOURCE/pres.f90	(revision 1929)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! ------------------
-!
+! Bugfix: weight_substep for initial call, replace by local variable
 !
 ! Former revisions:
@@ -162,4 +162,5 @@
     REAL(wp)     ::  threadsum      !<
     REAL(wp)     ::  weight_pres_l  !<
+    REAL(wp)     ::  weight_substep_l !<
 
     REAL(wp), DIMENSION(1:3)   ::  volume_flow_l       !<
@@ -178,9 +179,11 @@
 !
 !--    If pres is called before initial time step
-       weight_pres_l = 1.0_wp
-       d_weight_pres = 1.0_wp
+       weight_pres_l    = 1.0_wp
+       d_weight_pres    = 1.0_wp
+       weight_substep_l = 1.0_wp
     ELSE
-       weight_pres_l = weight_pres(intermediate_timestep_count)
-       d_weight_pres = 1.0_wp / weight_pres(intermediate_timestep_count)
+       weight_pres_l    = weight_pres(intermediate_timestep_count)
+       d_weight_pres    = 1.0_wp / weight_pres(intermediate_timestep_count)
+       weight_substep_l = weight_substep(intermediate_timestep_count)
     ENDIF
 
@@ -591,5 +594,5 @@
        !$OMP PARALLEL PRIVATE (i,j,k)
        !$OMP DO
-       !$acc kernels present( p, tend, weight_substep )
+       !$acc kernels present( p, tend, weight_substep_l )
        !$acc loop independent
        DO  i = nxl-1, nxr+1
@@ -599,5 +602,5 @@
              DO  k = nzb, nzt+1
                 p(k,j,i) = tend(k,j,i) * &
-                           weight_substep(intermediate_timestep_count)
+                           weight_substep_l
              ENDDO
           ENDDO
@@ -609,5 +612,5 @@
        !$OMP PARALLEL PRIVATE (i,j,k)
        !$OMP DO
-       !$acc kernels present( p, tend, weight_substep )
+       !$acc kernels present( p, tend, weight_substep_l )
        !$acc loop independent
        DO  i = nxl-1, nxr+1
@@ -617,5 +620,5 @@
              DO  k = nzb, nzt+1
                 p(k,j,i) = p(k,j,i) + tend(k,j,i) * &
-                           weight_substep(intermediate_timestep_count)
+                           weight_substep_l
              ENDDO
           ENDDO
Index: /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
===================================================================
--- /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 1928)
+++ /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 1929)
@@ -19,5 +19,6 @@
 ! Current revisions:
 ! ------------------
-! 
+! Bugfix: avoid segmentation fault in case one grid point is horizontally 
+! completely surrounded by topography 
 ! 
 ! Former revisions:
@@ -726,6 +727,8 @@
 !--             consequence would result in very large shear stresses and very
 !--             small momentum fluxes (both are generally unrealistic).
-                IF ( ( z_mo / ol(j,i) ) < zeta_min )  ol(j,i) = z_mo / zeta_min
-                IF ( ( z_mo / ol(j,i) ) > zeta_max )  ol(j,i) = z_mo / zeta_max
+                IF ( ( z_mo / ( ol(j,i) + 1E-30_wp ) ) < zeta_min )            &
+                   ol(j,i) = z_mo / zeta_min
+                IF ( ( z_mo / ( ol(j,i) + 1E-30_wp ) ) > zeta_max )            &
+                   ol(j,i) = z_mo / zeta_max
              ENDDO
           ENDDO
Index: /palm/trunk/SOURCE/wind_turbine_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/wind_turbine_model_mod.f90	(revision 1928)
+++ /palm/trunk/SOURCE/wind_turbine_model_mod.f90	(revision 1929)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix: added preprocessor directives for parallel and serial mode
 ! 
 ! Former revisions:
@@ -2058,5 +2058,6 @@
 
 !
-!--          Transfer of information to the other cpus            
+!--          Transfer of information to the other cpus
+#if defined( __parallel )         
              CALL MPI_ALLREDUCE( u_inflow_l, u_inflow, nturbines, MPI_REAL,    &
                                  MPI_SUM, comm2d, ierr )
@@ -2065,5 +2066,9 @@
              CALL MPI_ALLREDUCE( phi_yaw_l, phi_yaw, nturbines, MPI_REAL,      &
                                  MPI_SUM, comm2d, ierr )
-             
+#else
+             u_inflow = u_inflow_l
+             wdir     = wdir_l
+             phi_yaw  = phi_yaw_l
+#endif
              DO inot = 1, nturbines
 !             
@@ -2080,4 +2085,5 @@
 !              CALL MPI_ALLREDUCE( omega_gen, omega_gen_old, nturbines,        &
 !                                  MPI_REAL,MPI_SUM, comm2d, ierr )
+#if defined( __parallel )    
              CALL MPI_ALLREDUCE( torque_gen, torque_gen_old, nturbines,        &
                                  MPI_REAL, MPI_SUM, comm2d, ierr )
@@ -2086,4 +2092,9 @@
              CALL MPI_ALLREDUCE( omega_gen_f, omega_gen_f_old, nturbines,      &
                                  MPI_REAL, MPI_SUM, comm2d, ierr )
+#else
+             torque_gen_old  = torque_gen
+             omega_rot       = omega_rot_l
+             omega_gen_f_old = omega_gen_f
+#endif
             
           ENDIF
