Index: palm/trunk/SOURCE/advec_particles.f90
===================================================================
--- palm/trunk/SOURCE/advec_particles.f90	(revision 58)
+++ palm/trunk/SOURCE/advec_particles.f90	(revision 59)
@@ -4,7 +4,8 @@
 ! Actual revisions:
 ! -----------------
-! Particle reflection at vertical walls (by Jin Zhang), regarding vertical
-! walls in the SGS model,
-! + user interface user_advec_particles
+! Particle reflection at vertical walls implemented in new subroutine
+! particle_boundary_conds,
+! vertical walls are regarded in the SGS model,
+! + user_advec_particles
 ! TEST: PRINT statements on unit 9 (commented out)
 !
@@ -1834,684 +1835,7 @@
        ENDDO   ! advection loop
 
-
 !
 !--    Particle reflection from walls
-       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
-      
-      particles_x = particles(n)%x
-      particles_y = particles(n)%y
-      particles_z = particles(n)%z
-
- 
-    IF (k2<=nzb_s_inner(j2,i2).and.nzb_s_inner(j2,i2)/=0) THEN
-
-   old_positions_x = particles(n)%x - particles(n)%speed_x * dt_particle
-   old_positions_y = particles(n)%y - particles(n)%speed_y * dt_particle
-   old_positions_z = particles(n)%z - particles(n)%speed_z * dt_particle
-   i1= (old_positions_x +0.5*dx) * ddx
-   j1 =(old_positions_y +0.5*dy) * ddy
-   k1 = old_positions_z / dz
-    
-
-!-- CASE 1
-
-   IF ( particles(n)%x  >  old_positions_x .and. &
-         particles(n)%y  >  old_positions_y )THEN
- 
-     
-
-
-     t_index = 1
-   
-     Do i = i1, i2
-     xline = i*dx+0.5*dx
-     
-     t(t_index) =(xline-old_positions_x)/(particles(n)%x-old_positions_x)
-     t_index = t_index + 1
-     ENDDO
-
-     Do j = j1, j2
-     yline = j*dy+0.5*dy
-     
-     t(t_index) =(yline-old_positions_y)/(particles(n)%y-old_positions_y)
-     t_index = t_index + 1
-     ENDDO
-
-     IF ( particles(n)%z < old_positions_z ) THEN
-     Do k = k1, k2 , -1
-     zline = k*dz
-     t(t_index) = (old_positions_z-zline)/(old_positions_z-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 )  THEN
-                           EXIT
-                            ENDIF
-                         ENDDO
-                        t(jr) = tmp_t
-                       
-             ENDDO
-        ENDDO
-    
-
-   
-   
-    Do t_index = 1, t_index_number
-  
-    positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x)
-    positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y)
-    positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z)
-
-     
-    
-     i3 = (positions_x + 0.5*dx) * ddx   
-     j3 = (positions_y + 0.5*dy) * ddy
-     k3 =  positions_z / dz
-
-     i5 = positions_x * ddx
-     j5 = positions_y * ddy
-     k5 = positions_z / dz
-
-     
- 
-
-
-    IF (k5 <= nzb_s_inner(j5,i3).and.nzb_s_inner(j5,i3)/=0 ) THEN
-             
-         IF ( positions_z == nzb_s_inner(j5,i3)*dz.and. &
-          k3 == nzb_s_inner(j5,i3) ) THEN
-      
-          particles(n)%z = 2*positions_z - particles_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 ( positions_z == nzb_s_inner(j3,i5)*dz.and. &
-          k3 == nzb_s_inner(j3,i5) ) THEN
-      
-          particles(n)%z = 2*positions_z - particles_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 ( positions_z == nzb_s_inner(j3,i3)*dz.and. &
-          k3 == nzb_s_inner(j3,i3) ) THEN
-      
-          particles(n)%z = 2*positions_z - particles_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 ( positions_y == j3*dy-0.5*dy .and. &
-         positions_z < nzb_s_inner(j3,i3)*dz) THEN
-              
-          particles(n)%y = 2*positions_y - particles_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 ( positions_x == i3*dx-0.5*dx.and. &
-            positions_z < nzb_s_inner(j3,i3)*dz ) THEN
-              
-          particles(n)%x = 2*positions_x - particles_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
-          GOTO 999
-        ENDIF  
-        
-        ENDIF
-     
-     ENDIF 
-    ENDDO
-   ENDIF
-    
-!-- CASE 2
-  
-     IF ( particles(n)%x  >  old_positions_x .and. &
-         particles(n)%y  <  old_positions_y )THEN
- 
-
-     t_index = 1
-   
-     Do i = i1, i2
-     xline = i*dx+0.5*dx
-     
-     t(t_index) =(xline-old_positions_x)/(particles(n)%x-old_positions_x)
-     t_index = t_index + 1
-     ENDDO
-
-     Do j = j1, j2,-1
-     yline = j*dy-0.5*dy
-     
-     t(t_index) =(old_positions_y-yline)/(old_positions_y-particles(n)%y)
-     t_index = t_index + 1
-     ENDDO
-
-     IF ( particles(n)%z < old_positions_z ) THEN
-     Do k = k1, k2 , -1
-     zline = k*dz
-     t(t_index) = (old_positions_z-zline)/(old_positions_z-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 )  THEN
-                           EXIT
-                            ENDIF
-                         ENDDO
-                        t(jr) = tmp_t
-                       
-             ENDDO
-        ENDDO
-
-     
-  
-     
-    
-   
-    Do t_index = 1, t_index_number
-  
-    positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x)
-    positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y)
-    positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z)
-
-     
-    
-     i3 = (positions_x + 0.5*dx) * ddx   
-     j3 = (positions_y + 0.5*dy) * ddy
-     k3 =  positions_z / dz
-
-     i5 = positions_x * ddx
-     j5 = positions_y * ddy
-     k5 = positions_z / dz
-    
-
- 
-    IF (k5 <= nzb_s_inner(j3,i5).and.nzb_s_inner(j3,i5)/=0 ) THEN
-             
-         IF ( positions_z == nzb_s_inner(j3,i5)*dz.and. &
-          k3 == nzb_s_inner(j3,i5) ) THEN
-      
-          particles(n)%z = 2*positions_z - particles_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 ( positions_z == nzb_s_inner(j3,i3)*dz.and. &
-          k3 == nzb_s_inner(j3,i3) ) THEN
-      
-          particles(n)%z = 2*positions_z - particles_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 ( positions_x == i3*dx-0.5*dx.and.&
-       positions_z < nzb_s_inner(j3,i3)*dz ) THEN
-              
-          particles(n)%x = 2*positions_x - particles_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 ( positions_z == nzb_s_inner(j5,i3)*dz.and. &
-          k3 == nzb_s_inner(j5,i3) ) THEN
-      
-          particles(n)%z = 2*positions_z - particles_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 ( positions_y == j5 * dy + 0.5*dy.and.&
-        positions_z < nzb_s_inner(j5,i3)*dz ) THEN
-      
-          particles(n)%y = 2*positions_y - particles_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 < old_positions_x .and. &
-         particles(n)%y > old_positions_y )THEN
- 
-
-     t_index = 1
-   
-     Do i = i1, i2,-1
-     xline = i*dx-0.5*dx
-     
-     t(t_index) =(old_positions_x-xline)/(old_positions_x-particles(n)%x)
-     t_index = t_index + 1
-     ENDDO
-
-     Do j = j1, j2
-     yline = j*dy+0.5*dy
-     
-     t(t_index) =(yline-old_positions_y)/(particles(n)%y-old_positions_y)
-     t_index = t_index + 1
-     ENDDO
-     
-     IF ( particles(n)%z < old_positions_z ) THEN
-     Do k = k1, k2 , -1
-     zline = k*dz
-     t(t_index) = (old_positions_z-zline)/(old_positions_z-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 )  THEN
-                           EXIT
-                            ENDIF
-                         ENDDO
-                        t(jr) = tmp_t
-                       
-             ENDDO
-        ENDDO
-
-    
-   
-    Do t_index = 1, t_index_number
-  
-    positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x)
-    positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y)
-    positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z)
-
-     
-    
-     i3 = (positions_x + 0.5*dx) * ddx   
-     j3 = (positions_y + 0.5*dy) * ddy
-     k3 =  positions_z / dz
-
-     i5 = positions_x * ddx
-     j5 = positions_y * ddy
-     k5 = positions_z / dz
-
-
-
-     IF (k5 <= nzb_s_inner(j5,i3).and.nzb_s_inner(j5,i3)/=0 ) THEN
-             
-         IF ( positions_z == nzb_s_inner(j5,i3)*dz.and. &
-          k3 == nzb_s_inner(j5,i3) ) THEN
-      
-          particles(n)%z = 2*positions_z - particles_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 ( positions_z == nzb_s_inner(j3,i3)*dz.and. &
-          k3 == nzb_s_inner(j3,i3) ) THEN
-      
-          particles(n)%z = 2*positions_z - particles_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 ( positions_y == j3*dy-0.5*dy .and. &
-            positions_z < nzb_s_inner(j3,i3)*dz) THEN
-              
-          particles(n)%y = 2*positions_y - particles_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 ( positions_z == nzb_s_inner(j3,i5)*dz.and. &
-          k3 == nzb_s_inner(j3,i5) ) THEN
-      
-          particles(n)%z = 2*positions_z - particles_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 ( positions_x == i5 * dx + 0.5*dx .and. &
-               positions_z < nzb_s_inner(j3,i5)*dz ) THEN
-      
-          particles(n)%x= 2*positions_x - particles_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 < old_positions_x .and. &
-         particles(n)%y  < old_positions_y )THEN
- 
-
-     t_index = 1
-   
-     Do i = i1, i2,-1
-     xline = i*dx-0.5*dx
-     
-     t(t_index) =(old_positions_x-xline)/(old_positions_x-particles(n)%x)
-     t_index = t_index + 1
-     
-     
-     
-     ENDDO
-
-     Do j = j1, j2,-1
-     yline = j*dy-0.5*dy
-     
-     t(t_index) =(old_positions_y-yline)/(old_positions_y-particles(n)%y)
-     t_index = t_index + 1
-     
-     
-     ENDDO
-
-     IF ( particles(n)%z < old_positions_z ) THEN
-     Do k = k1, k2 , -1
-     zline = k*dz
-     t(t_index) = (old_positions_z-zline)/(old_positions_z-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 )  THEN
-                           EXIT
-                            ENDIF
-                         ENDDO
-                        t(jr) = tmp_t
-                       
-             ENDDO
-        ENDDO
-    
-    
-    Do t_index = 1, t_index_number
-  
-    positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x)
-    positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y)
-    positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z)
-
-     
-    
-     i3 = (positions_x + 0.5*dx) * ddx   
-     j3 = (positions_y + 0.5*dy) * ddy
-     k3 =  positions_z / dz
-
-     i5 = positions_x * ddx
-     j5 = positions_y * ddy
-     k5 = positions_z / dz
-
-
-
-
-     IF (k3 <= nzb_s_inner(j3,i3).and.nzb_s_inner(j3,i3)/=0 ) THEN
-
-       IF ( positions_z == nzb_s_inner(j3,i3)*dz.and. &
-          k3 == nzb_s_inner(j3,i3) ) THEN
-      
-          particles(n)%z = 2*positions_z - particles_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 ( positions_z == nzb_s_inner(j3,i5)*dz.and. &
-          k3 == nzb_s_inner(j3,i5) ) THEN
-      
-          particles(n)%z = 2*positions_z - particles_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 ( positions_x == i5 * dx + 0.5*dx .and. nzb_s_inner(j3,i5)/=0.and.&
-                 positions_z < nzb_s_inner(j3,i5)*dz) THEN
-      
-          particles(n)%x= 2*positions_x - particles_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 ( positions_z == nzb_s_inner(j5,i3)*dz.and. &
-          k5 == nzb_s_inner(j5,i3) ) THEN
-      
-          particles(n)%z = 2*positions_z - particles_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 ( positions_y == j5 * dy + 0.5*dy .and. nzb_s_inner(j5,i3)/=0.and.&
-                    positions_z < nzb_s_inner(j5,i3)*dz ) THEN
-      
-          particles(n)%y= 2*positions_y - particles_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   ! particle reflection from walls
-
-       CALL cpu_log( log_point_s(48), 'advec_part_refle', 'stop' )
+       CALL particle_boundary_conds
 
 !
