Index: /palm/trunk/SOURCE/Makefile
===================================================================
--- /palm/trunk/SOURCE/Makefile	(revision 57)
+++ /palm/trunk/SOURCE/Makefile	(revision 58)
@@ -6,5 +6,5 @@
 # default suffixes removed from the suffix list to avoid calling of m2c in
 # case of .mod files
-# +wall_fluxes
+# +wall_fluxes, +particle_boundary_conds
 #
 # Former revisions:
@@ -49,13 +49,14 @@
 	local_getenv.f90 local_stop.f90 local_system.f90 local_tremain.f90 \
 	local_tremain_ini.f90 modules.f90 netcdf.f90 package_parin.f90 \
-	palm.f90 parin.f90 poisfft.f90 poisfft_hybrid.f90 poismg.f90 \
-	prandtl_fluxes.f90 pres.f90 print_1d.f90 production_e.f90 \
-	prognostic_equations.f90 random_function.f90 random_gauss.f90 \
-	read_3d_binary.f90 read_var_list.f90 run_control.f90 singleton.f90 \
-	sor.f90 spline_x.f90 spline_y.f90 spline_z.f90 sum_up_3d_data.f90 \
-	swap_timelevel.f90 temperton_fft.f90 time_integration.f90 \
-	time_to_string.f90 timestep.f90 timestep_scheme_steering.f90 \
-	transpose.f90 user_interface.f90 wall_fluxes.f90 write_3d_binary.f90 \
-	write_compressed.f90 write_var_list.f90
+	palm.f90 parin.f90 particle_boundary_conds.f90 poisfft.f90 \
+	poisfft_hybrid.f90 poismg.f90 prandtl_fluxes.f90 pres.f90 print_1d.f90 \
+	production_e.f90 prognostic_equations.f90 random_function.f90 \
+	random_gauss.f90 read_3d_binary.f90 read_var_list.f90 run_control.f90 \
+	singleton.f90 sor.f90 spline_x.f90 spline_y.f90 spline_z.f90 \
+	sum_up_3d_data.f90 swap_timelevel.f90 temperton_fft.f90 \
+	time_integration.f90 time_to_string.f90 timestep.f90 \
+	timestep_scheme_steering.f90 transpose.f90 user_interface.f90 \
+	wall_fluxes.f90 write_3d_binary.f90 write_compressed.f90 \
+	write_var_list.f90
 
 OBJS =  advec_particles.o advec_s_bc.o advec_s_pw.o advec_s_up.o \
@@ -78,7 +79,7 @@
 	interaction_droplets_ptq.o local_getenv.o local_stop.o \
 	local_system.o local_tremain.o local_tremain_ini.o modules.o netcdf.o \
-	package_parin.o palm.o parin.o poisfft.o poisfft_hybrid.o poismg.o \
-	prandtl_fluxes.o pres.o print_1d.o production_e.o \
-	prognostic_equations.o random_function.o random_gauss.o \
+	package_parin.o palm.o parin.o particle_boundary_conds.o poisfft.o \
+	poisfft_hybrid.o poismg.o prandtl_fluxes.o pres.o print_1d.o \
+	production_e.o prognostic_equations.o random_function.o random_gauss.o \
 	read_3d_binary.o read_var_list.o run_control.o singleton.o sor.o \
 	spline_x.o spline_y.o spline_z.o sum_up_3d_data.o swap_timelevel.o \
@@ -187,4 +188,5 @@
 palm.o: modules.o
 parin.o: modules.o
+particle_boundary_conds.o: modules.o
 poisfft.o: modules.o fft_xy.o
 poisfft_hybrid.o: modules.o fft_xy.o
