Changeset 1359 for palm/trunk/SOURCE/lpm_droplet_collision.f90
- Timestamp:
- Apr 11, 2014 5:15:14 PM (10 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/lpm_droplet_collision.f90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lpm_droplet_collision.f90
r1323 r1359 1 SUBROUTINE lpm_droplet_collision 1 SUBROUTINE lpm_droplet_collision (i,j,k) 2 2 3 3 !--------------------------------------------------------------------------------! … … 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! New particle structure integrated. 23 ! Kind definition added to all floating point numbers. 23 24 ! 24 25 ! Former revisions: … … 75 76 !------------------------------------------------------------------------------! 76 77 78 77 79 USE arrays_3d, & 78 80 ONLY: diss, ql, ql_v, ql_vp, u, v, w, zu, zw … … 103 105 USE particle_attributes, & 104 106 ONLY: deleted_particles, dissipation_classes, hall_kernel, & 105 palm_kernel, particles, particle_mask, particle_type, & 106 prt_count, prt_start_index, use_kernel_tables, wang_kernel 107 palm_kernel, particles, particle_type, & 108 prt_count, use_kernel_tables, wang_kernel 109 110 USE pegrid 107 111 108 112 IMPLICIT NONE … … 124 128 INTEGER(iwp) :: rclass_s !: 125 129 130 INTEGER(iwp), DIMENSION(prt_count(k,j,i)) :: rclass_v !: 131 132 LOGICAL, SAVE :: first_flag = .TRUE. !: 133 134 TYPE(particle_type) :: tmp_particle !: 135 126 136 REAL(wp) :: aa !: 137 REAL(wp) :: auxn !: temporary variables 138 REAL(wp) :: auxs !: temporary variables 127 139 REAL(wp) :: bb !: 128 140 REAL(wp) :: cc !: … … 158 170 REAL(wp), DIMENSION(:), ALLOCATABLE :: weight !: 159 171 160 161 TYPE(particle_type) :: tmp_particle !:162 163 172 REAL, DIMENSION(prt_count(k,j,i)) :: ck 173 REAL, DIMENSION(prt_count(k,j,i)) :: r3v 174 REAL, DIMENSION(prt_count(k,j,i)) :: sum1v 175 REAL, DIMENSION(prt_count(k,j,i)) :: sum2v 164 176 165 177 CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'start' ) 166 178 167 DO i = nxl, nxr 168 DO j = nys, nyn 169 DO k = nzb+1, nzt 170 ! 171 !-- Collision requires at least two particles in the box 172 IF ( prt_count(k,j,i) > 1 ) THEN 173 ! 174 !-- First, sort particles within the gridbox by their size, 175 !-- using Shell's method (see Numerical Recipes) 176 !-- NOTE: In case of using particle tails, the re-sorting of 177 !-- ---- tails would have to be included here! 178 psi = prt_start_index(k,j,i) - 1 179 inc = 1 180 DO WHILE ( inc <= prt_count(k,j,i) ) 181 inc = 3 * inc + 1 179 ! 180 !-- Collision requires at least two particles in the box 181 IF ( prt_count(k,j,i) > 1 ) THEN 182 ! 183 !-- First, sort particles within the gridbox by their size, 184 !-- using Shell's method (see Numerical Recipes) 185 !-- NOTE: In case of using particle tails, the re-sorting of 186 !-- ---- tails would have to be included here! 187 IF ( .NOT. ( ( hall_kernel .OR. wang_kernel ) .AND. & 188 use_kernel_tables ) ) THEN 189 psi = 0 190 inc = 1 191 DO WHILE ( inc <= prt_count(k,j,i) ) 192 inc = 3 * inc + 1 193 ENDDO 194 195 DO WHILE ( inc > 1 ) 196 inc = inc / 3 197 DO is = inc+1, prt_count(k,j,i) 198 tmp_particle = particles(psi+is) 199 js = is 200 DO WHILE ( particles(psi+js-inc)%radius > & 201 tmp_particle%radius ) 202 particles(psi+js) = particles(psi+js-inc) 203 js = js - inc 204 IF ( js <= inc ) EXIT 182 205 ENDDO 183 184 DO WHILE ( inc > 1 ) 185 inc = inc / 3 186 DO is = inc+1, prt_count(k,j,i) 187 tmp_particle = particles(psi+is) 188 js = is 189 DO WHILE ( particles(psi+js-inc)%radius > & 190 tmp_particle%radius ) 191 particles(psi+js) = particles(psi+js-inc) 192 js = js - inc 193 IF ( js <= inc ) EXIT 194 ENDDO 195 particles(psi+js) = tmp_particle 196 ENDDO 206 particles(psi+js) = tmp_particle 207 ENDDO 208 ENDDO 209 ENDIF 210 211 psi = 1 212 pse = prt_count(k,j,i) 213 214 ! 215 !-- Now apply the different kernels 216 IF ( ( hall_kernel .OR. wang_kernel ) .AND. & 217 use_kernel_tables ) THEN 218 ! 219 !-- Fast method with pre-calculated efficiencies for 220 !-- discrete radius- and dissipation-classes. 221 !-- 222 !-- Determine dissipation class index of this gridbox 223 IF ( wang_kernel ) THEN 224 eclass = INT( diss(k,j,i) * 1.0E4_wp / 1000.0_wp * & 225 dissipation_classes ) + 1 226 epsilon = diss(k,j,i) 227 ELSE 228 epsilon = 0.0_wp 229 ENDIF 230 IF ( hall_kernel .OR. epsilon * 1.0E4_wp < 0.001_wp ) THEN 231 eclass = 0 ! Hall kernel is used 232 ELSE 233 eclass = MIN( dissipation_classes, eclass ) 234 ENDIF 235 236 ! 237 !-- Droplet collision are calculated using collision-coalescence 238 !-- formulation proposed by Wang (see PALM documentation) 239 !-- Since new radii after collision are defined by radii of all 240 !-- droplets before collision, temporary fields for new radii and 241 !-- weighting factors are needed 242 ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i))) 243 244 rad = 0.0_wp 245 weight = 0.0_wp 246 247 sum1v(1:prt_count(k,j,i)) = 0.0_wp 248 sum2v(1:prt_count(k,j,i)) = 0.0_wp 249 250 DO n = 1, prt_count(k,j,i) 251 252 rclass_l = particles(n)%class 253 ! 254 !-- Mass added due to collisions with smaller droplets 255 DO is = n+1, prt_count(k,j,i) 256 rclass_s = particles(is)%class 257 auxs = ckernel(rclass_l,rclass_s,eclass) * particles(is)%weight_factor 258 auxn = ckernel(rclass_s,rclass_l,eclass) * particles(n)%weight_factor 259 IF ( particles(is)%radius < particles(n)%radius ) THEN 260 sum1v(n) = sum1v(n) + particles(is)%radius**3 * auxs 261 sum2v(is) = sum2v(is) + auxn 262 ELSE 263 sum2v(n) = sum2v(n) + auxs 264 sum1v(is) = sum1v(is) + particles(n)%radius**3 * auxn 265 ENDIF 266 ENDDO 267 ENDDO 268 rclass_v = particles(1:prt_count(k,j,i))%class 269 DO n = 1, prt_count(k,j,i) 270 ck(n) = ckernel(rclass_v(n),rclass_v(n),eclass) 271 ENDDO 272 r3v = particles(1:prt_count(k,j,i))%radius**3 273 DO n = 1, prt_count(k,j,i) 274 sum3 = 0.0_wp 275 ddV = ddx * ddy / dz 276 ! 277 !-- Change of the current weighting factor 278 sum3 = 1 - dt_3d * ddV * & 279 ck(n) * & 280 ( particles(n)%weight_factor - 1 ) * 0.5_wp - & 281 dt_3d * ddV * sum2v(n) 282 weight(n) = particles(n)%weight_factor * sum3 283 ! 284 !-- Change of the current droplet radius 285 rad(n) = ( (r3v(n) + dt_3d * ddV * (sum1v(n) - sum2v(n) * r3v(n)) )/& 286 sum3 )**0.33333333333333_wp 287 288 ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n) & 289 * rad(n)**3 290 291 ENDDO 292 IF ( ANY(weight < 0.0_wp) ) THEN 293 WRITE( message_string, * ) 'negative weighting' 294 CALL message( 'lpm_droplet_collision', 'PA0028', & 295 2, 2, -1, 6, 1 ) 296 ENDIF 297 298 particles(psi:pse)%radius = rad(1:prt_count(k,j,i)) 299 particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i)) 300 301 DEALLOCATE(rad, weight) 302 303 ELSEIF ( ( hall_kernel .OR. wang_kernel ) .AND. & 304 .NOT. use_kernel_tables ) THEN 305 ! 306 !-- Collision efficiencies are calculated for every new 307 !-- grid box. First, allocate memory for kernel table. 308 !-- Third dimension is 1, because table is re-calculated for 309 !-- every new dissipation value. 310 ALLOCATE( ckernel(1:prt_count(k,j,i),1:prt_count(k,j,i),1:1) ) 311 ! 312 !-- Now calculate collision efficiencies for this box 313 CALL recalculate_kernel( i, j, k ) 314 315 ! 316 !-- Droplet collision are calculated using collision-coalescence 317 !-- formulation proposed by Wang (see PALM documentation) 318 !-- Since new radii after collision are defined by radii of all 319 !-- droplets before collision, temporary fields for new radii and 320 !-- weighting factors are needed 321 ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i))) 322 323 rad = 0.0_wp 324 weight = 0.0_wp 325 326 DO n = psi, pse, 1 327 328 sum1 = 0.0_wp 329 sum2 = 0.0_wp 330 sum3 = 0.0_wp 331 ! 332 !-- Mass added due to collisions with smaller droplets 333 DO is = psi, n-1 334 sum1 = sum1 + ( particles(is)%radius**3 * & 335 ckernel(n,is,1) * & 336 particles(is)%weight_factor ) 337 ENDDO 338 ! 339 !-- Rate of collisions with larger droplets 340 DO is = n+1, pse 341 sum2 = sum2 + ( ckernel(n,is,1) * & 342 particles(is)%weight_factor ) 343 ENDDO 344 345 r3 = particles(n)%radius**3 346 ddV = ddx * ddy / dz 347 is = 1 348 ! 349 !-- Change of the current weighting factor 350 sum3 = 1 - dt_3d * ddV * & 351 ckernel(n,n,1) * & 352 ( particles(n)%weight_factor - 1 ) * 0.5_wp - & 353 dt_3d * ddV * sum2 354 weight(n-is+1) = particles(n)%weight_factor * sum3 355 ! 356 !-- Change of the current droplet radius 357 rad(n-is+1) = ( (r3 + dt_3d * ddV * (sum1 - sum2 * r3) )/& 358 sum3 )**0.33333333333333_wp 359 360 IF ( weight(n-is+1) < 0.0_wp ) THEN 361 WRITE( message_string, * ) 'negative weighting', & 362 'factor: ', weight(n-is+1) 363 CALL message( 'lpm_droplet_collision', 'PA0037', & 364 2, 2, -1, 6, 1 ) 365 ENDIF 366 367 ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n-is+1) & 368 * rad(n-is+1)**3 369 370 ENDDO 371 372 particles(psi:pse)%radius = rad(1:prt_count(k,j,i)) 373 particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i)) 374 375 DEALLOCATE( rad, weight, ckernel ) 376 377 ELSEIF ( palm_kernel ) THEN 378 ! 379 !-- PALM collision kernel 380 ! 381 !-- Calculate the mean radius of all those particles which 382 !-- are of smaller size than the current particle and 383 !-- use this radius for calculating the collision efficiency 384 DO n = psi+prt_count(k,j,i)-1, psi+1, -1 385 386 sl_r3 = 0.0_wp 387 sl_r4 = 0.0_wp 388 389 DO is = n-1, psi, -1 390 IF ( particles(is)%radius < particles(n)%radius ) & 391 THEN 392 sl_r3 = sl_r3 + particles(is)%weight_factor & 393 * particles(is)%radius**3 394 sl_r4 = sl_r4 + particles(is)%weight_factor & 395 * particles(is)%radius**4 396 ENDIF 397 ENDDO 398 399 IF ( ( sl_r3 ) > 0.0_wp ) THEN 400 mean_r = ( sl_r4 ) / ( sl_r3 ) 401 402 CALL collision_efficiency_rogers( mean_r, & 403 particles(n)%radius, & 404 effective_coll_efficiency ) 405 406 ELSE 407 effective_coll_efficiency = 0.0_wp 408 ENDIF 409 410 IF ( effective_coll_efficiency > 1.0_wp .OR. & 411 effective_coll_efficiency < 0.0_wp ) & 412 THEN 413 WRITE( message_string, * ) 'collision_efficien' , & 414 'cy out of range:' ,effective_coll_efficiency 415 CALL message( 'lpm_droplet_collision', 'PA0145', 2, & 416 2, -1, 6, 1 ) 417 ENDIF 418 419 ! 420 !-- Interpolation of liquid water content 421 ii = particles(n)%x * ddx 422 jj = particles(n)%y * ddy 423 kk = ( particles(n)%z + 0.5_wp * dz ) / dz 424 425 x = particles(n)%x - ii * dx 426 y = particles(n)%y - jj * dy 427 aa = x**2 + y**2 428 bb = ( dx - x )**2 + y**2 429 cc = x**2 + ( dy - y )**2 430 dd = ( dx - x )**2 + ( dy - y )**2 431 gg = aa + bb + cc + dd 432 433 ql_int_l = ( (gg-aa) * ql(kk,jj,ii) + (gg-bb) * & 434 ql(kk,jj,ii+1) & 435 + (gg-cc) * ql(kk,jj+1,ii) + ( gg-dd ) * & 436 ql(kk,jj+1,ii+1) & 437 ) / ( 3.0_wp * gg ) 438 439 ql_int_u = ( (gg-aa) * ql(kk+1,jj,ii) + (gg-bb) * & 440 ql(kk+1,jj,ii+1) & 441 + (gg-cc) * ql(kk+1,jj+1,ii) + (gg-dd) * & 442 ql(kk+1,jj+1,ii+1) & 443 ) / ( 3.0_wp * gg ) 444 445 ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz *& 446 ( ql_int_u - ql_int_l ) 447 448 ! 449 !-- Interpolate u velocity-component 450 ii = ( particles(n)%x + 0.5_wp * dx ) * ddx 451 jj = particles(n)%y * ddy 452 kk = ( particles(n)%z + 0.5_wp * dz ) / dz ! only if equidistant 453 454 IF ( ( particles(n)%z - zu(kk) ) > ( 0.5_wp * dz ) ) kk = kk+1 455 456 x = particles(n)%x + ( 0.5_wp - ii ) * dx 457 y = particles(n)%y - jj * dy 458 aa = x**2 + y**2 459 bb = ( dx - x )**2 + y**2 460 cc = x**2 + ( dy - y )**2 461 dd = ( dx - x )**2 + ( dy - y )**2 462 gg = aa + bb + cc + dd 463 464 u_int_l = ( (gg-aa) * u(kk,jj,ii) + (gg-bb) * & 465 u(kk,jj,ii+1) & 466 + (gg-cc) * u(kk,jj+1,ii) + (gg-dd) * & 467 u(kk,jj+1,ii+1) & 468 ) / ( 3.0_wp * gg ) - u_gtrans 469 IF ( kk+1 == nzt+1 ) THEN 470 u_int = u_int_l 471 ELSE 472 u_int_u = ( (gg-aa) * u(kk+1,jj,ii) + (gg-bb) * & 473 u(kk+1,jj,ii+1) & 474 + (gg-cc) * u(kk+1,jj+1,ii) + (gg-dd) * & 475 u(kk+1,jj+1,ii+1) & 476 ) / ( 3.0_wp * gg ) - u_gtrans 477 u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz & 478 * ( u_int_u - u_int_l ) 479 ENDIF 480 481 ! 482 !-- Same procedure for interpolation of the v velocity-component 483 !-- (adopt index k from u velocity-component) 484 ii = particles(n)%x * ddx 485 jj = ( particles(n)%y + 0.5_wp * dy ) * ddy 486 487 x = particles(n)%x - ii * dx 488 y = particles(n)%y + ( 0.5_wp - jj ) * dy 489 aa = x**2 + y**2 490 bb = ( dx - x )**2 + y**2 491 cc = x**2 + ( dy - y )**2 492 dd = ( dx - x )**2 + ( dy - y )**2 493 gg = aa + bb + cc + dd 494 495 v_int_l = ( ( gg-aa ) * v(kk,jj,ii) + ( gg-bb ) * & 496 v(kk,jj,ii+1) & 497 + ( gg-cc ) * v(kk,jj+1,ii) + ( gg-dd ) * & 498 v(kk,jj+1,ii+1) & 499 ) / ( 3.0_wp * gg ) - v_gtrans 500 IF ( kk+1 == nzt+1 ) THEN 501 v_int = v_int_l 502 ELSE 503 v_int_u = ( (gg-aa) * v(kk+1,jj,ii) + (gg-bb) * & 504 v(kk+1,jj,ii+1) & 505 + (gg-cc) * v(kk+1,jj+1,ii) + (gg-dd) * & 506 v(kk+1,jj+1,ii+1) & 507 ) / ( 3.0_wp * gg ) - v_gtrans 508 v_int = v_int_l + ( particles(n)%z - zu(kk) ) / dz & 509 * ( v_int_u - v_int_l ) 510 ENDIF 511 512 ! 513 !-- Same procedure for interpolation of the w velocity-component 514 !-- (adopt index i from v velocity-component) 515 jj = particles(n)%y * ddy 516 kk = particles(n)%z / dz 517 518 x = particles(n)%x - ii * dx 519 y = particles(n)%y - jj * dy 520 aa = x**2 + y**2 521 bb = ( dx - x )**2 + y**2 522 cc = x**2 + ( dy - y )**2 523 dd = ( dx - x )**2 + ( dy - y )**2 524 gg = aa + bb + cc + dd 525 526 w_int_l = ( ( gg-aa ) * w(kk,jj,ii) + ( gg-bb ) * & 527 w(kk,jj,ii+1) & 528 + ( gg-cc ) * w(kk,jj+1,ii) + ( gg-dd ) * & 529 w(kk,jj+1,ii+1) & 530 ) / ( 3.0_wp * gg ) 531 IF ( kk+1 == nzt+1 ) THEN 532 w_int = w_int_l 533 ELSE 534 w_int_u = ( (gg-aa) * w(kk+1,jj,ii) + (gg-bb) * & 535 w(kk+1,jj,ii+1) & 536 + (gg-cc) * w(kk+1,jj+1,ii) + (gg-dd) * & 537 w(kk+1,jj+1,ii+1) & 538 ) / ( 3.0_wp * gg ) 539 w_int = w_int_l + ( particles(n)%z - zw(kk) ) / dz & 540 * ( w_int_u - w_int_l ) 541 ENDIF 542 543 ! 544 !-- Change in radius due to collision 545 delta_r = effective_coll_efficiency / 3.0_wp & 546 * pi * sl_r3 * ddx * ddy / dz & 547 * SQRT( ( u_int - particles(n)%speed_x )**2 & 548 + ( v_int - particles(n)%speed_y )**2 & 549 + ( w_int - particles(n)%speed_z )**2 & 550 ) * dt_3d 551 ! 552 !-- Change in volume due to collision 553 delta_v = particles(n)%weight_factor & 554 * ( ( particles(n)%radius + delta_r )**3 & 555 - particles(n)%radius**3 ) 556 557 ! 558 !-- Check if collected particles provide enough LWC for 559 !-- volume change of collector particle 560 IF ( delta_v >= sl_r3 .AND. sl_r3 > 0.0_wp ) THEN 561 562 delta_r = ( ( sl_r3/particles(n)%weight_factor ) & 563 + particles(n)%radius**3 )**( 1.0_wp / 3.0_wp ) & 564 - particles(n)%radius 565 566 DO is = n-1, psi, -1 567 IF ( particles(is)%radius < particles(n)%radius ) THEN 568 particles(is)%weight_factor = 0.0_wp 569 particles(is)%particle_mask = .FALSE. 570 deleted_particles = deleted_particles + 1 571 ENDIF 197 572 ENDDO 198 573 199 psi = prt_start_index(k,j,i) 200 pse = psi + prt_count(k,j,i)-1 201 202 ! 203 !-- Now apply the different kernels 204 IF ( ( hall_kernel .OR. wang_kernel ) .AND. & 205 use_kernel_tables ) THEN 206 ! 207 !-- Fast method with pre-calculated efficiencies for 208 !-- discrete radius- and dissipation-classes. 209 ! 210 !-- Determine dissipation class index of this gridbox 211 IF ( wang_kernel ) THEN 212 eclass = INT( diss(k,j,i) * 1.0E4_wp / 1000.0_wp * & 213 dissipation_classes ) + 1 214 epsilon = diss(k,j,i) 215 ELSE 216 epsilon = 0.0 574 ELSE IF ( delta_v < sl_r3 .AND. sl_r3 > 0.0_wp ) THEN 575 576 DO is = n-1, psi, -1 577 IF ( particles(is)%radius < particles(n)%radius & 578 .AND. sl_r3 > 0.0_wp ) THEN 579 particles(is)%weight_factor = & 580 ( ( particles(is)%weight_factor & 581 * ( particles(is)%radius**3 ) ) & 582 - ( delta_v & 583 * particles(is)%weight_factor & 584 * ( particles(is)%radius**3 ) & 585 / sl_r3 ) ) & 586 / ( particles(is)%radius**3 ) 587 588 IF ( particles(is)%weight_factor < 0.0_wp ) THEN 589 WRITE( message_string, * ) 'negative ', & 590 'weighting factor: ', & 591 particles(is)%weight_factor 592 CALL message( 'lpm_droplet_collision', & 593 'PA0039', & 594 2, 2, -1, 6, 1 ) 595 ENDIF 217 596 ENDIF 218 IF ( hall_kernel .OR. epsilon * 1.0E4_wp < 0.001 ) THEN 219 eclass = 0 ! Hall kernel is used 220 ELSE 221 eclass = MIN( dissipation_classes, eclass ) 222 ENDIF 223 224 ! 225 !-- Droplet collision are calculated using collision-coalescence 226 !-- formulation proposed by Wang (see PALM documentation) 227 !-- Since new radii after collision are defined by radii of all 228 !-- droplets before collision, temporary fields for new radii and 229 !-- weighting factors are needed 230 ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i))) 231 232 rad = 0.0 233 weight = 0.0 234 235 DO n = psi, pse, 1 236 237 sum1 = 0.0 238 sum2 = 0.0 239 sum3 = 0.0 240 241 rclass_l = particles(n)%class 242 ! 243 !-- Mass added due to collisions with smaller droplets 244 DO is = psi, n-1 245 246 rclass_s = particles(is)%class 247 sum1 = sum1 + ( particles(is)%radius**3 * & 248 ckernel(rclass_l,rclass_s,eclass) * & 249 particles(is)%weight_factor ) 250 251 ENDDO 252 ! 253 !-- Rate of collisions with larger droplets 254 DO is = n+1, pse 255 256 rclass_s = particles(is)%class 257 sum2 = sum2 + ( ckernel(rclass_l,rclass_s,eclass) * & 258 particles(is)%weight_factor ) 259 260 ENDDO 261 262 r3 = particles(n)%radius**3 263 ddV = ddx * ddy / dz 264 is = prt_start_index(k,j,i) 265 ! 266 !-- Change of the current weighting factor 267 sum3 = 1 - dt_3d * ddV * & 268 ckernel(rclass_l,rclass_l,eclass) * & 269 ( particles(n)%weight_factor - 1 ) * 0.5 - & 270 dt_3d * ddV * sum2 271 weight(n-is+1) = particles(n)%weight_factor * sum3 272 ! 273 !-- Change of the current droplet radius 274 rad(n-is+1) = ( (r3 + dt_3d * ddV * (sum1 - sum2 * r3) )/& 275 sum3 )**0.33333333333333_wp 276 277 IF ( weight(n-is+1) < 0.0 ) THEN 278 WRITE( message_string, * ) 'negative weighting', & 279 'factor: ', weight(n-is+1) 280 CALL message( 'lpm_droplet_collision', 'PA0028', & 281 2, 2, -1, 6, 1 ) 282 ENDIF 283 284 ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n-is+1) & 285 * rad(n-is+1)**3 286 287 ENDDO 288 289 particles(psi:pse)%radius = rad(1:prt_count(k,j,i)) 290 particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i)) 291 292 DEALLOCATE(rad, weight) 293 294 ELSEIF ( ( hall_kernel .OR. wang_kernel ) .AND. & 295 .NOT. use_kernel_tables ) THEN 296 ! 297 !-- Collision efficiencies are calculated for every new 298 !-- grid box. First, allocate memory for kernel table. 299 !-- Third dimension is 1, because table is re-calculated for 300 !-- every new dissipation value. 301 ALLOCATE( ckernel(prt_start_index(k,j,i): & 302 prt_start_index(k,j,i)+prt_count(k,j,i)-1, & 303 prt_start_index(k,j,i): & 304 prt_start_index(k,j,i)+prt_count(k,j,i)-1,1:1) ) 305 ! 306 !-- Now calculate collision efficiencies for this box 307 CALL recalculate_kernel( i, j, k ) 308 309 ! 310 !-- Droplet collision are calculated using collision-coalescence 311 !-- formulation proposed by Wang (see PALM documentation) 312 !-- Since new radii after collision are defined by radii of all 313 !-- droplets before collision, temporary fields for new radii and 314 !-- weighting factors are needed 315 ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i))) 316 317 rad = 0.0 318 weight = 0.0 319 320 DO n = psi, pse, 1 321 322 sum1 = 0.0 323 sum2 = 0.0 324 sum3 = 0.0 325 ! 326 !-- Mass added due to collisions with smaller droplets 327 DO is = psi, n-1 328 sum1 = sum1 + ( particles(is)%radius**3 * & 329 ckernel(n,is,1) * & 330 particles(is)%weight_factor ) 331 ENDDO 332 ! 333 !-- Rate of collisions with larger droplets 334 DO is = n+1, pse 335 sum2 = sum2 + ( ckernel(n,is,1) * & 336 particles(is)%weight_factor ) 337 ENDDO 338 339 r3 = particles(n)%radius**3 340 ddV = ddx * ddy / dz 341 is = prt_start_index(k,j,i) 342 ! 343 !-- Change of the current weighting factor 344 sum3 = 1 - dt_3d * ddV * & 345 ckernel(n,n,1) * & 346 ( particles(n)%weight_factor - 1 ) * 0.5 - & 347 dt_3d * ddV * sum2 348 weight(n-is+1) = particles(n)%weight_factor * sum3 349 ! 350 !-- Change of the current droplet radius 351 rad(n-is+1) = ( (r3 + dt_3d * ddV * (sum1 - sum2 * r3) )/& 352 sum3 )**0.33333333333333_wp 353 354 IF ( weight(n-is+1) < 0.0 ) THEN 355 WRITE( message_string, * ) 'negative weighting', & 356 'factor: ', weight(n-is+1) 357 CALL message( 'lpm_droplet_collision', 'PA0037', & 358 2, 2, -1, 6, 1 ) 359 ENDIF 360 361 ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n-is+1) & 362 * rad(n-is+1)**3 363 364 ENDDO 365 366 particles(psi:pse)%radius = rad(1:prt_count(k,j,i)) 367 particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i)) 368 369 DEALLOCATE( rad, weight, ckernel ) 370 371 ELSEIF ( palm_kernel ) THEN 372 ! 373 !-- PALM collision kernel 374 ! 375 !-- Calculate the mean radius of all those particles which 376 !-- are of smaller size than the current particle and 377 !-- use this radius for calculating the collision efficiency 378 DO n = psi+prt_count(k,j,i)-1, psi+1, -1 379 380 sl_r3 = 0.0 381 sl_r4 = 0.0 382 383 DO is = n-1, psi, -1 384 IF ( particles(is)%radius < particles(n)%radius ) & 385 THEN 386 sl_r3 = sl_r3 + particles(is)%weight_factor & 387 * particles(is)%radius**3 388 sl_r4 = sl_r4 + particles(is)%weight_factor & 389 * particles(is)%radius**4 390 ENDIF 391 ENDDO 392 393 IF ( ( sl_r3 ) > 0.0 ) THEN 394 mean_r = ( sl_r4 ) / ( sl_r3 ) 395 396 CALL collision_efficiency_rogers( mean_r, & 397 particles(n)%radius, & 398 effective_coll_efficiency ) 399 400 ELSE 401 effective_coll_efficiency = 0.0 402 ENDIF 403 404 IF ( effective_coll_efficiency > 1.0 .OR. & 405 effective_coll_efficiency < 0.0 ) & 406 THEN 407 WRITE( message_string, * ) 'collision_efficien' , & 408 'cy out of range:' ,effective_coll_efficiency 409 CALL message( 'lpm_droplet_collision', 'PA0145', 2, & 410 2, -1, 6, 1 ) 411 ENDIF 412 413 ! 414 !-- Interpolation of ... 415 ii = particles(n)%x * ddx 416 jj = particles(n)%y * ddy 417 kk = ( particles(n)%z + 0.5 * dz ) / dz 418 419 x = particles(n)%x - ii * dx 420 y = particles(n)%y - jj * dy 421 aa = x**2 + y**2 422 bb = ( dx - x )**2 + y**2 423 cc = x**2 + ( dy - y )**2 424 dd = ( dx - x )**2 + ( dy - y )**2 425 gg = aa + bb + cc + dd 426 427 ql_int_l = ( (gg-aa) * ql(kk,jj,ii) + (gg-bb) * & 428 ql(kk,jj,ii+1) & 429 + (gg-cc) * ql(kk,jj+1,ii) + ( gg-dd ) * & 430 ql(kk,jj+1,ii+1) & 431 ) / ( 3.0 * gg ) 432 433 ql_int_u = ( (gg-aa) * ql(kk+1,jj,ii) + (gg-bb) * & 434 ql(kk+1,jj,ii+1) & 435 + (gg-cc) * ql(kk+1,jj+1,ii) + (gg-dd) * & 436 ql(kk+1,jj+1,ii+1) & 437 ) / ( 3.0 * gg ) 438 439 ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz *& 440 ( ql_int_u - ql_int_l ) 441 442 ! 443 !-- Interpolate u velocity-component 444 ii = ( particles(n)%x + 0.5 * dx ) * ddx 445 jj = particles(n)%y * ddy 446 kk = ( particles(n)%z + 0.5 * dz ) / dz ! only if eqist 447 448 IF ( ( particles(n)%z - zu(kk) ) > (0.5*dz) ) kk = kk+1 449 450 x = particles(n)%x + ( 0.5 - ii ) * dx 451 y = particles(n)%y - jj * dy 452 aa = x**2 + y**2 453 bb = ( dx - x )**2 + y**2 454 cc = x**2 + ( dy - y )**2 455 dd = ( dx - x )**2 + ( dy - y )**2 456 gg = aa + bb + cc + dd 457 458 u_int_l = ( (gg-aa) * u(kk,jj,ii) + (gg-bb) * & 459 u(kk,jj,ii+1) & 460 + (gg-cc) * u(kk,jj+1,ii) + (gg-dd) * & 461 u(kk,jj+1,ii+1) & 462 ) / ( 3.0 * gg ) - u_gtrans 463 IF ( kk+1 == nzt+1 ) THEN 464 u_int = u_int_l 465 ELSE 466 u_int_u = ( (gg-aa) * u(kk+1,jj,ii) + (gg-bb) * & 467 u(kk+1,jj,ii+1) & 468 + (gg-cc) * u(kk+1,jj+1,ii) + (gg-dd) * & 469 u(kk+1,jj+1,ii+1) & 470 ) / ( 3.0 * gg ) - u_gtrans 471 u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz & 472 * ( u_int_u - u_int_l ) 473 ENDIF 474 475 ! 476 !-- Same procedure for interpolation of the v velocity-com- 477 !-- ponent (adopt index k from u velocity-component) 478 ii = particles(n)%x * ddx 479 jj = ( particles(n)%y + 0.5 * dy ) * ddy 480 481 x = particles(n)%x - ii * dx 482 y = particles(n)%y + ( 0.5 - jj ) * dy 483 aa = x**2 + y**2 484 bb = ( dx - x )**2 + y**2 485 cc = x**2 + ( dy - y )**2 486 dd = ( dx - x )**2 + ( dy - y )**2 487 gg = aa + bb + cc + dd 488 489 v_int_l = ( ( gg-aa ) * v(kk,jj,ii) + ( gg-bb ) * & 490 v(kk,jj,ii+1) & 491 + ( gg-cc ) * v(kk,jj+1,ii) + ( gg-dd ) * & 492 v(kk,jj+1,ii+1) & 493 ) / ( 3.0 * gg ) - v_gtrans 494 IF ( kk+1 == nzt+1 ) THEN 495 v_int = v_int_l 496 ELSE 497 v_int_u = ( (gg-aa) * v(kk+1,jj,ii) + (gg-bb) * & 498 v(kk+1,jj,ii+1) & 499 + (gg-cc) * v(kk+1,jj+1,ii) + (gg-dd) * & 500 v(kk+1,jj+1,ii+1) & 501 ) / ( 3.0 * gg ) - v_gtrans 502 v_int = v_int_l + ( particles(n)%z - zu(kk) ) / dz & 503 * ( v_int_u - v_int_l ) 504 ENDIF 505 506 ! 507 !-- Same procedure for interpolation of the w velocity-com- 508 !-- ponent (adopt index i from v velocity-component) 509 jj = particles(n)%y * ddy 510 kk = particles(n)%z / dz 511 512 x = particles(n)%x - ii * dx 513 y = particles(n)%y - jj * dy 514 aa = x**2 + y**2 515 bb = ( dx - x )**2 + y**2 516 cc = x**2 + ( dy - y )**2 517 dd = ( dx - x )**2 + ( dy - y )**2 518 gg = aa + bb + cc + dd 519 520 w_int_l = ( ( gg-aa ) * w(kk,jj,ii) + ( gg-bb ) * & 521 w(kk,jj,ii+1) & 522 + ( gg-cc ) * w(kk,jj+1,ii) + ( gg-dd ) * & 523 w(kk,jj+1,ii+1) & 524 ) / ( 3.0 * gg ) 525 IF ( kk+1 == nzt+1 ) THEN 526 w_int = w_int_l 527 ELSE 528 w_int_u = ( (gg-aa) * w(kk+1,jj,ii) + (gg-bb) * & 529 w(kk+1,jj,ii+1) & 530 + (gg-cc) * w(kk+1,jj+1,ii) + (gg-dd) * & 531 w(kk+1,jj+1,ii+1) & 532 ) / ( 3.0 * gg ) 533 w_int = w_int_l + ( particles(n)%z - zw(kk) ) / dz & 534 * ( w_int_u - w_int_l ) 535 ENDIF 536 537 ! 538 !-- Change in radius due to collision 539 delta_r = effective_coll_efficiency / 3.0_wp & 540 * pi * sl_r3 * ddx * ddy / dz & 541 * SQRT( ( u_int - particles(n)%speed_x )**2 & 542 + ( v_int - particles(n)%speed_y )**2 & 543 + ( w_int - particles(n)%speed_z )**2 & 544 ) * dt_3d 545 ! 546 !-- Change in volume due to collision 547 delta_v = particles(n)%weight_factor & 548 * ( ( particles(n)%radius + delta_r )**3 & 549 - particles(n)%radius**3 ) 550 551 ! 552 !-- Check if collected particles provide enough LWC for 553 !-- volume change of collector particle 554 IF ( delta_v >= sl_r3 .AND. sl_r3 > 0.0 ) THEN 555 556 delta_r = ( ( sl_r3/particles(n)%weight_factor ) & 557 + particles(n)%radius**3 )**( 1.0_wp/3.0_wp ) & 558 - particles(n)%radius 559 560 DO is = n-1, psi, -1 561 IF ( particles(is)%radius < & 562 particles(n)%radius ) THEN 563 particles(is)%weight_factor = 0.0 564 particle_mask(is) = .FALSE. 565 deleted_particles = deleted_particles + 1 566 ENDIF 567 ENDDO 568 569 ELSE IF ( delta_v < sl_r3 .AND. sl_r3 > 0.0 ) THEN 570 571 DO is = n-1, psi, -1 572 IF ( particles(is)%radius < particles(n)%radius & 573 .AND. sl_r3 > 0.0 ) THEN 574 particles(is)%weight_factor = & 575 ( ( particles(is)%weight_factor & 576 * ( particles(is)%radius**3 ) ) & 577 - ( delta_v & 578 * particles(is)%weight_factor & 579 * ( particles(is)%radius**3 ) & 580 / sl_r3 ) ) & 581 / ( particles(is)%radius**3 ) 582 583 IF ( particles(is)%weight_factor < 0.0 ) THEN 584 WRITE( message_string, * ) 'negative ', & 585 'weighting factor: ', & 586 particles(is)%weight_factor 587 CALL message( 'lpm_droplet_collision', & 588 'PA0039', & 589 2, 2, -1, 6, 1 ) 590 ENDIF 591 ENDIF 592 ENDDO 593 594 ENDIF 595 596 particles(n)%radius = particles(n)%radius + delta_r 597 ql_vp(k,j,i) = ql_vp(k,j,i) + & 598 particles(n)%weight_factor * & 599 ( particles(n)%radius**3 ) 600 ENDDO 601 602 ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor & 603 * particles(psi)%radius**3 604 605 ENDIF ! collision kernel 606 607 ELSE IF ( prt_count(k,j,i) == 1 ) THEN 608 609 psi = prt_start_index(k,j,i) 610 611 ! 612 !-- Calculate change of weighting factor due to self collision 613 IF ( ( hall_kernel .OR. wang_kernel ) .AND. & 614 use_kernel_tables ) THEN 615 616 IF ( wang_kernel ) THEN 617 eclass = INT( diss(k,j,i) * 1.0E4_wp / 1000.0_wp * & 618 dissipation_classes ) + 1 619 epsilon = diss(k,j,i) 620 ELSE 621 epsilon = 0.0 622 ENDIF 623 IF ( hall_kernel .OR. epsilon * 1.0E4_wp < 0.001 ) THEN 624 eclass = 0 ! Hall kernel is used 625 ELSE 626 eclass = MIN( dissipation_classes, eclass ) 627 ENDIF 628 629 ddV = ddx * ddy / dz 630 rclass_l = particles(psi)%class 631 sum3 = 1 - dt_3d * ddV * & 632 ( ckernel(rclass_l,rclass_l,eclass) * & 633 ( particles(psi)%weight_factor-1 ) * 0.5 ) 634 635 particles(psi)%radius = ( particles(psi)%radius**3 / & 636 sum3 )**0.33333333333333_wp 637 particles(psi)%weight_factor = particles(psi)%weight_factor & 638 * sum3 639 640 ELSE IF ( ( hall_kernel .OR. wang_kernel ) .AND. & 641 .NOT. use_kernel_tables ) THEN 642 ! 643 !-- Collision efficiencies are calculated for every new 644 !-- grid box. First, allocate memory for kernel table. 645 !-- Third dimension is 1, because table is re-calculated for 646 !-- every new dissipation value. 647 ALLOCATE( ckernel(psi:psi, psi:psi, 1:1) ) 648 ! 649 !-- Now calculate collision efficiencies for this box 650 CALL recalculate_kernel( i, j, k ) 651 652 ddV = ddx * ddy / dz 653 sum3 = 1 - dt_3d * ddV * ( ckernel(psi,psi,1) * & 654 ( particles(psi)%weight_factor - 1 ) * 0.5 ) 655 656 particles(psi)%radius = ( particles(psi)%radius**3 / & 657 sum3 )**0.33333333333333_wp 658 particles(psi)%weight_factor = particles(psi)%weight_factor & 659 * sum3 660 661 DEALLOCATE( ckernel ) 662 ENDIF 663 664 ql_vp(k,j,i) = particles(psi)%weight_factor * & 665 particles(psi)%radius**3 597 ENDDO 598 666 599 ENDIF 667 600 668 ! 669 !-- Check if condensation of LWC was conserved during collision 670 !-- process 671 IF ( ql_v(k,j,i) /= 0.0 ) THEN 672 IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001 .OR. & 673 ql_vp(k,j,i) / ql_v(k,j,i) <= 0.9999 ) THEN 674 WRITE( message_string, * ) 'LWC is not conserved during',& 675 ' collision! ', & 676 'LWC after condensation: ', & 677 ql_v(k,j,i), & 678 ' LWC after collision: ', & 679 ql_vp(k,j,i) 680 CALL message( 'lpm_droplet_collision', 'PA0040', & 681 2, 2, -1, 6, 1 ) 682 ENDIF 683 ENDIF 684 601 particles(n)%radius = particles(n)%radius + delta_r 602 ql_vp(k,j,i) = ql_vp(k,j,i) + & 603 particles(n)%weight_factor * & 604 ( particles(n)%radius**3 ) 685 605 ENDDO 686 ENDDO 687 ENDDO 606 607 ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor & 608 * particles(psi)%radius**3 609 610 ENDIF ! collision kernel 611 612 ELSE IF ( prt_count(k,j,i) == 1 ) THEN 613 614 psi = 1 615 616 ! 617 !-- Calculate change of weighting factor due to self collision 618 IF ( ( hall_kernel .OR. wang_kernel ) .AND. & 619 use_kernel_tables ) THEN 620 621 IF ( wang_kernel ) THEN 622 eclass = INT( diss(k,j,i) * 1.0E4_wp / 1000.0_wp * & 623 dissipation_classes ) + 1 624 epsilon = diss(k,j,i) 625 ELSE 626 epsilon = 0.0_wp 627 ENDIF 628 IF ( hall_kernel .OR. epsilon * 1.0E4_wp < 0.001_wp ) THEN 629 eclass = 0 ! Hall kernel is used 630 ELSE 631 eclass = MIN( dissipation_classes, eclass ) 632 ENDIF 633 634 ddV = ddx * ddy / dz 635 rclass_l = particles(psi)%class 636 sum3 = 1 - dt_3d * ddV * & 637 ( ckernel(rclass_l,rclass_l,eclass) * & 638 ( particles(psi)%weight_factor-1 ) * 0.5_wp ) 639 640 particles(psi)%radius = ( particles(psi)%radius**3 / & 641 sum3 )**0.33333333333333_wp 642 particles(psi)%weight_factor = particles(psi)%weight_factor & 643 * sum3 644 645 ELSE IF ( ( hall_kernel .OR. wang_kernel ) .AND. & 646 .NOT. use_kernel_tables ) THEN 647 ! 648 !-- Collision efficiencies are calculated for every new 649 !-- grid box. First, allocate memory for kernel table. 650 !-- Third dimension is 1, because table is re-calculated for 651 !-- every new dissipation value. 652 ALLOCATE( ckernel(psi:psi, psi:psi, 1:1) ) 653 ! 654 !-- Now calculate collision efficiencies for this box 655 CALL recalculate_kernel( i, j, k ) 656 657 ddV = ddx * ddy / dz 658 sum3 = 1 - dt_3d * ddV * ( ckernel(psi,psi,1) * & 659 ( particles(psi)%weight_factor - 1 ) * 0.5_wp ) 660 661 particles(psi)%radius = ( particles(psi)%radius**3 / & 662 sum3 )**0.33333333333333_wp 663 particles(psi)%weight_factor = particles(psi)%weight_factor & 664 * sum3 665 666 DEALLOCATE( ckernel ) 667 ENDIF 668 669 ql_vp(k,j,i) = particles(psi)%weight_factor * & 670 particles(psi)%radius**3 671 ENDIF 672 673 ! 674 !-- Check if condensation of LWC was conserved during collision process 675 IF ( ql_v(k,j,i) /= 0.0_wp ) THEN 676 IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001_wp .OR. & 677 ql_vp(k,j,i) / ql_v(k,j,i) <= 0.9999_wp ) THEN 678 WRITE( message_string, * ) 'LWC is not conserved during',& 679 ' collision! ', & 680 'LWC after condensation: ', & 681 ql_v(k,j,i), & 682 ' LWC after collision: ', & 683 ql_vp(k,j,i) 684 CALL message( 'lpm_droplet_collision', 'PA0040', & 685 2, 2, -1, 6, 1 ) 686 ENDIF 687 ENDIF 688 688 689 689 CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'stop' ) 690 690 691 692 691 END SUBROUTINE lpm_droplet_collision
Note: See TracChangeset
for help on using the changeset viewer.
