Index: palm/trunk/SOURCE/lpm_boundary_conds.f90
===================================================================
--- palm/trunk/SOURCE/particle_boundary_conds.f90	(revision 826)
+++ palm/trunk/SOURCE/lpm_boundary_conds.f90	(revision 849)
@@ -1,8 +1,9 @@
- SUBROUTINE particle_boundary_conds
+ SUBROUTINE lpm_boundary_conds( range )
 
 !------------------------------------------------------------------------------!
 ! Current revisions:
 ! -----------------
-! 
+! routine renamed lpm_boundary_conds, bottom and top boundary conditions
+! included
 !
 ! Former revisions:
@@ -20,12 +21,20 @@
 ! Description:
 ! ------------
-! Calculates the reflection of particles from vertical walls.
-! Routine developed by Jin Zhang (2006-2007)
+! Boundary conditions for the Lagrangian particles.
+! The routine consists of two different parts. One handles the bottom (flat)
+! and top boundary. In this part, also particles which exceeded their lifetime
+! are deleted.
+! The other part handles the reflection of particles from vertical walls.
+! This part was developed by Jin Zhang during 2006-2007.
+!
 ! To do: Code structure for finding the t_index values and for checking the
-!        reflection conditions is basically the same for all four cases, so it
+! -----  reflection conditions is basically the same for all four cases, so it
 !        should be possible to further simplify/shorten it.
-! THIS ROUTINE HAS NOT BEEN TESTED FOR OCEAN RUNS SO FAR! (see offset_ocean_*)
+!
+! THE WALLS PART OF THIS ROUTINE HAS NOT BEEN TESTED FOR OCEAN RUNS SO FAR!!!!
+! (see offset_ocean_*)
 !------------------------------------------------------------------------------!
 
+    USE arrays_3d
     USE control_parameters
     USE cpulog
@@ -38,6 +47,8 @@
     IMPLICIT NONE
 
+    CHARACTER (LEN=*) ::  range
+
     INTEGER ::  i, inc, ir, i1, i2, i3, i5, j, jr, j1, j2, j3, j5, k, k1, k2, &
-                k3, k5, n, t_index, t_index_number
+                k3, k5, n, nn, t_index, t_index_number
 
     LOGICAL ::  reflect_x, reflect_y, reflect_z
@@ -48,538 +59,631 @@
     REAL ::  t(1:200)
 
-    CALL cpu_log( log_point_s(48), 'advec_part_refle', 'start' )
-
-
-
-    reflect_x = .FALSE.
-    reflect_y = .FALSE.
-    reflect_z = .FALSE.
-
-    DO  n = 1, number_of_particles
-
-       dt_particle = particles(n)%age - particles(n)%age_m
-
-       i2 = ( particles(n)%x + 0.5 * dx ) * ddx
-       j2 = ( particles(n)%y + 0.5 * dy ) * ddy
-       k2 = particles(n)%z / dz + 1 + offset_ocean_nzt_m1
-
-       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 * dx ) * ddx
-          j1 = ( pos_y_old + 0.5 * 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 ) &
+
+
+    IF ( range == 'bottom/top' )  THEN
+
+!
+!--    Apply boundary conditions to those particles that have crossed the top or
+!--    bottom boundary and delete those particles, which are older than allowed
+       DO  n = 1, number_of_particles
+
+          nn = particles(n)%tail_id
+
+!
+!--       Stop if particles have moved further than the length of one 
+!--       PE subdomain (newly released particles have age = age_m!)
+          IF ( particles(n)%age /= particles(n)%age_m )  THEN
+             IF ( ABS(particles(n)%speed_x) >                                  &
+                  ((nxr-nxl+2)*dx)/(particles(n)%age-particles(n)%age_m)  .OR. &
+                  ABS(particles(n)%speed_y) >                                  &
+                  ((nyn-nys+2)*dy)/(particles(n)%age-particles(n)%age_m) )  THEN
+
+                  WRITE( message_string, * )  'particle too fast.  n = ',  n 
+                  CALL message( 'lpm_boundary_conds', 'PA0148', 2, 2, -1, 6, 1 )
+             ENDIF
+          ENDIF
+
+          IF ( particles(n)%age > particle_maximum_age  .AND.  &
+               particle_mask(n) )                              &
           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