Index: /palm/trunk/SOURCE/particle_boundary_conds.f90
===================================================================
--- /palm/trunk/SOURCE/particle_boundary_conds.f90	(revision 58)
+++ /palm/trunk/SOURCE/particle_boundary_conds.f90	(revision 58)
@@ -0,0 +1,599 @@
+ SUBROUTINE particle_boundary_conds
+
+!------------------------------------------------------------------------------!
+! Actual revisions:
+! -----------------
+! 
+!
+! Former revisions:
+! -----------------
+! $Id$
+! Initial version (2007/03/09)
+!
+! Description:
+! ------------
+! Calculates the reflection of particles from vertical walls.
+! Routine developed by Jin Zhang (2006-2007)
+!------------------------------------------------------------------------------!
+#if defined( __particles )
+
+    USE control_parameters
+    USE cpulog
+    USE grid_variables
+    USE indices
+    USE interfaces
+    USE particle_attributes
+    USE pegrid
+
+    IMPLICIT NONE
+
+    INTEGER ::  n
+
+
+    CALL cpu_log( log_point_s(48), 'advec_part_refle', 'start' )
+
+
+
+    DO  n = 1, number_of_particles
+
+       i2 = ( particles(n)%x + 0.5 * dx ) * ddx
+       j2 = ( particles(n)%y + 0.5 * dy ) * ddy
+       k2 = particles(n)%z / dz + 1
+
+       prt_x = particles(n)%x
+       prt_y = particles(n)%y
+       prt_z = particles(n)%z
+
+
+       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 * dx ) * ddx
+          j1 = ( pos_y_old + 0.5 * dy ) * ddy
+          k1 = pos_z_old / dz
+
+!
+!--       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 * dx
+                t(t_index) = ( xline - pos_x_old ) / &
+                             ( particles(n)%x - pos_x_old )
+                t_index    = t_index + 1
+             ENDDO
+
+             DO  j = j1, j2
+                yline      = j * dy + 0.5 * 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
+
+             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 * dx ) * ddx   
+                j3 = ( pos_y + 0.5 * dy ) * ddy
+                k3 = pos_z / dz
+
+                i5 = pos_x * ddx
+                j5 = pos_y * ddy
+                k5 = pos_z / dz
+
+                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
+                      particles(n)%z       = 2.0 * pos_z - prt_z
+                      particles(n)%speed_z = - particles(n)%speed_z
+
+                      IF ( use_sgs_for_particles )  THEN
+                         particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+                      ENDIF
+                      GOTO 999
+                   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
+                      particles(n)%z       = 2.0 * pos_z - prt_z
+                      particles(n)%speed_z = - particles(n)%speed_z
+
+                      IF ( use_sgs_for_particles )  THEN
+                         particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+                      ENDIF
+                      GOTO 999
+                   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
+                      particles(n)%z       = 2.0 * pos_z - prt_z
+                      particles(n)%speed_z = - particles(n)%speed_z
+
+                      IF ( use_sgs_for_particles )  THEN
+                         particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+                      ENDIF
+                      GOTO 999
+                   ENDIF
+
+                   IF ( pos_y == ( j3 * dy - 0.5 * dy )  .AND. &
+                        pos_z < nzb_s_inner(j3,i3) * dz )  THEN
+                      particles(n)%y       = 2.0 * pos_y - prt_y
+                      particles(n)%speed_y = - particles(n)%speed_y
+
+                      IF ( use_sgs_for_particles )  THEN
+                         particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs
+                      ENDIF
+                      GOTO 999
+                   ENDIF
+
+                   IF ( pos_x == ( i3 * dx - 0.5 * dx )  .AND. &
+                        pos_z < nzb_s_inner(j3,i3) * dz )  THEN
+                      particles(n)%x       = 2.0 * pos_x - prt_x
+                      particles(n)%speed_x = - particles(n)%speed_x
+
+                      IF ( use_sgs_for_particles )  THEN 
+                         particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs
+                      ENDIF
+                      GOTO 999
+                   ENDIF
+                ENDIF
+
+             ENDDO
+
+          ENDIF
+
+!
+!--       Case 2
+          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 * 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 * 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
+
+             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 * dx ) * ddx
+                j3 = ( pos_y + 0.5 * dy ) * ddy
+                k3 = pos_z / dz
+
+                i5 = pos_x * ddx
+                j5 = pos_y * ddy
+                k5 = pos_z / dz
+
+                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
+                      particles(n)%z       = 2.0 * pos_z - prt_z
+                      particles(n)%speed_z = - particles(n)%speed_z
+
+                      IF ( use_sgs_for_particles )  THEN
+                         particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+                      ENDIF
+                      GOTO 999
+                   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
+                      particles(n)%z       = 2.0 * pos_z - prt_z
+                      particles(n)%speed_z = - particles(n)%speed_z
+
+                      IF ( use_sgs_for_particles )  THEN
+                         particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+                      ENDIF
+                      GOTO 999
+                   ENDIF
+
+                   IF ( pos_x == ( i3 * dx - 0.5 * dx )  .AND. &
+                        pos_z < nzb_s_inner(j3,i3) * dz )  THEN
+                      particles(n)%x       = 2.0 * pos_x - prt_x
+                      particles(n)%speed_x = - particles(n)%speed_x
+
+                      IF ( use_sgs_for_particles )  THEN
+                         particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs
+                      ENDIF
+                      GOTO 999
+                   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
+                      particles(n)%z       = 2.0 * pos_z - prt_z
+                      particles(n)%speed_z = - particles(n)%speed_z
+                      IF ( use_sgs_for_particles )  THEN
+                         particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+                      ENDIF
+                      GOTO 999
+                   ENDIF
+
+                   IF ( pos_y == ( j5 * dy + 0.5 * dy )  .AND. &
+                        pos_z < nzb_s_inner(j5,i3) * dz )  THEN
+                      particles(n)%y       = 2.0 * pos_y - prt_y
+                      particles(n)%speed_y = - particles(n)%speed_y
+
+                      IF ( use_sgs_for_particles )  THEN
+                         particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs
+                      ENDIF
+                      GOTO 999
+                   ENDIF
+                ENDIF
+             ENDDO
+          ENDIF
+
+!
+!--       Case 3
+          IF ( 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 * 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 * 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
+
+             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 * dx ) * ddx
+                j3 = ( pos_y + 0.5 * dy ) * ddy
+                k3 =  pos_z / dz
+
+                i5 = pos_x * ddx
+                j5 = pos_y * ddy
+                k5 = pos_z / dz
+
+                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
+                      particles(n)%z       = 2.0 * pos_z - prt_z
+                      particles(n)%speed_z = - particles(n)%speed_z
+
+                      IF ( use_sgs_for_particles )  THEN  
+                         particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+                      ENDIF
+                      GOTO 999
+                   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
+                      particles(n)%z       = 2.0 * pos_z - prt_z
+                      particles(n)%speed_z = - particles(n)%speed_z
+
+                      IF ( use_sgs_for_particles )  THEN
+                         particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+                      ENDIF
+                      GOTO 999
+                   ENDIF
+
+                   IF ( pos_y == ( j3 * dy - 0.5 * dy )  .AND. &
+                        pos_z < nzb_s_inner(j3,i3) * dz )  THEN
+                      particles(n)%y       = 2.0 * pos_y - prt_y
+                      particles(n)%speed_y = - particles(n)%speed_y
+
+                      IF ( use_sgs_for_particles )  THEN
+                         particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs
+                      ENDIF
+                      GOTO 999
+                   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
+                      particles(n)%z       = 2.0 * pos_z - prt_z
+                      particles(n)%speed_z = - particles(n)%speed_z
+
+                      IF ( use_sgs_for_particles )  THEN
+                         particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+                      ENDIF
+                      GOTO 999
+                   ENDIF
+
+                   IF ( pos_x == ( i5 * dx + 0.5 * dx )  .AND. &
+                        pos_z < nzb_s_inner(j3,i5) * dz )  THEN
+                      particles(n)%x       = 2.0 * pos_x - prt_x
+                      particles(n)%speed_x = - particles(n)%speed_x
+
+                      IF ( use_sgs_for_particles )  THEN
+                         particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs
+                      ENDIF
+                      GOTO 999
+                   ENDIF
+                ENDIF 
+             ENDDO
+          ENDIF
+
+!
+!--       Case 4
+          IF ( 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 * 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 * 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
+
+             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 * dx ) * ddx   
+                j3 = ( pos_y + 0.5 * dy ) * ddy
+                k3 = pos_z / dz
+
+                i5 = pos_x * ddx
+                j5 = pos_y * ddy
+                k5 = pos_z / dz
+
+                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
+                      particles(n)%z       = 2.0 * pos_z - prt_z
+                      particles(n)%speed_z = - particles(n)%speed_z
+
+                      IF ( use_sgs_for_particles )  THEN
+                         particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+                      ENDIF
+                      GOTO 999
+                   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
+                      particles(n)%z       = 2.0 * pos_z - prt_z
+                      particles(n)%speed_z = - particles(n)%speed_z
+
+                      IF ( use_sgs_for_particles )  THEN 
+                         particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+                      ENDIF
+                      GOTO 999
+                   ENDIF
+
+                   IF ( pos_x == ( i5 * dx + 0.5 * dx )  .AND. &
+                        nzb_s_inner(j3,i5) /=0  .AND.          &
+                        pos_z < nzb_s_inner(j3,i5) * dz )  THEN
+                      particles(n)%x       = 2.0 * pos_x - prt_x
+                      particles(n)%speed_x = - particles(n)%speed_x
+
+                      IF ( use_sgs_for_particles )  THEN  
+                         particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs
+                      ENDIF
+                      GOTO 999
+                   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
+                      particles(n)%z       = 2.0 * pos_z - prt_z
+                      particles(n)%speed_z = - particles(n)%speed_z
+
+                      IF ( use_sgs_for_particles )  THEN 
+                         particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+                      ENDIF
+                      GOTO 999
+                   ENDIF
+
+                   IF ( pos_y == ( j5 * dy + 0.5 * dy )  .AND. &
+                        nzb_s_inner(j5,i3) /= 0  .AND.         &
+                        pos_z < nzb_s_inner(j5,i3) * dz )  THEN
+                      particles(n)%y       = 2.0 * pos_y - prt_y
+                      particles(n)%speed_y = - particles(n)%speed_y
+
+                      IF ( use_sgs_for_particles )  THEN
+                         particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs
+                      ENDIF
+                      GOTO 999
+                   ENDIF
+                ENDIF 
+             ENDDO
+          ENDIF
+
+     999  CONTINUE
+    
+       ENDIF
+
+    ENDDO
+
+    CALL cpu_log( log_point_s(48), 'advec_part_refle', 'stop' )
+
+ END SUBROUTINE particle_boundary_conds
