SUBROUTINE particle_boundary_conds !------------------------------------------------------------------------------! ! Actual revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: particle_boundary_conds.f90 58 2007-03-09 14:27:38Z raasch $ ! 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