-
-      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 * dx ) * ddx   
-                j3 = ( pos_y + 0.5 * 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
-
-                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 * dy )  .AND. &
-                        pos_z < nzb_s_inner(j3,i3) * dz )  THEN
-                      reflect_y = .TRUE.
-                      EXIT case1
-                   ENDIF
-
-                   IF ( pos_x == ( i3 * dx - 0.5 * 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 * 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
-
-      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 * dx ) * ddx
-                j3 = ( pos_y + 0.5 * 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 * 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 * 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 * 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
-
-      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 * dx ) * ddx
-                j3 = ( pos_y + 0.5 * 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 * 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 * 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 * 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
-
-      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 * dx ) * ddx   
-                j3 = ( pos_y + 0.5 * 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 * 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 * 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
+             particle_mask(n)  = .FALSE.
+             deleted_particles = deleted_particles + 1
+             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
+                tail_mask(nn) = .FALSE.
+                deleted_tails = deleted_tails + 1
+             ENDIF
+          ENDIF
+
+          IF ( particles(n)%z >= zu(nz)  .AND.  particle_mask(n) )  THEN
+             IF ( ibc_par_t == 1 )  THEN
+!
+!--             Particle absorption
+                particle_mask(n)  = .FALSE.
+                deleted_particles = deleted_particles + 1
+                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
+                   tail_mask(nn) = .FALSE.
+                   deleted_tails = deleted_tails + 1
+                ENDIF
+             ELSEIF ( ibc_par_t == 2 )  THEN
+!
+!--             Particle reflection
+                particles(n)%z       = 2.0 * zu(nz) - particles(n)%z
+                particles(n)%speed_z = -particles(n)%speed_z
+                IF ( use_sgs_for_particles  .AND. &
+                     particles(n)%rvar3 > 0.0 )  THEN
+                   particles(n)%rvar3 = -particles(n)%rvar3
+                ENDIF
+                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
+                   particle_tail_coordinates(1,3,nn) = 2.0 * zu(nz) - &
+                                               particle_tail_coordinates(1,3,nn)
+                ENDIF
+             ENDIF
+          ENDIF
+          IF ( particles(n)%z < zw(0)  .AND.  particle_mask(n) )  THEN
+             IF ( ibc_par_b == 1 )  THEN
+!
+!--             Particle absorption
+                particle_mask(n)  = .FALSE.
+                deleted_particles = deleted_particles + 1
+                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
+                   tail_mask(nn) = .FALSE.
+                   deleted_tails = deleted_tails + 1
+                ENDIF
+             ELSEIF ( ibc_par_b == 2 )  THEN
+!
+!--             Particle reflection
+                particles(n)%z       = 2.0 * zw(0) - particles(n)%z
+                particles(n)%speed_z = -particles(n)%speed_z
+                IF ( use_sgs_for_particles  .AND. &
+                     particles(n)%rvar3 < 0.0 )  THEN
+                   particles(n)%rvar3 = -particles(n)%rvar3
+                ENDIF
+                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
+                   particle_tail_coordinates(1,3,nn) = 2.0 * zu(nz) - &
+                                               particle_tail_coordinates(1,3,nn)
+                ENDIF
+                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
+                   particle_tail_coordinates(1,3,nn) = 2.0 * zw(0) - &
+                                               particle_tail_coordinates(1,3,nn)
+                ENDIF
+             ENDIF
+          ENDIF
+       ENDDO
+
+    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
+
+          dt_particle = particles(n)%age - particles(n)%age_m
+
+          i2 = ( particles(n)%x + 0.5 * dx ) * ddx
+          j2 = ( particles(n)%y + 0.5 * dy ) * ddy
+          k2 = particles(n)%z / dz + 1 + offset_ocean_nzt_m1
+
+          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 * dx ) * ddx
+             j1 = ( pos_y_old + 0.5 * 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 * 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
+
+         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 * dx ) * ddx   
+                   j3 = ( pos_y + 0.5 * 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
+
+                   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 * dy )  .AND. &
+                           pos_z < nzb_s_inner(j3,i3) * dz )  THEN
+                         reflect_y = .TRUE.
+                         EXIT case1
+                      ENDIF
+
+                      IF ( pos_x == ( i3 * dx - 0.5 * 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 * 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
+
+         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 * dx ) * ddx
+                   j3 = ( pos_y + 0.5 * 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 * 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 * 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 * 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
+
+         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 * dx ) * ddx
+                   j3 = ( pos_y + 0.5 * 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 * 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 * 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 * 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
+
+         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 * dx ) * ddx   
+                   j3 = ( pos_y + 0.5 * 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 * 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 * 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
-
-          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 * 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 * 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 * 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        
+                ENDDO case4
+
+             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 * 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 * 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 * 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        
     
-    ENDDO
-
-    CALL cpu_log( log_point_s(48), 'advec_part_refle', 'stop' )
-
-
- END SUBROUTINE particle_boundary_conds
+       ENDDO
+
+       CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'stop' )
+
+    ENDIF
+
+ END SUBROUTINE lpm_boundary_conds
