Changeset 799 for palm/trunk/SOURCE/advec_particles.f90
- Timestamp:
- Dec 21, 2011 5:48:03 PM (12 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/advec_particles.f90 (modified) (17 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_particles.f90
r792 r799 4 4 ! Current revisions: 5 5 ! ------------------ 6 ! particle arrays (particles, particles_temp) implemented as pointers in7 ! order to speed up sorting (see routine sort_particles)6 ! Implementation of Wang collision kernel and corresponding new parameter 7 ! wang_collision_kernel 8 8 ! 9 9 ! Former revisions: … … 11 11 ! $Id$ 12 12 ! 13 ! 792 2011-12-01 raasch 14 ! particle arrays (particles, particles_temp) implemented as pointers in 15 ! order to speed up sorting (see routine sort_particles) 16 ! 13 17 ! 759 2011-09-15 13:58:31Z raasch 14 18 ! Splitting of parallel I/O (routine write_particles) … … 17 21 ! Declaration of de_dx, de_dy, de_dz adapted to additional 18 22 ! ghost points. Furthermore the calls of exchange_horiz were modified. 19 !20 23 ! 21 24 ! 622 2010-12-10 08:08:13Z raasch … … 127 130 USE random_function_mod 128 131 USE statistics 132 USE wang_kernel_mod 129 133 130 134 IMPLICIT NONE 131 135 132 136 INTEGER :: agp, deleted_particles, deleted_tails, i, ie, ii, inc, is, j, & 133 jj, js, k, kk, kw, m, n, nc, n n, num_gp, psi, tlength, &134 t rlp_count, trlp_count_sum, trlp_count_recv,&137 jj, js, k, kk, kw, m, n, nc, nd, nn, num_gp, pse, psi, & 138 tlength, trlp_count, trlp_count_sum, trlp_count_recv, & 135 139 trlp_count_recv_sum, trlpt_count, trlpt_count_recv, & 136 140 trnp_count, trnp_count_sum, trnp_count_recv, & … … 139 143 trrp_count_recv_sum, trrpt_count, trrpt_count_recv, & 140 144 trsp_count, trsp_count_sum, trsp_count_recv, & 141 trsp_count_recv_sum, trspt_count, trspt_count_recv , nd145 trsp_count_recv_sum, trspt_count, trspt_count_recv 142 146 143 147 INTEGER :: gp_outside_of_building(1:8) … … 151 155 dt_particle, dt_particle_m, d_radius, d_sum, e_a, e_int, & 152 156 e_int_l, e_int_u, e_mean_int, e_s, exp_arg, exp_term, fs_int, & 153 gg, lagr_timescale, mean_r, new_r, p_int, pt_int, pt_int_l,&154 pt_int _u, q_int, q_int_l, q_int_u, ql_int, ql_int_l, ql_int_u,&155 random_gauss, sl_r3, sl_r4, t_int, u_int, u_int_l, u_int_u,&156 vv_int, v_int, v_int_l, v_int_u, w_int, w_int_l, w_int_u,&157 x, y157 gg, integral, lagr_timescale, lw_max, mean_r, new_r, p_int, & 158 pt_int, pt_int_l, pt_int_u, q_int, q_int_l, q_int_u, ql_int, & 159 ql_int_l, ql_int_u, random_gauss, sl_r3, sl_r4, t_int, u_int, & 160 u_int_l, u_int_u,vv_int, v_int, v_int_l, v_int_u, w_int, & 161 w_int_l, w_int_u, x, y 158 162 159 163 REAL, DIMENSION(1:30) :: de_dxi, de_dyi, de_dzi, dissi, d_gp_pl, ei 164 REAL, DIMENSION(:,:), ALLOCATABLE :: kern 160 165 161 166 REAL :: location(1:30,1:3) … … 399 404 CALL cpu_log( log_point_s(43), 'advec_part_coll', 'start' ) 400 405 406 IF ( wang_collision_kernel ) THEN 407 408 DO i = nxl, nxr 409 DO j = nys, nyn 410 DO k = nzb+1, nzt 411 ! 412 !-- Collision requires at least two particles in the box 413 IF ( prt_count(k,j,i) > 1 ) THEN 414 ! 415 !-- First, sort particles within the gridbox by their size, 416 !-- using Shell's method (see Numerical Recipes) 417 !-- NOTE: In case of using particle tails, the re-sorting of 418 !-- ---- tails would have to be included here! 419 psi = prt_start_index(k,j,i) - 1 420 inc = 1 421 DO WHILE ( inc <= prt_count(k,j,i) ) 422 inc = 3 * inc + 1 423 ENDDO 424 425 DO WHILE ( inc > 1 ) 426 inc = inc / 3 427 DO is = inc+1, prt_count(k,j,i) 428 tmp_particle = particles(psi+is) 429 js = is 430 DO WHILE ( particles(psi+js-inc)%radius > & 431 tmp_particle%radius ) 432 particles(psi+js) = particles(psi+js-inc) 433 js = js - inc 434 IF ( js <= inc ) EXIT 435 ENDDO 436 particles(psi+js) = tmp_particle 437 ENDDO 438 ENDDO 439 440 psi = prt_start_index(k,j,i) 441 pse = psi + prt_count(k,j,i)-1 442 443 ALLOCATE( kern(psi:pse,psi:pse) ) 444 445 ! 446 !-- Calculate collision kernel for all particles in grid box 447 CALL colker( i, j, k, kern ) 448 ! 449 !-- collison kernel is calculated in cm**3/s but needed in m**3/s 450 kern = kern * 1.0E-6 451 452 DO n = pse, psi+1, -1 453 454 integral = 0.0 455 lw_max = 0.0 456 457 ! 458 !-- Calculate growth of collector particle radius using all 459 !-- droplets smaller than current droplet 460 DO is = psi, n-1 461 462 integral = integral + & 463 ( particles(is)%radius**3 * kern(n,is) * & 464 particles(is)%weight_factor ) 465 ! 466 !-- Calculate volume of liquid water of the collected 467 !-- droplets which is the maximum liquid water available 468 !-- for droplet growth 469 lw_max = lw_max + ( particles(is)%radius**3 * & 470 particles(is)%weight_factor ) 471 472 ENDDO 473 474 ! 475 !-- Change in radius of collector droplet due to collision 476 delta_r = 1.0 / ( 3.0 * particles(n)%radius**2 ) * & 477 integral * dt_3d * ddx * ddy / dz 478 479 ! 480 !-- Change in volume of collector droplet due to collision 481 delta_v = particles(n)%weight_factor & 482 * ( ( particles(n)%radius + delta_r )**3 & 483 - particles(n)%radius**3 ) 484 485 IF ( lw_max < delta_v .AND. delta_v > 0.0 ) THEN 486 487 PRINT*, 'Particle has grown to much because the & 488 change of volume of particle is larger than & 489 liquid water available!' 490 491 ELSEIF ( lw_max == delta_v .AND. delta_v > 0.0 ) THEN 492 493 DO is = psi, n-1 494 495 particles(is)%weight_factor = 0.0 496 particle_mask(is) = .FALSE. 497 deleted_particles = deleted_particles + 1 498 499 ENDDO 500 501 ELSEIF ( lw_max > delta_v .AND. delta_v > 0.0 ) THEN 502 ! 503 !-- Calculate new weighting factor of collected droplets 504 DO is = psi, n-1 505 506 particles(is)%weight_factor = & 507 particles(is)%weight_factor - & 508 ( ( kern(n,is) * particles(is)%weight_factor / & 509 integral ) * delta_v ) 510 511 IF ( particles(is)%weight_factor < 0.0 ) THEN 512 WRITE( message_string, * ) 'negative ', & 513 'weighting factor: ', & 514 particles(is)%weight_factor 515 CALL message( 'advec_particles', '', 2, 2, & 516 -1, 6, 1 ) 517 518 ELSEIF ( particles(is)%weight_factor == 0.0 ) THEN 519 520 particles(is)%weight_factor = 0.0 521 particle_mask(is) = .FALSE. 522 deleted_particles = deleted_particles + 1 523 524 ENDIF 525 526 ENDDO 527 528 ENDIF 529 530 particles(n)%radius = particles(n)%radius + delta_r 531 ql_vp(k,j,i) = ql_vp(k,j,i) + particles(n)%weight_factor & 532 * particles(n)%radius**3 533 534 ENDDO 535 536 ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor & 537 * particles(psi)%radius**3 538 539 DEALLOCATE( kern ) 540 541 ELSE IF ( prt_count(k,j,i) == 1 ) THEN 542 543 psi = prt_start_index(k,j,i) 544 ql_vp(k,j,i) = particles(psi)%weight_factor * & 545 particles(psi)%radius**3 546 ENDIF 547 548 ! 549 !-- Check if condensation of LWC was conserved during collision 550 !-- process 551 IF ( ql_v(k,j,i) /= 0.0 ) THEN 552 IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001 .OR. & 553 ql_vp(k,j,i) / ql_v(k,j,i) <= 0.9999 ) THEN 554 WRITE( message_string, * ) 'LWC is not conserved during',& 555 ' collision! ', & 556 'LWC after condensation: ', & 557 ql_v(k,j,i), & 558 ' LWC after collision: ', & 559 ql_vp(k,j,i) 560 CALL message( 'advec_particles', '', 2, 2, -1, 6, & 561 1 ) 562 ENDIF 563 ENDIF 564 565 ENDDO 566 ENDDO 567 ENDDO 568 569 ELSE 570 401 571 DO i = nxl, nxr 402 572 DO j = nys, nyn … … 683 853 ENDDO 684 854 ENDDO 855 856 ENDIF 685 857 686 858 CALL cpu_log( log_point_s(43), 'advec_part_coll', 'stop' ) … … 1200 1372 de_dzi(num_gp) = de_dz(k+1,j+1,i) 1201 1373 ENDIF 1202 1374 1203 1375 IF ( k > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) & 1204 1376 THEN … … 1620 1792 de_dzi(num_gp) = 0.0 1621 1793 ENDIF 1622 1794 1623 1795 ! 1624 1796 !-- If wall between gridpoint 5 and gridpoint 7, then … … 1971 2143 #else 1972 2144 dt_3d_reached = dt_3d_reached_l 1973 #endif 2145 #endif 1974 2146 1975 2147 CALL cpu_log( log_point_s(44), 'advec_part_advec', 'stop' ) … … 2162 2334 !-- data 2163 2335 DO n = 1, number_of_particles 2164 2336 i = ( particles(n)%x + 0.5 * dx ) * ddx 2165 2337 ! 2166 2338 !-- Above calculation does not work for indices less than zero 2167 2339 IF ( particles(n)%x < -0.5 * dx ) i = -1 2168 2340 2169 2341 IF ( i < nxl ) THEN … … 2225 2397 nn = particles(n)%tail_id 2226 2398 2227 2399 i = ( particles(n)%x + 0.5 * dx ) * ddx 2228 2400 ! 2229 2401 !-- Above calculation does not work for indices less than zero 2230 2402 IF ( particles(n)%x < - 0.5 * dx ) i = -1 2231 2403 2232 2404 IF ( i < nxl ) THEN … … 2602 2774 DO n = 1, number_of_particles 2603 2775 IF ( particle_mask(n) ) THEN 2604 2776 j = ( particles(n)%y + 0.5 * dy ) * ddy 2605 2777 ! 2606 2778 !-- Above calculation does not work for indices less than zero 2607 2779 IF ( particles(n)%y < -0.5 * dy ) j = -1 2608 2780 2609 2781 IF ( j < nys ) THEN … … 2669 2841 !-- moved. 2670 2842 IF ( particle_mask(n) ) THEN 2671 2843 j = ( particles(n)%y + 0.5 * dy ) * ddy 2672 2844 ! 2673 2845 !-- Above calculation does not work for indices less than zero 2674 2846 IF ( particles(n)%y < -0.5 * dy ) j = -1 2675 2847 2676 2848 IF ( j < nys ) THEN … … 3815 3987 3816 3988 CHARACTER (LEN=10) :: particle_binary_version 3817 3818 3989 INTEGER :: i 3819 3990 … … 3879 4050 3880 4051 END SUBROUTINE write_particles 4052 3881 4053 3882 4054
Note: See TracChangeset
for help on using the changeset viewer.
