SUBROUTINE particle_boundary_conds !------------------------------------------------------------------------------! ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: particle_boundary_conds.f90 484 2010-02-05 07:36:54Z maronga $ ! ! 150 2008-02-29 08:19:58Z raasch ! Vertical index calculations adjusted for ocean runs. ! ! Initial version (2007/03/09) ! ! Description: ! ------------ ! Calculates the reflection of particles from vertical walls. ! Routine developed by Jin Zhang (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 ! should be possible to further simplify/shorten it. ! THIS ROUTINE HAS NOT BEEN TESTED FOR OCEAN RUNS SO FAR! (see offset_ocean_*) !------------------------------------------------------------------------------! USE control_parameters USE cpulog USE grid_variables USE indices USE interfaces USE particle_attributes USE pegrid IMPLICIT NONE 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 LOGICAL :: reflect_x, reflect_y, reflect_z REAL :: dt_particle, pos_x, pos_x_old, pos_y, pos_y_old, pos_z, & pos_z_old, prt_x, prt_y, prt_z, tmp_t, xline, yline, zline 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 ) & 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)%speed_z_sgs = - particles(n)%speed_z_sgs 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)%speed_y_sgs = - particles(n)%speed_y_sgs 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)%speed_x_sgs = - particles(n)%speed_x_sgs ENDIF reflect_x = .FALSE. ENDIF ENDDO CALL cpu_log( log_point_s(48), 'advec_part_refle', 'stop' ) END SUBROUTINE particle_boundary_conds