Changeset 2698 for palm/trunk/SOURCE/lpm_boundary_conds.f90
 Timestamp:
 Dec 14, 2017 6:46:24 PM (4 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/lpm_boundary_conds.f90
r2696 r2698 20 20 ! Current revisions: 21 21 !  22 ! 22 ! Particle reflections at downwardfacing walls implemented. Moreover, 23 ! reflections are adjusted to revised particle grid box location. 24 ! (responsible Philipp Thiele) 23 25 ! 24 26 ! Former revisions: … … 92 94 !> (see offset_ocean_*) 93 95 !! 94 SUBROUTINE lpm_boundary_conds( location )96 SUBROUTINE lpm_boundary_conds( location , i, j, k ) 95 97 96 98 … … 108 110 109 111 USE indices, & 110 ONLY: nxl, nxr, nyn, nys, nz, nzb 112 ONLY: nxl, nxr, nyn, nys, nz, nzb, wall_flags_0,nyng,nysg 111 113 112 114 USE kinds … … 119 121 USE pegrid 120 122 121 USE surface_mod, &122 ONLY: get_topography_top_index123 124 123 IMPLICIT NONE 125 124 126 125 CHARACTER (LEN=*) :: location !< 126 127 INTEGER(iwp), INTENT(IN) :: i !< 128 INTEGER(iwp), INTENT(IN) :: j !< 129 INTEGER(iwp), INTENT(IN) :: k !< 127 130 128 131 INTEGER(iwp) :: inc !< dummy for sorting algorithmus … … 133 136 INTEGER(iwp) :: jr !< dummy for sorting algorithmus 134 137 INTEGER(iwp) :: j1 !< grid index (y) of old particle position 135 INTEGER(iwp) :: j2 !< grid index (x) of current particle position 136 INTEGER(iwp) :: j3 !< grid index (x) of intermediate particle position 137 INTEGER(iwp) :: k_wall !< vertical index of topography top 138 INTEGER(iwp) :: j2 !< grid index (y) of current particle position 139 INTEGER(iwp) :: j3 !< grid index (y) of intermediate particle position 140 INTEGER(iwp) :: k1 !< grid index (z) of old particle position 141 INTEGER(iwp) :: k2 !< grid index (z) of current particle position 142 INTEGER(iwp) :: k3 !< grid index (z) of intermediate particle position 138 143 INTEGER(iwp) :: n !< particle number 139 144 INTEGER(iwp) :: t_index !< running index for intermediate particle timesteps in reflection algorithmus 140 145 INTEGER(iwp) :: t_index_number !< number of intermediate particle timesteps in reflection algorithmus 141 INTEGER(iwp) :: tmp_x !< dummy for sorting algorithmus 142 INTEGER(iwp) :: tmp_y !< dummy for sorting algorithmus 146 INTEGER(iwp) :: tmp_x !< dummy for sorting algorithm 147 INTEGER(iwp) :: tmp_y !< dummy for sorting algorithm 148 INTEGER(iwp) :: tmp_z !< dummy for sorting algorithm 143 149 144 150 INTEGER(iwp), DIMENSION(0:10) :: x_ind(0:10) = 0 !< index array (x) of intermediate particle positions 145 INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (x) of intermediate particle positions 151 INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (y) of intermediate particle positions 152 INTEGER(iwp), DIMENSION(0:10) :: z_ind(0:10) = 0 !< index array (z) of intermediate particle positions 146 153 147 154 LOGICAL :: cross_wall_x !< flag to check if particle reflection along x is necessary 148 155 LOGICAL :: cross_wall_y !< flag to check if particle reflection along y is necessary 149 LOGICAL :: downwards !< flag to check if particle reflection along z is necessary (only if particle move downwards)156 LOGICAL :: cross_wall_z !< flag to check if particle reflection along z is necessary 150 157 LOGICAL :: reflect_x !< flag to check if particle is already reflected along x 151 158 LOGICAL :: reflect_y !< flag to check if particle is already reflected along y … … 156 163 LOGICAL :: x_wall_reached !< flag to check if particle has already reached wall 157 164 LOGICAL :: y_wall_reached !< flag to check if particle has already reached wall 165 LOGICAL :: z_wall_reached !< flag to check if particle has already reached wall 158 166 159 167 LOGICAL, DIMENSION(0:10) :: reach_x !< flag to check if particle is at a yzwall … … 177 185 REAL(wp) :: xwall !< location of wall in x 178 186 REAL(wp) :: ywall !< location of wall in y 179 REAL(wp) :: zwall1 !< location of wall in z (old grid box) 180 REAL(wp) :: zwall2 !< location of wall in z (current grid box) 181 REAL(wp) :: zwall3 !< location of wall in z (old y, current x) 182 REAL(wp) :: zwall4 !< location of wall in z (current y, old x) 187 REAL(wp) :: zwall !< location of wall in z 183 188 184 189 REAL(wp), DIMENSION(0:10) :: t !< reflection time … … 263 268 i2 = particles(n)%x * ddx 264 269 j2 = particles(n)%y * ddy 270 IF (zw(k) < particles(n)%z ) k2 = k + 1 271 IF (zw(k1) > particles(n)%z ) k2 = k  1 265 272 ! 266 273 ! Save current particle positions … … 275 282 ! 276 283 ! Obtain x/y indices for old particle positions 277 i1 = pos_x_old * ddx 278 j1 = pos_y_old * ddy 284 i1 = i 285 j1 = j 286 k1 = k 279 287 ! 280 288 ! Determine horizontal as well as vertical walls at which particle can … … 283 291 ! Wall to the right 284 292 IF ( prt_x > pos_x_old ) THEN 285 xwall = ( i1 + 0.5_wp) * dx293 xwall = ( i1 + 1 ) * dx 286 294 ! 287 295 ! Wall to the left 288 296 ELSE 289 xwall = ( i1  0.5_wp )* dx297 xwall = i1 * dx 290 298 ENDIF 291 299 ! … … 293 301 ! Wall to the north 294 302 IF ( prt_y > pos_y_old ) THEN 295 ywall = ( j1 + 0.5_wp) * dy303 ywall = ( j1 +1 ) * dy 296 304 ! Wall to the south 297 305 ELSE 298 ywall = ( j1  0.5_wp ) * dy 299 ENDIF 300 ! 301 ! Walls aligned in xy layer at which particle can be possiblly reflected. 302 ! The construct of MERGE and BTEST is used to determine the topography 303 ! top index (former nzb_s_inner). 304 zwall1 = zw( get_topography_top_index( j2, i2, 's' ) ) 305 zwall2 = zw( get_topography_top_index( j1, i1, 's' ) ) 306 zwall3 = zw( get_topography_top_index( j1, i2, 's' ) ) 307 zwall4 = zw( get_topography_top_index( j2, i1, 's' ) ) 306 ywall = j1 * dy 307 ENDIF 308 309 IF ( prt_z > pos_z_old ) THEN 310 zwall = zw(k) 311 ELSE 312 zwall = zw(k1) 313 ENDIF 308 314 ! 309 315 ! Initialize flags to check if particle reflection is necessary 310 downwards = .FALSE.311 316 cross_wall_x = .FALSE. 312 317 cross_wall_y = .FALSE. 318 cross_wall_z = .FALSE. 313 319 ! 314 320 ! Initialize flags to check if a wall is reached … … 318 324 ! 319 325 ! Initialize flags to check if a particle was already reflected 320 reflect_x = .FALSE.321 reflect_y = .FALSE.322 reflect_z = .FALSE.323 ! 324 ! Initialize flags to check if a verticalwall is already crossed.326 reflect_x = .FALSE. 327 reflect_y = .FALSE. 328 reflect_z = .FALSE. 329 ! 330 ! Initialize flags to check if a wall is already crossed. 325 331 ! ( Required to obtain correct indices. ) 326 332 x_wall_reached = .FALSE. 327 333 y_wall_reached = .FALSE. 334 z_wall_reached = .FALSE. 328 335 ! 329 336 ! Initialize time array … … 342 349 x_ind(t_index) = i2 343 350 y_ind(t_index) = j1 351 z_ind(t_index) = k1 344 352 reach_x(t_index) = .TRUE. 345 353 reach_y(t_index) = .FALSE. … … 360 368 x_ind(t_index) = i1 361 369 y_ind(t_index) = j2 370 z_ind(t_index) = k1 362 371 reach_x(t_index) = .FALSE. 363 372 reach_y(t_index) = .TRUE. … … 369 378 ! 370 379 ! zdirection 371 ! At first, check if particle moves downwards. Only in this case a 372 ! particle can be reflected vertically. 373 IF ( prt_z < pos_z_old ) THEN 374 375 downwards = .TRUE. 376 dum = 1.0_wp / MERGE( MAX( prt_z  pos_z_old, 1E30_wp ), & 377 MIN( prt_z  pos_z_old, 1E30_wp ), & 378 prt_z > pos_z_old ) 379 380 t(t_index) = ( zwall1  pos_z_old ) * dum 381 x_ind(t_index) = i2 382 y_ind(t_index) = j2 383 reach_x(t_index) = .FALSE. 384 reach_y(t_index) = .FALSE. 385 reach_z(t_index) = .TRUE. 386 IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) & 387 t_index = t_index + 1 388 389 reach_x(t_index) = .FALSE. 390 reach_y(t_index) = .FALSE. 391 reach_z(t_index) = .TRUE. 392 t(t_index) = ( zwall2  pos_z_old ) * dum 393 x_ind(t_index) = i1 394 y_ind(t_index) = j1 395 IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) & 396 t_index = t_index + 1 397 398 reach_x(t_index) = .FALSE. 399 reach_y(t_index) = .FALSE. 400 reach_z(t_index) = .TRUE. 401 t(t_index) = ( zwall3  pos_z_old ) * dum 402 x_ind(t_index) = i2 403 y_ind(t_index) = j1 404 IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) & 405 t_index = t_index + 1 406 407 reach_x(t_index) = .FALSE. 408 reach_y(t_index) = .FALSE. 409 reach_z(t_index) = .TRUE. 410 t(t_index) = ( zwall4  pos_z_old ) * dum 411 x_ind(t_index) = i1 412 y_ind(t_index) = j2 413 IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) & 414 t_index = t_index + 1 415 416 END IF 380 t(t_index) = (zwall  pos_z_old ) & 381 / MERGE( MAX( prt_z  pos_z_old, 1E30_wp ), & 382 MIN( prt_z  pos_z_old, 1E30_wp ), & 383 prt_z > pos_z_old ) 384 385 x_ind(t_index) = i1 386 y_ind(t_index) = j1 387 z_ind(t_index) = k2 388 reach_x(t_index) = .FALSE. 389 reach_y(t_index) = .FALSE. 390 reach_z(t_index) = .TRUE. 391 IF( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp) THEN 392 t_index = t_index + 1 393 cross_wall_z = .TRUE. 394 ENDIF 395 417 396 t_index_number = t_index  1 418 397 ! 419 398 ! Carry out reflection only if particle reaches any wall 420 IF ( cross_wall_x .OR. cross_wall_y .OR. downwards) THEN399 IF ( cross_wall_x .OR. cross_wall_y .OR. cross_wall_z ) THEN 421 400 ! 422 401 ! Sort fractional timesteps in ascending order. Also sort the … … 435 414 tmp_x = x_ind(ir) 436 415 tmp_y = y_ind(ir) 416 tmp_z = z_ind(ir) 437 417 tmp_reach_x = reach_x(ir) 438 418 tmp_reach_y = reach_y(ir) … … 443 423 x_ind(jr) = x_ind(jrinc) 444 424 y_ind(jr) = y_ind(jrinc) 425 z_ind(jr) = z_ind(jrinc) 445 426 reach_x(jr) = reach_x(jrinc) 446 427 reach_y(jr) = reach_y(jrinc) … … 452 433 x_ind(jr) = tmp_x 453 434 y_ind(jr) = tmp_y 435 z_ind(jr) = tmp_z 454 436 reach_x(jr) = tmp_reach_x 455 437 reach_y(jr) = tmp_reach_y … … 480 462 i3 = x_ind(t_index) 481 463 j3 = y_ind(t_index) 464 k3 = z_ind(t_index) 482 465 ! 483 466 ! Check which wall is already reached 484 467 IF ( .NOT. x_wall_reached ) x_wall_reached = reach_x(t_index) 485 IF ( .NOT. y_wall_reached ) y_wall_reached = reach_y(t_index) 468 IF ( .NOT. y_wall_reached ) y_wall_reached = reach_y(t_index) 469 IF ( .NOT. z_wall_reached ) z_wall_reached = reach_z(t_index) 486 470 ! 487 471 ! Check if a particle needs to be reflected at any yzwall. If 488 472 ! necessary, carry out reflection. Please note, a security 489 ! constant is required, as the particle position do not473 ! constant is required, as the particle position does not 490 474 ! necessarily exactly match the wall location due to rounding 491 ! errors. At first, determine index of topography top at (j3,i3) 492 k_wall = get_topography_top_index( j3, i3, 's' ) 493 IF ( ABS( pos_x  xwall ) < eps .AND. & 494 pos_z <= zw(k_wall) .AND. & 495 reach_x(t_index) .AND. & 475 ! errors. 476 IF ( reach_x(t_index) .AND. & 477 ABS( pos_x  xwall ) < eps .AND. & 478 .NOT. BTEST(wall_flags_0(k3,j3,i3),0) .AND. & 496 479 .NOT. reflect_x ) THEN 497 ! 480 ! 481 ! 498 482 ! Reflection in xdirection. 499 483 ! Ensure correct reflection by MIN/MAX functions, depending on 500 484 ! direction of particle transport. 501 ! Due to rounding errors pos_x do not exactly matchesthe wall485 ! Due to rounding errors pos_x does not exactly match the wall 502 486 ! location, leading to erroneous reflection. 503 487 pos_x = MERGE( MIN( 2.0_wp * xwall  pos_x, xwall ), & … … 508 492 particles(n)%speed_x =  particles(n)%speed_x 509 493 ! 510 ! Change alsosign of subgridscale particle speed494 ! Also change sign of subgridscale particle speed 511 495 particles(n)%rvar1 =  particles(n)%rvar1 512 496 ! … … 514 498 reflect_x = .TRUE. 515 499 ! 516 ! As particle do not crosses any further yzwall during500 ! As the particle does not cross any further yzwall during 517 501 ! this timestep, set further xindices to the current one. 518 502 x_ind(t_index:t_index_number) = i1 … … 522 506 ELSEIF ( x_wall_reached .AND. .NOT. reflect_x ) THEN 523 507 x_ind(t_index:t_index_number) = i2 524 ENDIF 508 ENDIF !particle reflection in x direction done 509 525 510 ! 526 511 ! Check if a particle needs to be reflected at any xzwall. If 527 ! necessary, carry out reflection. At first, determine index of 528 ! topography top at (j3,i3) 529 k_wall = get_topography_top_index( j3, i3, 's' ) 530 IF ( ABS( pos_y  ywall ) < eps .AND. & 531 pos_z <= zw(k_wall) .AND. & 532 reach_y(t_index) .AND. & 533 .NOT. reflect_y ) THEN 534 512 ! necessary, carry out reflection. Please note, a security 513 ! constant is required, as the particle position does not 514 ! necessarily exactly match the wall location due to rounding 515 ! errors. 516 WRITE(9,*) 'wall_test_y nysg:',nysg ,' y: ',particles(n)%y,' j3: ',j3,'nyng:', nyng 517 FLUSH(9) 518 WRITE(9,*) BTEST(wall_flags_0(k3,j3,i3),0) 519 FLUSH(9) 520 IF ( reach_y(t_index) .AND. & 521 ABS( pos_y  ywall ) < eps .AND. & 522 .NOT. BTEST(wall_flags_0(k3,j3,i3),0) .AND. & 523 .NOT. reflect_y ) THEN 524 ! 525 ! 526 ! Reflection in ydirection. 527 ! Ensure correct reflection by MIN/MAX functions, depending on 528 ! direction of particle transport. 529 ! Due to rounding errors pos_y does not exactly match the wall 530 ! location, leading to erroneous reflection. 535 531 pos_y = MERGE( MIN( 2.0_wp * ywall  pos_y, ywall ), & 536 532 MAX( 2.0_wp * ywall  pos_y, ywall ), & 537 particles(n)%y > ywall ) 538 539 particles(n)%speed_y =  particles(n)%speed_y 540 particles(n)%rvar2 =  particles(n)%rvar2 541 542 reflect_y = .TRUE. 533 particles(n)%y > ywall ) 534 ! 535 ! Change sign of particle speed 536 particles(n)%speed_y =  particles(n)%speed_y 537 ! 538 ! Also change sign of subgridscale particle speed 539 particles(n)%rvar2 =  particles(n)%rvar2 540 ! 541 ! Set flag that reflection along y is already done 542 reflect_y = .TRUE. 543 ! 544 ! As the particle does not cross any further xzwall during 545 ! this timestep, set further yindices to the current one. 543 546 y_ind(t_index:t_index_number) = j1 544 547 ! 548 ! If particle already reached the wall but was not reflected, 549 ! set further yindices to the new one. 545 550 ELSEIF ( y_wall_reached .AND. .NOT. reflect_y ) THEN 546 y_ind(t_index:t_index_number) = j2 547 ENDIF 551 y_ind(t_index:t_index_number) = j2 552 ENDIF !particle reflection in y direction done 553 548 554 ! 549 555 ! Check if a particle needs to be reflected at any xywall. If 550 ! necessary, carry out reflection. 551 IF ( downwards .AND. reach_z(t_index) .AND. & 556 ! necessary, carry out reflection. Please note, a security 557 ! constant is required, as the particle position does not 558 ! necessarily exactly match the wall location due to rounding 559 ! errors. 560 IF ( reach_z(t_index) .AND. & 561 ABS( pos_z  zwall ) < eps .AND. & 562 .NOT. BTEST(wall_flags_0(k3,j3,i3),0) .AND. & 552 563 .NOT. reflect_z ) THEN 553 ! 554 ! Determine index of topography top at (j3,i3) and chick if 555 ! particle is below. 556 k_wall = get_topography_top_index( j3, i3, 's' ) 557 IF ( pos_z  zw(k_wall) < eps ) THEN 558 559 pos_z = MAX( 2.0_wp * zw(k_wall)  pos_z, & 560 zw(k_wall) ) 561 562 particles(n)%speed_z =  particles(n)%speed_z 563 particles(n)%rvar3 =  particles(n)%rvar3 564 565 reflect_z = .TRUE. 566 567 ENDIF 568 569 ENDIF 564 ! 565 ! 566 ! Reflection in zdirection. 567 ! Ensure correct reflection by MIN/MAX functions, depending on 568 ! direction of particle transport. 569 ! Due to rounding errors pos_z does not exactly match the wall 570 ! location, leading to erroneous reflection. 571 pos_z = MERGE( MIN( 2.0_wp * zwall  pos_z, zwall ), & 572 MAX( 2.0_wp * zwall  pos_z, zwall ), & 573 particles(n)%z > zwall ) 574 ! 575 ! Change sign of particle speed 576 particles(n)%speed_z =  particles(n)%speed_z 577 ! 578 ! Also change sign of subgridscale particle speed 579 particles(n)%rvar3 =  particles(n)%rvar3 580 ! 581 ! Set flag that reflection along z is already done 582 reflect_z = .TRUE. 583 ! 584 ! As the particle does not cross any further xywall during 585 ! this timestep, set further zindices to the current one. 586 z_ind(t_index:t_index_number) = k1 587 ! 588 ! If particle already reached the wall but was not reflected, 589 ! set further zindices to the new one. 590 ELSEIF ( z_wall_reached .AND. .NOT. reflect_z ) THEN 591 z_ind(t_index:t_index_number) = k2 592 ENDIF !particle reflection in z direction done 593 570 594 ! 571 595 ! Swap time
Note: See TracChangeset
for help on using the changeset viewer.