SUBROUTINE lpm_boundary_conds( range ) !--------------------------------------------------------------------------------! ! This file is part of PALM. ! ! PALM is free software: you can redistribute it and/or modify it under the terms ! of the GNU General Public License as published by the Free Software Foundation, ! either version 3 of the License, or (at your option) any later version. ! ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License along with ! PALM. If not, see . ! ! Copyright 1997-2014 Leibniz Universitaet Hannover !--------------------------------------------------------------------------------! ! ! Current revisions: ! ----------------- ! ! Former revisions: ! ----------------- ! $Id: lpm_boundary_conds.f90 1321 2014-03-20 09:40:40Z witha $ ! ! 1320 2014-03-20 08:40:49Z raasch ! ONLY-attribute added to USE-statements, ! kind-parameters added to all INTEGER and REAL declaration statements, ! kinds are defined in new module kinds, ! revision history before 2012 removed, ! comment fields (!:) to be used for variable explanations added to ! all variable declaration statements ! ! 1036 2012-10-22 13:43:42Z raasch ! code put under GPL (PALM 3.9) ! ! 849 2012-03-15 10:35:09Z raasch ! routine renamed lpm_boundary_conds, bottom and top boundary conditions ! included (former part of advec_particles) ! ! 824 2012-02-17 09:09:57Z raasch ! particle attributes speed_x|y|z_sgs renamed rvar1|2|3 ! ! Initial version (2007/03/09) ! ! Description: ! ------------ ! 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 ! should be possible to further simplify/shorten it. ! ! THE WALLS PART OF THIS ROUTINE HAS NOT BEEN TESTED FOR OCEAN RUNS SO FAR!!!! ! (see offset_ocean_*) !------------------------------------------------------------------------------! USE arrays_3d, & ONLY: zu, zw USE control_parameters, & ONLY: dz, message_string, particle_maximum_age USE cpulog, & ONLY: cpu_log, log_point_s USE grid_variables, & ONLY: ddx, dx, ddy, dy USE indices, & ONLY: nxl, nxr, nyn, nys, nz, nzb_s_inner USE kinds USE particle_attributes, & ONLY: deleted_particles, deleted_tails, ibc_par_b, ibc_par_t, & number_of_particles, particles, particle_mask, & particle_tail_coordinates, particle_type, offset_ocean_nzt_m1, & tail_mask, use_particle_tails, use_sgs_for_particles USE pegrid IMPLICIT NONE CHARACTER (LEN=*) :: range !: INTEGER(iwp) :: i !: INTEGER(iwp) :: inc !: INTEGER(iwp) :: ir !: INTEGER(iwp) :: i1 !: INTEGER(iwp) :: i2 !: INTEGER(iwp) :: i3 !: INTEGER(iwp) :: i5 !: INTEGER(iwp) :: j !: INTEGER(iwp) :: jr !: INTEGER(iwp) :: j1 !: INTEGER(iwp) :: j2 !: INTEGER(iwp) :: j3 !: INTEGER(iwp) :: j5 !: INTEGER(iwp) :: k !: INTEGER(iwp) :: k1 !: INTEGER(iwp) :: k2 !: INTEGER(iwp) :: k3 !: INTEGER(iwp) :: k5 !: INTEGER(iwp) :: n !: INTEGER(iwp) :: nn !: INTEGER(iwp) :: t_index !: INTEGER(iwp) :: t_index_number !: LOGICAL :: reflect_x !: LOGICAL :: reflect_y !: LOGICAL :: reflect_z !: REAL(wp) :: dt_particle !: REAL(wp) :: pos_x !: REAL(wp) :: pos_x_old !: REAL(wp) :: pos_y !: REAL(wp) :: pos_y_old !: REAL(wp) :: pos_z !: REAL(wp) :: pos_z_old !: REAL(wp) :: prt_x !: REAL(wp) :: prt_y !: REAL(wp) :: prt_z !: REAL(wp) :: t(1:200) !: REAL(wp) :: tmp_t !: REAL(wp) :: xline !: REAL(wp) :: yline !: REAL(wp) :: zline !: 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 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 CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'stop' ) ENDIF END SUBROUTINE lpm_boundary_conds