Changeset 849 for palm/trunk/SOURCE/lpm_boundary_conds.f90
- Timestamp:
- Mar 15, 2012 10:35:09 AM (12 years ago)
- File:
-
- 1 moved
-
palm/trunk/SOURCE/lpm_boundary_conds.f90 (moved) (moved from palm/trunk/SOURCE/particle_boundary_conds.f90) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lpm_boundary_conds.f90
r848 r849 1 SUBROUTINE particle_boundary_conds1 SUBROUTINE lpm_boundary_conds( range ) 2 2 3 3 !------------------------------------------------------------------------------! 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! routine renamed lpm_boundary_conds, bottom and top boundary conditions 7 ! included 7 8 ! 8 9 ! Former revisions: … … 20 21 ! Description: 21 22 ! ------------ 22 ! Calculates the reflection of particles from vertical walls. 23 ! Routine developed by Jin Zhang (2006-2007) 23 ! Boundary conditions for the Lagrangian particles. 24 ! The routine consists of two different parts. One handles the bottom (flat) 25 ! and top boundary. In this part, also particles which exceeded their lifetime 26 ! are deleted. 27 ! The other part handles the reflection of particles from vertical walls. 28 ! This part was developed by Jin Zhang during 2006-2007. 29 ! 24 30 ! To do: Code structure for finding the t_index values and for checking the 25 ! 31 ! ----- reflection conditions is basically the same for all four cases, so it 26 32 ! should be possible to further simplify/shorten it. 27 ! THIS ROUTINE HAS NOT BEEN TESTED FOR OCEAN RUNS SO FAR! (see offset_ocean_*) 33 ! 34 ! THE WALLS PART OF THIS ROUTINE HAS NOT BEEN TESTED FOR OCEAN RUNS SO FAR!!!! 35 ! (see offset_ocean_*) 28 36 !------------------------------------------------------------------------------! 29 37 38 USE arrays_3d 30 39 USE control_parameters 31 40 USE cpulog … … 38 47 IMPLICIT NONE 39 48 49 CHARACTER (LEN=*) :: range 50 40 51 INTEGER :: i, inc, ir, i1, i2, i3, i5, j, jr, j1, j2, j3, j5, k, k1, k2, & 41 k3, k5, n, t_index, t_index_number52 k3, k5, n, nn, t_index, t_index_number 42 53 43 54 LOGICAL :: reflect_x, reflect_y, reflect_z … … 48 59 REAL :: t(1:200) 49 60 50 CALL cpu_log( log_point_s(48), 'advec_part_refle', 'start' ) 51 52 53 54 reflect_x = .FALSE. 55 reflect_y = .FALSE. 56 reflect_z = .FALSE. 57 58 DO n = 1, number_of_particles 59 60 dt_particle = particles(n)%age - particles(n)%age_m 61 62 i2 = ( particles(n)%x + 0.5 * dx ) * ddx 63 j2 = ( particles(n)%y + 0.5 * dy ) * ddy 64 k2 = particles(n)%z / dz + 1 + offset_ocean_nzt_m1 65 66 prt_x = particles(n)%x 67 prt_y = particles(n)%y 68 prt_z = particles(n)%z 69 70 ! 71 !-- If the particle position is below the surface, it has to be reflected 72 IF ( k2 <= nzb_s_inner(j2,i2) .AND. nzb_s_inner(j2,i2) /=0 ) THEN 73 74 pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle 75 pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle 76 pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle 77 i1 = ( pos_x_old + 0.5 * dx ) * ddx 78 j1 = ( pos_y_old + 0.5 * dy ) * ddy 79 k1 = pos_z_old / dz + offset_ocean_nzt_m1 80 81 ! 82 !-- Case 1 83 IF ( particles(n)%x > pos_x_old .AND. particles(n)%y > pos_y_old ) & 61 62 63 IF ( range == 'bottom/top' ) THEN 64 65 ! 66 !-- Apply boundary conditions to those particles that have crossed the top or 67 !-- bottom boundary and delete those particles, which are older than allowed 68 DO n = 1, number_of_particles 69 70 nn = particles(n)%tail_id 71 72 ! 73 !-- Stop if particles have moved further than the length of one 74 !-- PE subdomain (newly released particles have age = age_m!) 75 IF ( particles(n)%age /= particles(n)%age_m ) THEN 76 IF ( ABS(particles(n)%speed_x) > & 77 ((nxr-nxl+2)*dx)/(particles(n)%age-particles(n)%age_m) .OR. & 78 ABS(particles(n)%speed_y) > & 79 ((nyn-nys+2)*dy)/(particles(n)%age-particles(n)%age_m) ) THEN 80 81 WRITE( message_string, * ) 'particle too fast. n = ', n 82 CALL message( 'lpm_boundary_conds', 'PA0148', 2, 2, -1, 6, 1 ) 83 ENDIF 84 ENDIF 85 86 IF ( particles(n)%age > particle_maximum_age .AND. & 87 particle_mask(n) ) & 84 88 THEN 85 t_index = 1 86 87 DO i = i1, i2 88 xline = i * dx + 0.5 * dx 89 t(t_index) = ( xline - pos_x_old ) / & 90 ( particles(n)%x - pos_x_old ) 91 t_index = t_index + 1 92 ENDDO 93 94 DO j = j1, j2 95 yline = j * dy + 0.5 * dy 96 t(t_index) = ( yline - pos_y_old ) / & 97 ( particles(n)%y - pos_y_old ) 98 t_index = t_index + 1 99 ENDDO 100 101 IF ( particles(n)%z < pos_z_old ) THEN 102 DO k = k1, k2 , -1 103 zline = k * dz 104 t(t_index) = ( pos_z_old - zline ) / & 105 ( pos_z_old - particles(n)%z ) 106 t_index = t_index + 1 107 ENDDO 108 ENDIF 109 110 t_index_number = t_index - 1 111 112 ! 113 !-- Sorting t 114 inc = 1 115 jr = 1 116 DO WHILE ( inc <= t_index_number ) 117 inc = 3 * inc + 1 118 ENDDO 119 120 DO WHILE ( inc > 1 ) 121 inc = inc / 3 122 DO ir = inc+1, t_index_number 123 tmp_t = t(ir) 124 jr = ir 125 DO WHILE ( t(jr-inc) > tmp_t ) 126 t(jr) = t(jr-inc) 127 jr = jr - inc 128 IF ( jr <= inc ) EXIT 129 ENDDO 130 t(jr) = tmp_t 131 ENDDO 132 ENDDO 133 134 case1: DO t_index = 1, t_index_number 135 136 pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old ) 137 pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old ) 138 pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old ) 139 140 i3 = ( pos_x + 0.5 * dx ) * ddx 141 j3 = ( pos_y + 0.5 * dy ) * ddy 142 k3 = pos_z / dz + offset_ocean_nzt_m1 143 144 i5 = pos_x * ddx 145 j5 = pos_y * ddy 146 k5 = pos_z / dz + offset_ocean_nzt_m1 147 148 IF ( k5 <= nzb_s_inner(j5,i3) .AND. & 149 nzb_s_inner(j5,i3) /= 0 ) THEN 150 151 IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. & 152 k3 == nzb_s_inner(j5,i3) ) THEN 153 reflect_z = .TRUE. 154 EXIT case1 155 ENDIF 156 157 ENDIF 158 159 IF ( k5 <= nzb_s_inner(j3,i5) .AND. & 160 nzb_s_inner(j3,i5) /= 0 ) THEN 161 162 IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. & 163 k3 == nzb_s_inner(j3,i5) ) THEN 164 reflect_z = .TRUE. 165 EXIT case1 166 ENDIF 167 168 ENDIF 169 170 IF ( k3 <= nzb_s_inner(j3,i3) .AND. & 171 nzb_s_inner(j3,i3) /= 0 ) THEN 172 173 IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. & 174 k3 == nzb_s_inner(j3,i3) ) THEN 175 reflect_z = .TRUE. 176 EXIT case1 177 ENDIF 178 179 IF ( pos_y == ( j3 * dy - 0.5 * dy ) .AND. & 180 pos_z < nzb_s_inner(j3,i3) * dz ) THEN 181 reflect_y = .TRUE. 182 EXIT case1 183 ENDIF 184 185 IF ( pos_x == ( i3 * dx - 0.5 * dx ) .AND. & 186 pos_z < nzb_s_inner(j3,i3) * dz ) THEN 187 reflect_x = .TRUE. 188 EXIT case1 189 ENDIF 190 191 ENDIF 192 193 ENDDO case1 194 195 ! 196 !-- Case 2 197 ELSEIF ( particles(n)%x > pos_x_old .AND. & 198 particles(n)%y < pos_y_old ) THEN 199 200 t_index = 1 201 202 DO i = i1, i2 203 xline = i * dx + 0.5 * dx 204 t(t_index) = ( xline - pos_x_old ) / & 205 ( particles(n)%x - pos_x_old ) 206 t_index = t_index + 1 207 ENDDO 208 209 DO j = j1, j2, -1 210 yline = j * dy - 0.5 * dy 211 t(t_index) = ( pos_y_old - yline ) / & 212 ( pos_y_old - particles(n)%y ) 213 t_index = t_index + 1 214 ENDDO 215 216 IF ( particles(n)%z < pos_z_old ) THEN 217 DO k = k1, k2 , -1 218 zline = k * dz 219 t(t_index) = ( pos_z_old - zline ) / & 220 ( pos_z_old - particles(n)%z ) 221 t_index = t_index + 1 222 ENDDO 223 ENDIF 224 t_index_number = t_index-1 225 226 ! 227 !-- Sorting t 228 inc = 1 229 jr = 1 230 DO WHILE ( inc <= t_index_number ) 231 inc = 3 * inc + 1 232 ENDDO 233 234 DO WHILE ( inc > 1 ) 235 inc = inc / 3 236 DO ir = inc+1, t_index_number 237 tmp_t = t(ir) 238 jr = ir 239 DO WHILE ( t(jr-inc) > tmp_t ) 240 t(jr) = t(jr-inc) 241 jr = jr - inc 242 IF ( jr <= inc ) EXIT 243 ENDDO 244 t(jr) = tmp_t 245 ENDDO 246 ENDDO 247 248 case2: DO t_index = 1, t_index_number 249 250 pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old ) 251 pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old ) 252 pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old ) 253 254 i3 = ( pos_x + 0.5 * dx ) * ddx 255 j3 = ( pos_y + 0.5 * dy ) * ddy 256 k3 = pos_z / dz + offset_ocean_nzt_m1 257 258 i5 = pos_x * ddx 259 j5 = pos_y * ddy 260 k5 = pos_z / dz + offset_ocean_nzt_m1 261 262 IF ( k5 <= nzb_s_inner(j3,i5) .AND. & 263 nzb_s_inner(j3,i5) /= 0 ) THEN 264 265 IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. & 266 k3 == nzb_s_inner(j3,i5) ) THEN 267 reflect_z = .TRUE. 268 EXIT case2 269 ENDIF 270 271 ENDIF 272 273 IF ( k3 <= nzb_s_inner(j3,i3) .AND. & 274 nzb_s_inner(j3,i3) /= 0 ) THEN 275 276 IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. & 277 k3 == nzb_s_inner(j3,i3) ) THEN 278 reflect_z = .TRUE. 279 EXIT case2 280 ENDIF 281 282 IF ( pos_x == ( i3 * dx - 0.5 * dx ) .AND. & 283 pos_z < nzb_s_inner(j3,i3) * dz ) THEN 284 reflect_x = .TRUE. 285 EXIT case2 286 ENDIF 287 288 ENDIF 289 290 IF ( k5 <= nzb_s_inner(j5,i3) .AND. & 291 nzb_s_inner(j5,i3) /= 0 ) THEN 292 293 IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. & 294 k3 == nzb_s_inner(j5,i3) ) THEN 295 reflect_z = .TRUE. 296 EXIT case2 297 ENDIF 298 299 IF ( pos_y == ( j5 * dy + 0.5 * dy ) .AND. & 300 pos_z < nzb_s_inner(j5,i3) * dz ) THEN 301 reflect_y = .TRUE. 302 EXIT case2 303 ENDIF 304 305 ENDIF 306 307 ENDDO case2 308 309 ! 310 !-- Case 3 311 ELSEIF ( particles(n)%x < pos_x_old .AND. & 312 particles(n)%y > pos_y_old ) THEN 313 314 t_index = 1 315 316 DO i = i1, i2, -1 317 xline = i * dx - 0.5 * dx 318 t(t_index) = ( pos_x_old - xline ) / & 319 ( pos_x_old - particles(n)%x ) 320 t_index = t_index + 1 321 ENDDO 322 323 DO j = j1, j2 324 yline = j * dy + 0.5 * dy 325 t(t_index) = ( yline - pos_y_old ) / & 326 ( particles(n)%y - pos_y_old ) 327 t_index = t_index + 1 328 ENDDO 329 330 IF ( particles(n)%z < pos_z_old ) THEN 331 DO k = k1, k2 , -1 332 zline = k * dz 333 t(t_index) = ( pos_z_old - zline ) / & 334 ( pos_z_old - particles(n)%z ) 335 t_index = t_index + 1 336 ENDDO 337 ENDIF 338 t_index_number = t_index - 1 339 340 ! 341 !-- Sorting t 342 inc = 1 343 jr = 1 344 345 DO WHILE ( inc <= t_index_number ) 346 inc = 3 * inc + 1 347 ENDDO 348 349 DO WHILE ( inc > 1 ) 350 inc = inc / 3 351 DO ir = inc+1, t_index_number 352 tmp_t = t(ir) 353 jr = ir 354 DO WHILE ( t(jr-inc) > tmp_t ) 355 t(jr) = t(jr-inc) 356 jr = jr - inc 357 IF ( jr <= inc ) EXIT 358 ENDDO 359 t(jr) = tmp_t 360 ENDDO 361 ENDDO 362 363 case3: DO t_index = 1, t_index_number 364 365 pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old ) 366 pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old ) 367 pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old ) 368 369 i3 = ( pos_x + 0.5 * dx ) * ddx 370 j3 = ( pos_y + 0.5 * dy ) * ddy 371 k3 = pos_z / dz + offset_ocean_nzt_m1 372 373 i5 = pos_x * ddx 374 j5 = pos_y * ddy 375 k5 = pos_z / dz + offset_ocean_nzt_m1 376 377 IF ( k5 <= nzb_s_inner(j5,i3) .AND. & 378 nzb_s_inner(j5,i3) /= 0 ) THEN 379 380 IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. & 381 k3 == nzb_s_inner(j5,i3) ) THEN 382 reflect_z = .TRUE. 383 EXIT case3 384 ENDIF 385 386 ENDIF 387 388 IF ( k3 <= nzb_s_inner(j3,i3) .AND. & 389 nzb_s_inner(j3,i3) /= 0 ) THEN 390 391 IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. & 392 k3 == nzb_s_inner(j3,i3) ) THEN 393 reflect_z = .TRUE. 394 EXIT case3 395 ENDIF 396 397 IF ( pos_y == ( j3 * dy - 0.5 * dy ) .AND. & 398 pos_z < nzb_s_inner(j3,i3) * dz ) THEN 399 reflect_y = .TRUE. 400 EXIT case3 401 ENDIF 402 403 ENDIF 404 405 IF ( k5 <= nzb_s_inner(j3,i5) .AND. & 406 nzb_s_inner(j3,i5) /= 0 ) THEN 407 408 IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. & 409 k3 == nzb_s_inner(j3,i5) ) THEN 410 reflect_z = .TRUE. 411 EXIT case3 412 ENDIF 413 414 IF ( pos_x == ( i5 * dx + 0.5 * dx ) .AND. & 415 pos_z < nzb_s_inner(j3,i5) * dz ) THEN 416 reflect_x = .TRUE. 417 EXIT case3 418 ENDIF 419 420 ENDIF 421 422 ENDDO case3 423 424 ! 425 !-- Case 4 426 ELSEIF ( particles(n)%x < pos_x_old .AND. & 427 particles(n)%y < pos_y_old ) THEN 428 429 t_index = 1 430 431 DO i = i1, i2, -1 432 xline = i * dx - 0.5 * dx 433 t(t_index) = ( pos_x_old - xline ) / & 434 ( pos_x_old - particles(n)%x ) 435 t_index = t_index + 1 436 ENDDO 437 438 DO j = j1, j2, -1 439 yline = j * dy - 0.5 * dy 440 t(t_index) = ( pos_y_old - yline ) / & 441 ( pos_y_old - particles(n)%y ) 442 t_index = t_index + 1 443 ENDDO 444 445 IF ( particles(n)%z < pos_z_old ) THEN 446 DO k = k1, k2 , -1 447 zline = k * dz 448 t(t_index) = ( pos_z_old - zline ) / & 449 ( pos_z_old-particles(n)%z ) 450 t_index = t_index + 1 451 ENDDO 452 ENDIF 453 t_index_number = t_index-1 454 455 ! 456 !-- Sorting t 457 inc = 1 458 jr = 1 459 460 DO WHILE ( inc <= t_index_number ) 461 inc = 3 * inc + 1 462 ENDDO 463 464 DO WHILE ( inc > 1 ) 465 inc = inc / 3 466 DO ir = inc+1, t_index_number 467 tmp_t = t(ir) 468 jr = ir 469 DO WHILE ( t(jr-inc) > tmp_t ) 470 t(jr) = t(jr-inc) 471 jr = jr - inc 472 IF ( jr <= inc ) EXIT 473 ENDDO 474 t(jr) = tmp_t 475 ENDDO 476 ENDDO 477 478 case4: DO t_index = 1, t_index_number 479 480 pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old ) 481 pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old ) 482 pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old ) 483 484 i3 = ( pos_x + 0.5 * dx ) * ddx 485 j3 = ( pos_y + 0.5 * dy ) * ddy 486 k3 = pos_z / dz + offset_ocean_nzt_m1 487 488 i5 = pos_x * ddx 489 j5 = pos_y * ddy 490 k5 = pos_z / dz + offset_ocean_nzt_m1 491 492 IF ( k3 <= nzb_s_inner(j3,i3) .AND. & 493 nzb_s_inner(j3,i3) /= 0 ) THEN 494 495 IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. & 496 k3 == nzb_s_inner(j3,i3) ) THEN 497 reflect_z = .TRUE. 498 EXIT case4 499 ENDIF 500 501 ENDIF 502 503 IF ( k5 <= nzb_s_inner(j3,i5) .AND. & 504 nzb_s_inner(j3,i5) /= 0 ) THEN 505 506 IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. & 507 k3 == nzb_s_inner(j3,i5) ) THEN 508 reflect_z = .TRUE. 509 EXIT case4 510 ENDIF 511 512 IF ( pos_x == ( i5 * dx + 0.5 * dx ) .AND. & 513 nzb_s_inner(j3,i5) /=0 .AND. & 514 pos_z < nzb_s_inner(j3,i5) * dz ) THEN 515 reflect_x = .TRUE. 516 EXIT case4 517 ENDIF 518 519 ENDIF 520 521 IF ( k5 <= nzb_s_inner(j5,i3) .AND. & 522 nzb_s_inner(j5,i3) /= 0 ) THEN 523 524 IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. & 525 k5 == nzb_s_inner(j5,i3) ) THEN 526 reflect_z = .TRUE. 527 EXIT case4 528 ENDIF 529 530 IF ( pos_y == ( j5 * dy + 0.5 * dy ) .AND. & 531 nzb_s_inner(j5,i3) /= 0 .AND. & 532 pos_z < nzb_s_inner(j5,i3) * dz ) THEN 533 reflect_y = .TRUE. 534 EXIT case4 535 ENDIF 536 537 ENDIF 89 particle_mask(n) = .FALSE. 90 deleted_particles = deleted_particles + 1 91 IF ( use_particle_tails .AND. nn /= 0 ) THEN 92 tail_mask(nn) = .FALSE. 93 deleted_tails = deleted_tails + 1 94 ENDIF 95 ENDIF 96 97 IF ( particles(n)%z >= zu(nz) .AND. particle_mask(n) ) THEN 98 IF ( ibc_par_t == 1 ) THEN 99 ! 100 !-- Particle absorption 101 particle_mask(n) = .FALSE. 102 deleted_particles = deleted_particles + 1 103 IF ( use_particle_tails .AND. nn /= 0 ) THEN 104 tail_mask(nn) = .FALSE. 105 deleted_tails = deleted_tails + 1 106 ENDIF 107 ELSEIF ( ibc_par_t == 2 ) THEN 108 ! 109 !-- Particle reflection 110 particles(n)%z = 2.0 * zu(nz) - particles(n)%z 111 particles(n)%speed_z = -particles(n)%speed_z 112 IF ( use_sgs_for_particles .AND. & 113 particles(n)%rvar3 > 0.0 ) THEN 114 particles(n)%rvar3 = -particles(n)%rvar3 115 ENDIF 116 IF ( use_particle_tails .AND. nn /= 0 ) THEN 117 particle_tail_coordinates(1,3,nn) = 2.0 * zu(nz) - & 118 particle_tail_coordinates(1,3,nn) 119 ENDIF 120 ENDIF 121 ENDIF 122 IF ( particles(n)%z < zw(0) .AND. particle_mask(n) ) THEN 123 IF ( ibc_par_b == 1 ) THEN 124 ! 125 !-- Particle absorption 126 particle_mask(n) = .FALSE. 127 deleted_particles = deleted_particles + 1 128 IF ( use_particle_tails .AND. nn /= 0 ) THEN 129 tail_mask(nn) = .FALSE. 130 deleted_tails = deleted_tails + 1 131 ENDIF 132 ELSEIF ( ibc_par_b == 2 ) THEN 133 ! 134 !-- Particle reflection 135 particles(n)%z = 2.0 * zw(0) - particles(n)%z 136 particles(n)%speed_z = -particles(n)%speed_z 137 IF ( use_sgs_for_particles .AND. & 138 particles(n)%rvar3 < 0.0 ) THEN 139 particles(n)%rvar3 = -particles(n)%rvar3 140 ENDIF 141 IF ( use_particle_tails .AND. nn /= 0 ) THEN 142 particle_tail_coordinates(1,3,nn) = 2.0 * zu(nz) - & 143 particle_tail_coordinates(1,3,nn) 144 ENDIF 145 IF ( use_particle_tails .AND. nn /= 0 ) THEN 146 particle_tail_coordinates(1,3,nn) = 2.0 * zw(0) - & 147 particle_tail_coordinates(1,3,nn) 148 ENDIF 149 ENDIF 150 ENDIF 151 ENDDO 152 153 ELSEIF ( range == 'walls' ) THEN 154 155 CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'start' ) 156 157 reflect_x = .FALSE. 158 reflect_y = .FALSE. 159 reflect_z = .FALSE. 160 161 DO n = 1, number_of_particles 162 163 dt_particle = particles(n)%age - particles(n)%age_m 164 165 i2 = ( particles(n)%x + 0.5 * dx ) * ddx 166 j2 = ( particles(n)%y + 0.5 * dy ) * ddy 167 k2 = particles(n)%z / dz + 1 + offset_ocean_nzt_m1 168 169 prt_x = particles(n)%x 170 prt_y = particles(n)%y 171 prt_z = particles(n)%z 172 173 ! 174 !-- If the particle position is below the surface, it has to be reflected 175 IF ( k2 <= nzb_s_inner(j2,i2) .AND. nzb_s_inner(j2,i2) /=0 ) THEN 176 177 pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle 178 pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle 179 pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle 180 i1 = ( pos_x_old + 0.5 * dx ) * ddx 181 j1 = ( pos_y_old + 0.5 * dy ) * ddy 182 k1 = pos_z_old / dz + offset_ocean_nzt_m1 183 184 ! 185 !-- Case 1 186 IF ( particles(n)%x > pos_x_old .AND. particles(n)%y > pos_y_old )& 187 THEN 188 t_index = 1 189 190 DO i = i1, i2 191 xline = i * dx + 0.5 * dx 192 t(t_index) = ( xline - pos_x_old ) / & 193 ( particles(n)%x - pos_x_old ) 194 t_index = t_index + 1 195 ENDDO 196 197 DO j = j1, j2 198 yline = j * dy + 0.5 * dy 199 t(t_index) = ( yline - pos_y_old ) / & 200 ( particles(n)%y - pos_y_old ) 201 t_index = t_index + 1 202 ENDDO 203 204 IF ( particles(n)%z < pos_z_old ) THEN 205 DO k = k1, k2 , -1 206 zline = k * dz 207 t(t_index) = ( pos_z_old - zline ) / & 208 ( pos_z_old - particles(n)%z ) 209 t_index = t_index + 1 210 ENDDO 211 ENDIF 212 213 t_index_number = t_index - 1 214 215 ! 216 !-- Sorting t 217 inc = 1 218 jr = 1 219 DO WHILE ( inc <= t_index_number ) 220 inc = 3 * inc + 1 221 ENDDO 222 223 DO WHILE ( inc > 1 ) 224 inc = inc / 3 225 DO ir = inc+1, t_index_number 226 tmp_t = t(ir) 227 jr = ir 228 DO WHILE ( t(jr-inc) > tmp_t ) 229 t(jr) = t(jr-inc) 230 jr = jr - inc 231 IF ( jr <= inc ) EXIT 232 ENDDO 233 t(jr) = tmp_t 234 ENDDO 235 ENDDO 236 237 case1: DO t_index = 1, t_index_number 238 239 pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old ) 240 pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old ) 241 pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old ) 242 243 i3 = ( pos_x + 0.5 * dx ) * ddx 244 j3 = ( pos_y + 0.5 * dy ) * ddy 245 k3 = pos_z / dz + offset_ocean_nzt_m1 246 247 i5 = pos_x * ddx 248 j5 = pos_y * ddy 249 k5 = pos_z / dz + offset_ocean_nzt_m1 250 251 IF ( k5 <= nzb_s_inner(j5,i3) .AND. & 252 nzb_s_inner(j5,i3) /= 0 ) THEN 253 254 IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. & 255 k3 == nzb_s_inner(j5,i3) ) THEN 256 reflect_z = .TRUE. 257 EXIT case1 258 ENDIF 259 260 ENDIF 261 262 IF ( k5 <= nzb_s_inner(j3,i5) .AND. & 263 nzb_s_inner(j3,i5) /= 0 ) THEN 264 265 IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. & 266 k3 == nzb_s_inner(j3,i5) ) THEN 267 reflect_z = .TRUE. 268 EXIT case1 269 ENDIF 270 271 ENDIF 272 273 IF ( k3 <= nzb_s_inner(j3,i3) .AND. & 274 nzb_s_inner(j3,i3) /= 0 ) THEN 275 276 IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. & 277 k3 == nzb_s_inner(j3,i3) ) THEN 278 reflect_z = .TRUE. 279 EXIT case1 280 ENDIF 281 282 IF ( pos_y == ( j3 * dy - 0.5 * dy ) .AND. & 283 pos_z < nzb_s_inner(j3,i3) * dz ) THEN 284 reflect_y = .TRUE. 285 EXIT case1 286 ENDIF 287 288 IF ( pos_x == ( i3 * dx - 0.5 * dx ) .AND. & 289 pos_z < nzb_s_inner(j3,i3) * dz ) THEN 290 reflect_x = .TRUE. 291 EXIT case1 292 ENDIF 293 294 ENDIF 295 296 ENDDO case1 297 298 ! 299 !-- Case 2 300 ELSEIF ( particles(n)%x > pos_x_old .AND. & 301 particles(n)%y < pos_y_old ) THEN 302 303 t_index = 1 304 305 DO i = i1, i2 306 xline = i * dx + 0.5 * dx 307 t(t_index) = ( xline - pos_x_old ) / & 308 ( particles(n)%x - pos_x_old ) 309 t_index = t_index + 1 310 ENDDO 311 312 DO j = j1, j2, -1 313 yline = j * dy - 0.5 * dy 314 t(t_index) = ( pos_y_old - yline ) / & 315 ( pos_y_old - particles(n)%y ) 316 t_index = t_index + 1 317 ENDDO 318 319 IF ( particles(n)%z < pos_z_old ) THEN 320 DO k = k1, k2 , -1 321 zline = k * dz 322 t(t_index) = ( pos_z_old - zline ) / & 323 ( pos_z_old - particles(n)%z ) 324 t_index = t_index + 1 325 ENDDO 326 ENDIF 327 t_index_number = t_index-1 328 329 ! 330 !-- Sorting t 331 inc = 1 332 jr = 1 333 DO WHILE ( inc <= t_index_number ) 334 inc = 3 * inc + 1 335 ENDDO 336 337 DO WHILE ( inc > 1 ) 338 inc = inc / 3 339 DO ir = inc+1, t_index_number 340 tmp_t = t(ir) 341 jr = ir 342 DO WHILE ( t(jr-inc) > tmp_t ) 343 t(jr) = t(jr-inc) 344 jr = jr - inc 345 IF ( jr <= inc ) EXIT 346 ENDDO 347 t(jr) = tmp_t 348 ENDDO 349 ENDDO 350 351 case2: DO t_index = 1, t_index_number 352 353 pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old ) 354 pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old ) 355 pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old ) 356 357 i3 = ( pos_x + 0.5 * dx ) * ddx 358 j3 = ( pos_y + 0.5 * dy ) * ddy 359 k3 = pos_z / dz + offset_ocean_nzt_m1 360 361 i5 = pos_x * ddx 362 j5 = pos_y * ddy 363 k5 = pos_z / dz + offset_ocean_nzt_m1 364 365 IF ( k5 <= nzb_s_inner(j3,i5) .AND. & 366 nzb_s_inner(j3,i5) /= 0 ) THEN 367 368 IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. & 369 k3 == nzb_s_inner(j3,i5) ) THEN 370 reflect_z = .TRUE. 371 EXIT case2 372 ENDIF 373 374 ENDIF 375 376 IF ( k3 <= nzb_s_inner(j3,i3) .AND. & 377 nzb_s_inner(j3,i3) /= 0 ) THEN 378 379 IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. & 380 k3 == nzb_s_inner(j3,i3) ) THEN 381 reflect_z = .TRUE. 382 EXIT case2 383 ENDIF 384 385 IF ( pos_x == ( i3 * dx - 0.5 * dx ) .AND. & 386 pos_z < nzb_s_inner(j3,i3) * dz ) THEN 387 reflect_x = .TRUE. 388 EXIT case2 389 ENDIF 390 391 ENDIF 392 393 IF ( k5 <= nzb_s_inner(j5,i3) .AND. & 394 nzb_s_inner(j5,i3) /= 0 ) THEN 395 396 IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. & 397 k3 == nzb_s_inner(j5,i3) ) THEN 398 reflect_z = .TRUE. 399 EXIT case2 400 ENDIF 401 402 IF ( pos_y == ( j5 * dy + 0.5 * dy ) .AND. & 403 pos_z < nzb_s_inner(j5,i3) * dz ) THEN 404 reflect_y = .TRUE. 405 EXIT case2 406 ENDIF 407 408 ENDIF 409 410 ENDDO case2 411 412 ! 413 !-- Case 3 414 ELSEIF ( particles(n)%x < pos_x_old .AND. & 415 particles(n)%y > pos_y_old ) THEN 416 417 t_index = 1 418 419 DO i = i1, i2, -1 420 xline = i * dx - 0.5 * dx 421 t(t_index) = ( pos_x_old - xline ) / & 422 ( pos_x_old - particles(n)%x ) 423 t_index = t_index + 1 424 ENDDO 425 426 DO j = j1, j2 427 yline = j * dy + 0.5 * dy 428 t(t_index) = ( yline - pos_y_old ) / & 429 ( particles(n)%y - pos_y_old ) 430 t_index = t_index + 1 431 ENDDO 432 433 IF ( particles(n)%z < pos_z_old ) THEN 434 DO k = k1, k2 , -1 435 zline = k * dz 436 t(t_index) = ( pos_z_old - zline ) / & 437 ( pos_z_old - particles(n)%z ) 438 t_index = t_index + 1 439 ENDDO 440 ENDIF 441 t_index_number = t_index - 1 442 443 ! 444 !-- Sorting t 445 inc = 1 446 jr = 1 447 448 DO WHILE ( inc <= t_index_number ) 449 inc = 3 * inc + 1 450 ENDDO 451 452 DO WHILE ( inc > 1 ) 453 inc = inc / 3 454 DO ir = inc+1, t_index_number 455 tmp_t = t(ir) 456 jr = ir 457 DO WHILE ( t(jr-inc) > tmp_t ) 458 t(jr) = t(jr-inc) 459 jr = jr - inc 460 IF ( jr <= inc ) EXIT 461 ENDDO 462 t(jr) = tmp_t 463 ENDDO 464 ENDDO 465 466 case3: DO t_index = 1, t_index_number 467 468 pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old ) 469 pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old ) 470 pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old ) 471 472 i3 = ( pos_x + 0.5 * dx ) * ddx 473 j3 = ( pos_y + 0.5 * dy ) * ddy 474 k3 = pos_z / dz + offset_ocean_nzt_m1 475 476 i5 = pos_x * ddx 477 j5 = pos_y * ddy 478 k5 = pos_z / dz + offset_ocean_nzt_m1 479 480 IF ( k5 <= nzb_s_inner(j5,i3) .AND. & 481 nzb_s_inner(j5,i3) /= 0 ) THEN 482 483 IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. & 484 k3 == nzb_s_inner(j5,i3) ) THEN 485 reflect_z = .TRUE. 486 EXIT case3 487 ENDIF 488 489 ENDIF 490 491 IF ( k3 <= nzb_s_inner(j3,i3) .AND. & 492 nzb_s_inner(j3,i3) /= 0 ) THEN 493 494 IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. & 495 k3 == nzb_s_inner(j3,i3) ) THEN 496 reflect_z = .TRUE. 497 EXIT case3 498 ENDIF 499 500 IF ( pos_y == ( j3 * dy - 0.5 * dy ) .AND. & 501 pos_z < nzb_s_inner(j3,i3) * dz ) THEN 502 reflect_y = .TRUE. 503 EXIT case3 504 ENDIF 505 506 ENDIF 507 508 IF ( k5 <= nzb_s_inner(j3,i5) .AND. & 509 nzb_s_inner(j3,i5) /= 0 ) THEN 510 511 IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. & 512 k3 == nzb_s_inner(j3,i5) ) THEN 513 reflect_z = .TRUE. 514 EXIT case3 515 ENDIF 516 517 IF ( pos_x == ( i5 * dx + 0.5 * dx ) .AND. & 518 pos_z < nzb_s_inner(j3,i5) * dz ) THEN 519 reflect_x = .TRUE. 520 EXIT case3 521 ENDIF 522 523 ENDIF 524 525 ENDDO case3 526 527 ! 528 !-- Case 4 529 ELSEIF ( particles(n)%x < pos_x_old .AND. & 530 particles(n)%y < pos_y_old ) THEN 531 532 t_index = 1 533 534 DO i = i1, i2, -1 535 xline = i * dx - 0.5 * dx 536 t(t_index) = ( pos_x_old - xline ) / & 537 ( pos_x_old - particles(n)%x ) 538 t_index = t_index + 1 539 ENDDO 540 541 DO j = j1, j2, -1 542 yline = j * dy - 0.5 * dy 543 t(t_index) = ( pos_y_old - yline ) / & 544 ( pos_y_old - particles(n)%y ) 545 t_index = t_index + 1 546 ENDDO 547 548 IF ( particles(n)%z < pos_z_old ) THEN 549 DO k = k1, k2 , -1 550 zline = k * dz 551 t(t_index) = ( pos_z_old - zline ) / & 552 ( pos_z_old-particles(n)%z ) 553 t_index = t_index + 1 554 ENDDO 555 ENDIF 556 t_index_number = t_index-1 557 558 ! 559 !-- Sorting t 560 inc = 1 561 jr = 1 562 563 DO WHILE ( inc <= t_index_number ) 564 inc = 3 * inc + 1 565 ENDDO 566 567 DO WHILE ( inc > 1 ) 568 inc = inc / 3 569 DO ir = inc+1, t_index_number 570 tmp_t = t(ir) 571 jr = ir 572 DO WHILE ( t(jr-inc) > tmp_t ) 573 t(jr) = t(jr-inc) 574 jr = jr - inc 575 IF ( jr <= inc ) EXIT 576 ENDDO 577 t(jr) = tmp_t 578 ENDDO 579 ENDDO 580 581 case4: DO t_index = 1, t_index_number 582 583 pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old ) 584 pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old ) 585 pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old ) 586 587 i3 = ( pos_x + 0.5 * dx ) * ddx 588 j3 = ( pos_y + 0.5 * dy ) * ddy 589 k3 = pos_z / dz + offset_ocean_nzt_m1 590 591 i5 = pos_x * ddx 592 j5 = pos_y * ddy 593 k5 = pos_z / dz + offset_ocean_nzt_m1 594 595 IF ( k3 <= nzb_s_inner(j3,i3) .AND. & 596 nzb_s_inner(j3,i3) /= 0 ) THEN 597 598 IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. & 599 k3 == nzb_s_inner(j3,i3) ) THEN 600 reflect_z = .TRUE. 601 EXIT case4 602 ENDIF 603 604 ENDIF 605 606 IF ( k5 <= nzb_s_inner(j3,i5) .AND. & 607 nzb_s_inner(j3,i5) /= 0 ) THEN 608 609 IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. & 610 k3 == nzb_s_inner(j3,i5) ) THEN 611 reflect_z = .TRUE. 612 EXIT case4 613 ENDIF 614 615 IF ( pos_x == ( i5 * dx + 0.5 * dx ) .AND. & 616 nzb_s_inner(j3,i5) /=0 .AND. & 617 pos_z < nzb_s_inner(j3,i5) * dz ) THEN 618 reflect_x = .TRUE. 619 EXIT case4 620 ENDIF 621 622 ENDIF 623 624 IF ( k5 <= nzb_s_inner(j5,i3) .AND. & 625 nzb_s_inner(j5,i3) /= 0 ) THEN 626 627 IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. & 628 k5 == nzb_s_inner(j5,i3) ) THEN 629 reflect_z = .TRUE. 630 EXIT case4 631 ENDIF 632 633 IF ( pos_y == ( j5 * dy + 0.5 * dy ) .AND. & 634 nzb_s_inner(j5,i3) /= 0 .AND. & 635 pos_z < nzb_s_inner(j5,i3) * dz ) THEN 636 reflect_y = .TRUE. 637 EXIT case4 638 ENDIF 639 640 ENDIF 538 641 539 ENDDO case4540 541 ENDIF542 543 ENDIF ! Check, if particle position is below the surface544 545 ! 546 !-- Do the respective reflection, in case that one of the above conditions547 !-- is found to be true548 IF ( reflect_z ) THEN549 550 particles(n)%z = 2.0 * pos_z - prt_z551 particles(n)%speed_z = - particles(n)%speed_z552 553 IF ( use_sgs_for_particles ) THEN554 particles(n)%rvar3 = - particles(n)%rvar3555 ENDIF556 reflect_z = .FALSE.557 558 ELSEIF ( reflect_y ) THEN559 560 particles(n)%y = 2.0 * pos_y - prt_y561 particles(n)%speed_y = - particles(n)%speed_y562 563 IF ( use_sgs_for_particles ) THEN564 particles(n)%rvar2 = - particles(n)%rvar2565 ENDIF566 reflect_y = .FALSE.567 568 ELSEIF ( reflect_x ) THEN569 570 particles(n)%x = 2.0 * pos_x - prt_x571 particles(n)%speed_x = - particles(n)%speed_x572 573 IF ( use_sgs_for_particles ) THEN574 particles(n)%rvar1 = - particles(n)%rvar1575 ENDIF576 reflect_x = .FALSE.577 578 ENDIF642 ENDDO case4 643 644 ENDIF 645 646 ENDIF ! Check, if particle position is below the surface 647 648 ! 649 !-- Do the respective reflection, in case that one of the above conditions 650 !-- is found to be true 651 IF ( reflect_z ) THEN 652 653 particles(n)%z = 2.0 * pos_z - prt_z 654 particles(n)%speed_z = - particles(n)%speed_z 655 656 IF ( use_sgs_for_particles ) THEN 657 particles(n)%rvar3 = - particles(n)%rvar3 658 ENDIF 659 reflect_z = .FALSE. 660 661 ELSEIF ( reflect_y ) THEN 662 663 particles(n)%y = 2.0 * pos_y - prt_y 664 particles(n)%speed_y = - particles(n)%speed_y 665 666 IF ( use_sgs_for_particles ) THEN 667 particles(n)%rvar2 = - particles(n)%rvar2 668 ENDIF 669 reflect_y = .FALSE. 670 671 ELSEIF ( reflect_x ) THEN 672 673 particles(n)%x = 2.0 * pos_x - prt_x 674 particles(n)%speed_x = - particles(n)%speed_x 675 676 IF ( use_sgs_for_particles ) THEN 677 particles(n)%rvar1 = - particles(n)%rvar1 678 ENDIF 679 reflect_x = .FALSE. 680 681 ENDIF 579 682 580 ENDDO 581 582 CALL cpu_log( log_point_s(48), 'advec_part_refle', 'stop' ) 583 584 585 END SUBROUTINE particle_boundary_conds 683 ENDDO 684 685 CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'stop' ) 686 687 ENDIF 688 689 END SUBROUTINE lpm_boundary_conds
Note: See TracChangeset
for help on using the changeset viewer.
