Changeset 1929 for palm/trunk/SOURCE/lpm_boundary_conds.f90
- Timestamp:
- Jun 9, 2016 4:25:25 PM (8 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/lpm_boundary_conds.f90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lpm_boundary_conds.f90
r1823 r1929 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Rewritten wall reflection 22 22 ! 23 23 ! Former revisions: … … 102 102 CHARACTER (LEN=*) :: range !< 103 103 104 INTEGER(iwp) :: i !< 105 INTEGER(iwp) :: inc !< 106 INTEGER(iwp) :: ir !< 107 INTEGER(iwp) :: i1 !< 108 INTEGER(iwp) :: i2 !< 109 INTEGER(iwp) :: i3 !< 110 INTEGER(iwp) :: i5 !< 111 INTEGER(iwp) :: j !< 112 INTEGER(iwp) :: jr !< 113 INTEGER(iwp) :: j1 !< 114 INTEGER(iwp) :: j2 !< 115 INTEGER(iwp) :: j3 !< 116 INTEGER(iwp) :: j5 !< 117 INTEGER(iwp) :: k !< 118 INTEGER(iwp) :: k1 !< 119 INTEGER(iwp) :: k2 !< 120 INTEGER(iwp) :: k3 !< 121 INTEGER(iwp) :: k5 !< 122 INTEGER(iwp) :: n !< 123 INTEGER(iwp) :: t_index !< 124 INTEGER(iwp) :: t_index_number !< 104 INTEGER(iwp) :: inc !< dummy for sorting algorithmus 105 INTEGER(iwp) :: ir !< dummy for sorting algorithmus 106 INTEGER(iwp) :: i1 !< grid index (x) of old particle position 107 INTEGER(iwp) :: i2 !< grid index (x) of current particle position 108 INTEGER(iwp) :: i3 !< grid index (x) of intermediate particle position 109 INTEGER(iwp) :: jr !< dummy for sorting algorithmus 110 INTEGER(iwp) :: j1 !< grid index (y) of old particle position 111 INTEGER(iwp) :: j2 !< grid index (x) of current particle position 112 INTEGER(iwp) :: j3 !< grid index (x) of intermediate particle position 113 INTEGER(iwp) :: n !< particle number 114 INTEGER(iwp) :: t_index !< running index for intermediate particle timesteps in reflection algorithmus 115 INTEGER(iwp) :: t_index_number !< number of intermediate particle timesteps in reflection algorithmus 116 INTEGER(iwp) :: tmp_x !< dummy for sorting algorithmus 117 INTEGER(iwp) :: tmp_y !< dummy for sorting algorithmus 118 119 INTEGER(iwp), DIMENSION(0:10) :: x_ind(0:10) = 0 !< index array (x) of intermediate particle positions 120 INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (x) of intermediate particle positions 125 121 126 LOGICAL :: reflect_x !< 127 LOGICAL :: reflect_y !< 128 LOGICAL :: reflect_z !< 129 130 REAL(wp) :: dt_particle !< 131 REAL(wp) :: pos_x !< 132 REAL(wp) :: pos_x_old !< 133 REAL(wp) :: pos_y !< 134 REAL(wp) :: pos_y_old !< 135 REAL(wp) :: pos_z !< 136 REAL(wp) :: pos_z_old !< 137 REAL(wp) :: prt_x !< 138 REAL(wp) :: prt_y !< 139 REAL(wp) :: prt_z !< 140 REAL(wp) :: t(1:200) !< 141 REAL(wp) :: tmp_t !< 142 REAL(wp) :: xline !< 143 REAL(wp) :: yline !< 144 REAL(wp) :: zline !< 122 LOGICAL :: cross_wall_x !< flag to check if particle reflection along x is necessary 123 LOGICAL :: cross_wall_y !< flag to check if particle reflection along y is necessary 124 LOGICAL :: downwards !< flag to check if particle reflection along z is necessary (only if particle move downwards) 125 LOGICAL :: reflect_x !< flag to check if particle is already reflected along x 126 LOGICAL :: reflect_y !< flag to check if particle is already reflected along y 127 LOGICAL :: reflect_z !< flag to check if particle is already reflected along z 128 LOGICAL :: tmp_reach_x !< dummy for sorting algorithmus 129 LOGICAL :: tmp_reach_y !< dummy for sorting algorithmus 130 LOGICAL :: tmp_reach_z !< dummy for sorting algorithmus 131 LOGICAL :: x_wall_reached !< flag to check if particle has already reached wall 132 LOGICAL :: y_wall_reached !< flag to check if particle has already reached wall 133 134 LOGICAL, DIMENSION(0:10) :: reach_x !< flag to check if particle is at a yz-wall 135 LOGICAL, DIMENSION(0:10) :: reach_y !< flag to check if particle is at a xz-wall 136 LOGICAL, DIMENSION(0:10) :: reach_z !< flag to check if particle is at a xy-wall 137 138 REAL(wp) :: dt_particle !< particle timestep 139 REAL(wp) :: dum !< dummy argument 140 REAL(wp) :: eps = 1E-10_wp !< security number to check if particle has reached a wall 141 REAL(wp) :: pos_x !< intermediate particle position (x) 142 REAL(wp) :: pos_x_old !< particle position (x) at previous particle timestep 143 REAL(wp) :: pos_y !< intermediate particle position (y) 144 REAL(wp) :: pos_y_old !< particle position (y) at previous particle timestep 145 REAL(wp) :: pos_z !< intermediate particle position (z) 146 REAL(wp) :: pos_z_old !< particle position (z) at previous particle timestep 147 REAL(wp) :: prt_x !< current particle position (x) 148 REAL(wp) :: prt_y !< current particle position (y) 149 REAL(wp) :: prt_z !< current particle position (z) 150 REAL(wp) :: t_old !< previous reflection time 151 REAL(wp) :: tmp_t !< dummy for sorting algorithmus 152 REAL(wp) :: xwall !< location of wall in x 153 REAL(wp) :: ywall !< location of wall in y 154 REAL(wp) :: zwall1 !< location of wall in z (old grid box) 155 REAL(wp) :: zwall2 !< location of wall in z (current grid box) 156 REAL(wp) :: zwall3 !< location of wall in z (old y, current x) 157 REAL(wp) :: zwall4 !< location of wall in z (current y, old x) 158 159 REAL(wp), DIMENSION(0:10) :: t !< reflection time 160 145 161 146 162 IF ( range == 'bottom/top' ) THEN … … 211 227 ELSEIF ( range == 'walls' ) THEN 212 228 229 213 230 CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'start' ) 214 231 215 reflect_x = .FALSE.216 reflect_y = .FALSE.217 reflect_z = .FALSE.218 219 232 DO n = 1, number_of_particles 220 233 ! 234 !-- Recalculate particle timestep 221 235 dt_particle = particles(n)%age - particles(n)%age_m 222 236 ! 237 !-- Obtain x/y indices for current particle position 223 238 i2 = ( particles(n)%x + 0.5_wp * dx ) * ddx 224 239 j2 = ( particles(n)%y + 0.5_wp * dy ) * ddy 225 k2 = particles(n)%z / dz + 1 + offset_ocean_nzt_m1 226 240 ! 241 !-- Save current particle positions 227 242 prt_x = particles(n)%x 228 243 prt_y = particles(n)%y 229 244 prt_z = particles(n)%z 230 231 ! 232 !-- If the particle position is below the surface, it has to be reflected 233 IF ( k2 <= nzb_s_inner(j2,i2) .AND. nzb_s_inner(j2,i2) /=0 ) THEN 234 235 pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle 236 pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle 237 pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle 238 i1 = ( pos_x_old + 0.5_wp * dx ) * ddx 239 j1 = ( pos_y_old + 0.5_wp * dy ) * ddy 240 k1 = pos_z_old / dz + offset_ocean_nzt_m1 241 242 ! 243 !-- Case 1 244 IF ( particles(n)%x > pos_x_old .AND. particles(n)%y > pos_y_old )& 245 THEN 246 t_index = 1 247 248 DO i = i1, i2 249 xline = i * dx + 0.5_wp * dx 250 t(t_index) = ( xline - pos_x_old ) / & 251 ( particles(n)%x - pos_x_old ) 252 t_index = t_index + 1 245 ! 246 !-- Recalculate old particle positions 247 pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle 248 pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle 249 pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle 250 ! 251 !-- Obtain x/y indices for old particle positions 252 i1 = ( pos_x_old + 0.5_wp * dx ) * ddx 253 j1 = ( pos_y_old + 0.5_wp * dy ) * ddy 254 ! 255 !-- Determine horizontal as well as vertical walls at which particle can 256 !-- be potentially reflected. 257 !-- Start with walls aligned in yz layer. 258 !-- Wall to the right 259 IF ( prt_x > pos_x_old ) THEN 260 xwall = ( i1 + 0.5_wp ) * dx 261 ! 262 !-- Wall to the left 263 ELSE 264 xwall = ( i1 - 0.5_wp ) * dx 265 ENDIF 266 ! 267 !-- Walls aligned in xz layer 268 !-- Wall to the north 269 IF ( prt_y > pos_y_old ) THEN 270 ywall = ( j1 + 0.5_wp ) * dy 271 !-- Wall to the south 272 ELSE 273 ywall = ( j1 - 0.5_wp ) * dy 274 ENDIF 275 ! 276 !-- Walls aligned in xy layer at which particle can be possiblly reflected 277 zwall1 = zw(nzb_s_inner(j2,i2)) 278 zwall2 = zw(nzb_s_inner(j1,i1)) 279 zwall3 = zw(nzb_s_inner(j1,i2)) 280 zwall4 = zw(nzb_s_inner(j2,i1)) 281 ! 282 !-- Initialize flags to check if particle reflection is necessary 283 downwards = .FALSE. 284 cross_wall_x = .FALSE. 285 cross_wall_y = .FALSE. 286 ! 287 !-- Initialize flags to check if a wall is reached 288 reach_x = .FALSE. 289 reach_y = .FALSE. 290 reach_z = .FALSE. 291 ! 292 !-- Initialize flags to check if a particle was already reflected 293 reflect_x = .FALSE. 294 reflect_y = .FALSE. 295 reflect_z = .FALSE. 296 ! 297 !-- Initialize flags to check if a vertical wall is already crossed. 298 !-- ( Required to obtain correct indices. ) 299 x_wall_reached = .FALSE. 300 y_wall_reached = .FALSE. 301 ! 302 !-- Initialize time array 303 t = 0.0_wp 304 ! 305 !-- Check if particle can reach any wall. This case, calculate the 306 !-- fractional time needed to reach this wall. Store this fractional 307 !-- timestep in array t. Moreover, store indices for these grid 308 !-- boxes where the respective wall belongs to. 309 !-- Start with x-direction. 310 t_index = 1 311 t(t_index) = ( xwall - pos_x_old ) & 312 / MERGE( MAX( prt_x - pos_x_old, 1E-30_wp ), & 313 MIN( prt_x - pos_x_old, -1E-30_wp ), & 314 prt_x > pos_x_old ) 315 x_ind(t_index) = i2 316 y_ind(t_index) = j1 317 reach_x(t_index) = .TRUE. 318 reach_y(t_index) = .FALSE. 319 reach_z(t_index) = .FALSE. 320 ! 321 !-- Store these values only if particle really reaches any wall. t must 322 !-- be in a interval between [0:1]. 323 IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) THEN 324 t_index = t_index + 1 325 cross_wall_x = .TRUE. 326 ENDIF 327 ! 328 !-- y-direction 329 t(t_index) = ( ywall - pos_y_old ) & 330 / MERGE( MAX( prt_y - pos_y_old, 1E-30_wp ), & 331 MIN( prt_y - pos_y_old, -1E-30_wp ), & 332 prt_y > pos_y_old ) 333 x_ind(t_index) = i1 334 y_ind(t_index) = j2 335 reach_x(t_index) = .FALSE. 336 reach_y(t_index) = .TRUE. 337 reach_z(t_index) = .FALSE. 338 IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) THEN 339 t_index = t_index + 1 340 cross_wall_y = .TRUE. 341 ENDIF 342 ! 343 !-- z-direction 344 !-- At first, check if particle moves downwards. Only in this case a 345 !-- particle can be reflected vertically. 346 IF ( prt_z < pos_z_old ) THEN 347 348 downwards = .TRUE. 349 dum = 1.0_wp / MERGE( MAX( prt_z - pos_z_old, 1E-30_wp ), & 350 MIN( prt_z - pos_z_old, -1E-30_wp ), & 351 prt_z > pos_z_old ) 352 353 t(t_index) = ( zwall1 - pos_z_old ) * dum 354 x_ind(t_index) = i2 355 y_ind(t_index) = j2 356 reach_x(t_index) = .FALSE. 357 reach_y(t_index) = .FALSE. 358 reach_z(t_index) = .TRUE. 359 IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) & 360 t_index = t_index + 1 361 362 reach_x(t_index) = .FALSE. 363 reach_y(t_index) = .FALSE. 364 reach_z(t_index) = .TRUE. 365 t(t_index) = ( zwall2 - pos_z_old ) * dum 366 x_ind(t_index) = i1 367 y_ind(t_index) = j1 368 IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) & 369 t_index = t_index + 1 370 371 reach_x(t_index) = .FALSE. 372 reach_y(t_index) = .FALSE. 373 reach_z(t_index) = .TRUE. 374 t(t_index) = ( zwall3 - pos_z_old ) * dum 375 x_ind(t_index) = i2 376 y_ind(t_index) = j1 377 IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) & 378 t_index = t_index + 1 379 380 reach_x(t_index) = .FALSE. 381 reach_y(t_index) = .FALSE. 382 reach_z(t_index) = .TRUE. 383 t(t_index) = ( zwall4 - pos_z_old ) * dum 384 x_ind(t_index) = i1 385 y_ind(t_index) = j2 386 IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) & 387 t_index = t_index + 1 388 389 END IF 390 t_index_number = t_index - 1 391 ! 392 !-- Carry out reflection only if particle reaches any wall 393 IF ( cross_wall_x .OR. cross_wall_y .OR. downwards ) THEN 394 ! 395 !-- Sort fractional timesteps in ascending order. Also sort the 396 !-- corresponding indices and flag according to the time interval a 397 !-- particle reaches the respective wall. 398 inc = 1 399 jr = 1 400 DO WHILE ( inc <= t_index_number ) 401 inc = 3 * inc + 1 402 ENDDO 403 404 DO WHILE ( inc > 1 ) 405 inc = inc / 3 406 DO ir = inc+1, t_index_number 407 tmp_t = t(ir) 408 tmp_x = x_ind(ir) 409 tmp_y = y_ind(ir) 410 tmp_reach_x = reach_x(ir) 411 tmp_reach_y = reach_y(ir) 412 tmp_reach_z = reach_z(ir) 413 jr = ir 414 DO WHILE ( t(jr-inc) > tmp_t ) 415 t(jr) = t(jr-inc) 416 x_ind(jr) = x_ind(jr-inc) 417 y_ind(jr) = y_ind(jr-inc) 418 reach_x(jr) = reach_x(jr-inc) 419 reach_y(jr) = reach_y(jr-inc) 420 reach_z(jr) = reach_z(jr-inc) 421 jr = jr - inc 422 IF ( jr <= inc ) EXIT 423 ENDDO 424 t(jr) = tmp_t 425 x_ind(jr) = tmp_x 426 y_ind(jr) = tmp_y 427 reach_x(jr) = tmp_reach_x 428 reach_y(jr) = tmp_reach_y 429 reach_z(jr) = tmp_reach_z 253 430 ENDDO 254 255 DO j = j1, j2 256 yline = j * dy + 0.5_wp * dy 257 t(t_index) = ( yline - pos_y_old ) / & 258 ( particles(n)%y - pos_y_old ) 259 t_index = t_index + 1 260 ENDDO 261 262 IF ( particles(n)%z < pos_z_old ) THEN 263 DO k = k1, k2 , -1 264 zline = k * dz 265 t(t_index) = ( pos_z_old - zline ) / & 266 ( pos_z_old - particles(n)%z ) 267 t_index = t_index + 1 268 ENDDO 431 ENDDO 432 ! 433 !-- Initialize temporary particle positions 434 pos_x = pos_x_old 435 pos_y = pos_y_old 436 pos_z = pos_z_old 437 ! 438 !-- Loop over all times a particle possibly moves into a new grid box 439 t_old = 0.0_wp 440 DO t_index = 1, t_index_number 441 ! 442 !-- Calculate intermediate particle position according to the 443 !-- timesteps a particle reaches any wall. 444 pos_x = pos_x + ( t(t_index) - t_old ) * dt_particle & 445 * particles(n)%speed_x 446 pos_y = pos_y + ( t(t_index) - t_old ) * dt_particle & 447 * particles(n)%speed_y 448 pos_z = pos_z + ( t(t_index) - t_old ) * dt_particle & 449 * particles(n)%speed_z 450 ! 451 !-- Obtain x/y grid indices for intermediate particle position from 452 !-- sorted index array 453 i3 = x_ind(t_index) 454 j3 = y_ind(t_index) 455 ! 456 !-- Check which wall is already reached 457 IF ( .NOT. x_wall_reached ) x_wall_reached = reach_x(t_index) 458 IF ( .NOT. y_wall_reached ) y_wall_reached = reach_y(t_index) 459 ! 460 !-- Check if a particle needs to be reflected at any yz-wall. If 461 !-- necessary, carry out reflection. Please note, a security 462 !-- constant is required, as the particle position do not 463 !-- necessarily exactly match the wall location due to rounding 464 !-- errors. 465 IF ( ABS( pos_x - xwall ) < eps .AND. & 466 pos_z <= zw(nzb_s_inner(j3,i3)) .AND. & 467 reach_x(t_index) .AND. & 468 .NOT. reflect_x ) THEN 469 ! 470 !-- Reflection in x-direction. 471 !-- Ensure correct reflection by MIN/MAX functions, depending on 472 !-- direction of particle transport. 473 !-- Due to rounding errors pos_x do not exactly matches the wall 474 !-- location, leading to erroneous reflection. 475 pos_x = MERGE( MIN( 2.0_wp * xwall - pos_x, xwall ), & 476 MAX( 2.0_wp * xwall - pos_x, xwall ), & 477 particles(n)%x > xwall ) 478 ! 479 !-- Change sign of particle speed 480 particles(n)%speed_x = - particles(n)%speed_x 481 ! 482 !-- Change also sign of subgrid-scale particle speed 483 particles(n)%rvar1 = - particles(n)%rvar1 484 ! 485 !-- Set flag that reflection along x is already done 486 reflect_x = .TRUE. 487 ! 488 !-- As particle do not crosses any further yz-wall during 489 !-- this timestep, set further x-indices to the current one. 490 x_ind(t_index:t_index_number) = i1 491 ! 492 !-- If particle already reached the wall but was not reflected, 493 !-- set further x-indices to the new one. 494 ELSEIF ( x_wall_reached .AND. .NOT. reflect_x ) THEN 495 x_ind(t_index:t_index_number) = i2 269 496 ENDIF 270 271 t_index_number = t_index - 1 272 273 ! 274 !-- Sorting t 275 inc = 1 276 jr = 1 277 DO WHILE ( inc <= t_index_number ) 278 inc = 3 * inc + 1 279 ENDDO 280 281 DO WHILE ( inc > 1 ) 282 inc = inc / 3 283 DO ir = inc+1, t_index_number 284 tmp_t = t(ir) 285 jr = ir 286 DO WHILE ( t(jr-inc) > tmp_t ) 287 t(jr) = t(jr-inc) 288 jr = jr - inc 289 IF ( jr <= inc ) EXIT 290 ENDDO 291 t(jr) = tmp_t 292 ENDDO 293 ENDDO 294 295 case1: DO t_index = 1, t_index_number 296 297 pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old ) 298 pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old ) 299 pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old ) 300 301 i3 = ( pos_x + 0.5_wp * dx ) * ddx 302 j3 = ( pos_y + 0.5_wp * dy ) * ddy 303 k3 = pos_z / dz + offset_ocean_nzt_m1 304 305 i5 = pos_x * ddx 306 j5 = pos_y * ddy 307 k5 = pos_z / dz + offset_ocean_nzt_m1 308 309 IF ( k5 <= nzb_s_inner(j5,i3) .AND. & 310 nzb_s_inner(j5,i3) /= 0 ) THEN 311 312 IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. & 313 k3 == nzb_s_inner(j5,i3) ) THEN 314 reflect_z = .TRUE. 315 EXIT case1 316 ENDIF 497 ! 498 !-- Check if a particle needs to be reflected at any xz-wall. If 499 !-- necessary, carry out reflection. 500 IF ( ABS( pos_y - ywall ) < eps .AND. & 501 pos_z <= zw(nzb_s_inner(j3,i3)) .AND. & 502 reach_y(t_index) .AND. & 503 .NOT. reflect_y ) THEN 504 505 pos_y = MERGE( MIN( 2.0_wp * ywall - pos_y, ywall ), & 506 MAX( 2.0_wp * ywall - pos_y, ywall ), & 507 particles(n)%y > ywall ) 508 509 particles(n)%speed_y = - particles(n)%speed_y 510 particles(n)%rvar2 = - particles(n)%rvar2 511 512 reflect_y = .TRUE. 513 y_ind(t_index:t_index_number) = j1 514 515 ELSEIF ( y_wall_reached .AND. .NOT. reflect_y ) THEN 516 y_ind(t_index:t_index_number) = j2 517 ENDIF 518 ! 519 !-- Check if a particle needs to be reflected at any xy-wall. If 520 !-- necessary, carry out reflection. 521 IF ( downwards .AND. reach_z(t_index) .AND. & 522 .NOT. reflect_z ) THEN 523 IF ( pos_z - zw(nzb_s_inner(j3,i3)) < eps ) THEN 524 525 pos_z = MAX( 2.0_wp * zw(nzb_s_inner(j3,i3)) - pos_z, & 526 zw(nzb_s_inner(j3,i3)) ) 527 528 particles(n)%speed_z = - particles(n)%speed_z 529 particles(n)%rvar3 = - particles(n)%rvar3 530 531 reflect_z = .TRUE. 317 532 318 533 ENDIF 319 534 320 IF ( k5 <= nzb_s_inner(j3,i5) .AND. &321 nzb_s_inner(j3,i5) /= 0 ) THEN322 323 IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. &324 k3 == nzb_s_inner(j3,i5) ) THEN325 reflect_z = .TRUE.326 EXIT case1327 ENDIF328 329 ENDIF330 331 IF ( k3 <= nzb_s_inner(j3,i3) .AND. &332 nzb_s_inner(j3,i3) /= 0 ) THEN333 334 IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. &335 k3 == nzb_s_inner(j3,i3) ) THEN336 reflect_z = .TRUE.337 EXIT case1338 ENDIF339 340 IF ( pos_y == ( j3 * dy - 0.5_wp * dy ) .AND. &341 pos_z < nzb_s_inner(j3,i3) * dz ) THEN342 reflect_y = .TRUE.343 EXIT case1344 ENDIF345 346 IF ( pos_x == ( i3 * dx - 0.5_wp * dx ) .AND. &347 pos_z < nzb_s_inner(j3,i3) * dz ) THEN348 reflect_x = .TRUE.349 EXIT case1350 ENDIF351 352 ENDIF353 354 ENDDO case1355 356 !357 !-- Case 2358 ELSEIF ( particles(n)%x > pos_x_old .AND. &359 particles(n)%y < pos_y_old ) THEN360 361 t_index = 1362 363 DO i = i1, i2364 xline = i * dx + 0.5_wp * dx365 t(t_index) = ( xline - pos_x_old ) / &366 ( particles(n)%x - pos_x_old )367 t_index = t_index + 1368 ENDDO369 370 DO j = j1, j2, -1371 yline = j * dy - 0.5_wp * dy372 t(t_index) = ( pos_y_old - yline ) / &373 ( pos_y_old - particles(n)%y )374 t_index = t_index + 1375 ENDDO376 377 IF ( particles(n)%z < pos_z_old ) THEN378 DO k = k1, k2 , -1379 zline = k * dz380 t(t_index) = ( pos_z_old - zline ) / &381 ( pos_z_old - particles(n)%z )382 t_index = t_index + 1383 ENDDO384 535 ENDIF 385 t_index_number = t_index-1 386 387 ! 388 !-- Sorting t 389 inc = 1 390 jr = 1 391 DO WHILE ( inc <= t_index_number ) 392 inc = 3 * inc + 1 393 ENDDO 394 395 DO WHILE ( inc > 1 ) 396 inc = inc / 3 397 DO ir = inc+1, t_index_number 398 tmp_t = t(ir) 399 jr = ir 400 DO WHILE ( t(jr-inc) > tmp_t ) 401 t(jr) = t(jr-inc) 402 jr = jr - inc 403 IF ( jr <= inc ) EXIT 404 ENDDO 405 t(jr) = tmp_t 406 ENDDO 407 ENDDO 408 409 case2: DO t_index = 1, t_index_number 410 411 pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old ) 412 pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old ) 413 pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old ) 414 415 i3 = ( pos_x + 0.5_wp * dx ) * ddx 416 j3 = ( pos_y + 0.5_wp * dy ) * ddy 417 k3 = pos_z / dz + offset_ocean_nzt_m1 418 419 i5 = pos_x * ddx 420 j5 = pos_y * ddy 421 k5 = pos_z / dz + offset_ocean_nzt_m1 422 423 IF ( k5 <= nzb_s_inner(j3,i5) .AND. & 424 nzb_s_inner(j3,i5) /= 0 ) THEN 425 426 IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. & 427 k3 == nzb_s_inner(j3,i5) ) THEN 428 reflect_z = .TRUE. 429 EXIT case2 430 ENDIF 431 432 ENDIF 433 434 IF ( k3 <= nzb_s_inner(j3,i3) .AND. & 435 nzb_s_inner(j3,i3) /= 0 ) THEN 436 437 IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. & 438 k3 == nzb_s_inner(j3,i3) ) THEN 439 reflect_z = .TRUE. 440 EXIT case2 441 ENDIF 442 443 IF ( pos_x == ( i3 * dx - 0.5_wp * dx ) .AND. & 444 pos_z < nzb_s_inner(j3,i3) * dz ) THEN 445 reflect_x = .TRUE. 446 EXIT case2 447 ENDIF 448 449 ENDIF 450 451 IF ( k5 <= nzb_s_inner(j5,i3) .AND. & 452 nzb_s_inner(j5,i3) /= 0 ) THEN 453 454 IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. & 455 k3 == nzb_s_inner(j5,i3) ) THEN 456 reflect_z = .TRUE. 457 EXIT case2 458 ENDIF 459 460 IF ( pos_y == ( j5 * dy + 0.5_wp * dy ) .AND. & 461 pos_z < nzb_s_inner(j5,i3) * dz ) THEN 462 reflect_y = .TRUE. 463 EXIT case2 464 ENDIF 465 466 ENDIF 467 468 ENDDO case2 469 470 ! 471 !-- Case 3 472 ELSEIF ( particles(n)%x < pos_x_old .AND. & 473 particles(n)%y > pos_y_old ) THEN 474 475 t_index = 1 476 477 DO i = i1, i2, -1 478 xline = i * dx - 0.5_wp * dx 479 t(t_index) = ( pos_x_old - xline ) / & 480 ( pos_x_old - particles(n)%x ) 481 t_index = t_index + 1 482 ENDDO 483 484 DO j = j1, j2 485 yline = j * dy + 0.5_wp * dy 486 t(t_index) = ( yline - pos_y_old ) / & 487 ( particles(n)%y - pos_y_old ) 488 t_index = t_index + 1 489 ENDDO 490 491 IF ( particles(n)%z < pos_z_old ) THEN 492 DO k = k1, k2 , -1 493 zline = k * dz 494 t(t_index) = ( pos_z_old - zline ) / & 495 ( pos_z_old - particles(n)%z ) 496 t_index = t_index + 1 497 ENDDO 498 ENDIF 499 t_index_number = t_index - 1 500 501 ! 502 !-- Sorting t 503 inc = 1 504 jr = 1 505 506 DO WHILE ( inc <= t_index_number ) 507 inc = 3 * inc + 1 508 ENDDO 509 510 DO WHILE ( inc > 1 ) 511 inc = inc / 3 512 DO ir = inc+1, t_index_number 513 tmp_t = t(ir) 514 jr = ir 515 DO WHILE ( t(jr-inc) > tmp_t ) 516 t(jr) = t(jr-inc) 517 jr = jr - inc 518 IF ( jr <= inc ) EXIT 519 ENDDO 520 t(jr) = tmp_t 521 ENDDO 522 ENDDO 523 524 case3: DO t_index = 1, t_index_number 525 526 pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old ) 527 pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old ) 528 pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old ) 529 530 i3 = ( pos_x + 0.5_wp * dx ) * ddx 531 j3 = ( pos_y + 0.5_wp * dy ) * ddy 532 k3 = pos_z / dz + offset_ocean_nzt_m1 533 534 i5 = pos_x * ddx 535 j5 = pos_y * ddy 536 k5 = pos_z / dz + offset_ocean_nzt_m1 537 538 IF ( k5 <= nzb_s_inner(j5,i3) .AND. & 539 nzb_s_inner(j5,i3) /= 0 ) THEN 540 541 IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. & 542 k3 == nzb_s_inner(j5,i3) ) THEN 543 reflect_z = .TRUE. 544 EXIT case3 545 ENDIF 546 547 ENDIF 548 549 IF ( k3 <= nzb_s_inner(j3,i3) .AND. & 550 nzb_s_inner(j3,i3) /= 0 ) THEN 551 552 IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. & 553 k3 == nzb_s_inner(j3,i3) ) THEN 554 reflect_z = .TRUE. 555 EXIT case3 556 ENDIF 557 558 IF ( pos_y == ( j3 * dy - 0.5_wp * dy ) .AND. & 559 pos_z < nzb_s_inner(j3,i3) * dz ) THEN 560 reflect_y = .TRUE. 561 EXIT case3 562 ENDIF 563 564 ENDIF 565 566 IF ( k5 <= nzb_s_inner(j3,i5) .AND. & 567 nzb_s_inner(j3,i5) /= 0 ) THEN 568 569 IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. & 570 k3 == nzb_s_inner(j3,i5) ) THEN 571 reflect_z = .TRUE. 572 EXIT case3 573 ENDIF 574 575 IF ( pos_x == ( i5 * dx + 0.5_wp * dx ) .AND. & 576 pos_z < nzb_s_inner(j3,i5) * dz ) THEN 577 reflect_x = .TRUE. 578 EXIT case3 579 ENDIF 580 581 ENDIF 582 583 ENDDO case3 584 585 ! 586 !-- Case 4 587 ELSEIF ( particles(n)%x < pos_x_old .AND. & 588 particles(n)%y < pos_y_old ) THEN 589 590 t_index = 1 591 592 DO i = i1, i2, -1 593 xline = i * dx - 0.5_wp * dx 594 t(t_index) = ( pos_x_old - xline ) / & 595 ( pos_x_old - particles(n)%x ) 596 t_index = t_index + 1 597 ENDDO 598 599 DO j = j1, j2, -1 600 yline = j * dy - 0.5_wp * dy 601 t(t_index) = ( pos_y_old - yline ) / & 602 ( pos_y_old - particles(n)%y ) 603 t_index = t_index + 1 604 ENDDO 605 606 IF ( particles(n)%z < pos_z_old ) THEN 607 DO k = k1, k2 , -1 608 zline = k * dz 609 t(t_index) = ( pos_z_old - zline ) / & 610 ( pos_z_old-particles(n)%z ) 611 t_index = t_index + 1 612 ENDDO 613 ENDIF 614 t_index_number = t_index-1 615 616 ! 617 !-- Sorting t 618 inc = 1 619 jr = 1 620 621 DO WHILE ( inc <= t_index_number ) 622 inc = 3 * inc + 1 623 ENDDO 624 625 DO WHILE ( inc > 1 ) 626 inc = inc / 3 627 DO ir = inc+1, t_index_number 628 tmp_t = t(ir) 629 jr = ir 630 DO WHILE ( t(jr-inc) > tmp_t ) 631 t(jr) = t(jr-inc) 632 jr = jr - inc 633 IF ( jr <= inc ) EXIT 634 ENDDO 635 t(jr) = tmp_t 636 ENDDO 637 ENDDO 638 639 case4: DO t_index = 1, t_index_number 640 641 pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old ) 642 pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old ) 643 pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old ) 644 645 i3 = ( pos_x + 0.5_wp * dx ) * ddx 646 j3 = ( pos_y + 0.5_wp * dy ) * ddy 647 k3 = pos_z / dz + offset_ocean_nzt_m1 648 649 i5 = pos_x * ddx 650 j5 = pos_y * ddy 651 k5 = pos_z / dz + offset_ocean_nzt_m1 652 653 IF ( k3 <= nzb_s_inner(j3,i3) .AND. & 654 nzb_s_inner(j3,i3) /= 0 ) THEN 655 656 IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. & 657 k3 == nzb_s_inner(j3,i3) ) THEN 658 reflect_z = .TRUE. 659 EXIT case4 660 ENDIF 661 662 ENDIF 663 664 IF ( k5 <= nzb_s_inner(j3,i5) .AND. & 665 nzb_s_inner(j3,i5) /= 0 ) THEN 666 667 IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. & 668 k3 == nzb_s_inner(j3,i5) ) THEN 669 reflect_z = .TRUE. 670 EXIT case4 671 ENDIF 672 673 IF ( pos_x == ( i5 * dx + 0.5_wp * dx ) .AND. & 674 nzb_s_inner(j3,i5) /=0 .AND. & 675 pos_z < nzb_s_inner(j3,i5) * dz ) THEN 676 reflect_x = .TRUE. 677 EXIT case4 678 ENDIF 679 680 ENDIF 681 682 IF ( k5 <= nzb_s_inner(j5,i3) .AND. & 683 nzb_s_inner(j5,i3) /= 0 ) THEN 684 685 IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. & 686 k5 == nzb_s_inner(j5,i3) ) THEN 687 reflect_z = .TRUE. 688 EXIT case4 689 ENDIF 690 691 IF ( pos_y == ( j5 * dy + 0.5_wp * dy ) .AND. & 692 nzb_s_inner(j5,i3) /= 0 .AND. & 693 pos_z < nzb_s_inner(j5,i3) * dz ) THEN 694 reflect_y = .TRUE. 695 EXIT case4 696 ENDIF 697 698 ENDIF 699 700 ENDDO case4 536 ! 537 !-- Swap time 538 t_old = t(t_index) 539 540 ENDDO 541 ! 542 !-- If a particle was reflected, calculate final position from last 543 !-- intermediate position. 544 IF ( reflect_x .OR. reflect_y .OR. reflect_z ) THEN 545 546 particles(n)%x = pos_x + ( 1.0_wp - t_old ) * dt_particle & 547 * particles(n)%speed_x 548 particles(n)%y = pos_y + ( 1.0_wp - t_old ) * dt_particle & 549 * particles(n)%speed_y 550 particles(n)%z = pos_z + ( 1.0_wp - t_old ) * dt_particle & 551 * particles(n)%speed_z 701 552 702 553 ENDIF 703 554 704 ENDIF ! Check, if particle position is below the surface 705 706 ! 707 !-- Do the respective reflection, in case that one of the above conditions 708 !-- is found to be true 709 IF ( reflect_z ) THEN 710 711 particles(n)%z = 2.0_wp * pos_z - prt_z 712 particles(n)%speed_z = - particles(n)%speed_z 713 714 IF ( use_sgs_for_particles ) THEN 715 particles(n)%rvar3 = - particles(n)%rvar3 716 ENDIF 717 reflect_z = .FALSE. 718 719 ELSEIF ( reflect_y ) THEN 720 721 particles(n)%y = 2.0_wp * pos_y - prt_y 722 particles(n)%speed_y = - particles(n)%speed_y 723 724 IF ( use_sgs_for_particles ) THEN 725 particles(n)%rvar2 = - particles(n)%rvar2 726 ENDIF 727 reflect_y = .FALSE. 728 729 ELSEIF ( reflect_x ) THEN 730 731 particles(n)%x = 2.0_wp * pos_x - prt_x 732 particles(n)%speed_x = - particles(n)%speed_x 733 734 IF ( use_sgs_for_particles ) THEN 735 particles(n)%rvar1 = - particles(n)%rvar1 736 ENDIF 737 reflect_x = .FALSE. 738 739 ENDIF 740 555 ENDIF 556 741 557 ENDDO 742 558
Note: See TracChangeset
for help on using the changeset viewer.
