Changeset 799 for palm/trunk/SOURCE
- Timestamp:
- Dec 21, 2011 5:48:03 PM (13 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 4 edited
-
advec_particles.f90 (modified) (17 diffs)
-
compute_vpt.f90 (modified) (2 diffs)
-
interaction_droplets_ptq.f90 (modified) (3 diffs)
-
wang_kernel.f90 (modified) (25 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 -
palm/trunk/SOURCE/compute_vpt.f90
r484 r799 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! Bugfix: ql is now included in calculation of vpt in case of 7 ! cloud droplets 7 8 ! 8 9 ! Former revisions: … … 32 33 INTEGER :: k 33 34 34 IF ( .NOT. cloud_physics ) THEN35 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 35 36 vpt = pt * ( 1.0 + 0.61 * q ) 36 ELSE 37 ELSE IF (cloud_physics) THEN 37 38 DO k = nzb, nzt+1 38 39 vpt(k,:,:) = ( pt(k,:,:) + pt_d_t(k) * l_d_cp * ql(k,:,:) ) * & 39 40 ( 1.0 + 0.61 * q(k,:,:) - 1.61 * ql(k,:,:) ) 40 41 ENDDO 42 ELSE 43 vpt = pt * ( 1.0 + 0.61 * q - 1.61 * ql ) 41 44 ENDIF 42 45 -
palm/trunk/SOURCE/interaction_droplets_ptq.f90
r484 r799 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! Bugfix: pt_d_t(k) was missing in calculation of pt_p 7 7 ! 8 8 ! Former revisions: … … 53 53 DO k = nzb_2d(j,i)+1, nzt 54 54 q_p(k,j,i) = q_p(k,j,i) - ql_c(k,j,i) 55 pt_p(k,j,i) = pt_p(k,j,i) + l_d_cp * ql_c(k,j,i) 55 pt_p(k,j,i) = pt_p(k,j,i) + l_d_cp * ql_c(k,j,i) * pt_d_t(k) 56 56 ENDDO 57 57 ENDDO … … 80 80 DO k = nzb_2d(j,i)+1, nzt 81 81 q_p(k,j,i) = q_p(k,j,i) - ql_c(k,j,i) 82 pt_p(k,j,i) = pt_p(k,j,i) + l_d_cp * ql_c(k,j,i) 82 pt_p(k,j,i) = pt_p(k,j,i) + l_d_cp * ql_c(k,j,i) * pt_d_t(k) 83 83 ENDDO 84 84 -
palm/trunk/SOURCE/wang_kernel.f90
r792 r799 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! speed optimizations + formatting 6 ! speed optimizations and formatting 7 ! Bugfix: iq=1 is not allowed (routine effic) 8 ! Bugfix: replaced stop by ec=0.0 in case of very small ec (routine effic) 7 9 ! 8 10 ! Former revisions: … … 169 171 ! of droplets, air density, air viscosity, turbulent dissipation rate, 170 172 ! taylor microscale reynolds number, gravitational acceleration 173 171 174 USE constants 172 175 USE cloud_parameters … … 176 179 IMPLICIT NONE 177 180 178 REAL :: airvisc, airdens, waterdens, gravity, anu,Relamda, &181 REAL :: Relamda, & 179 182 Tl, Lf, tauk, eta, vk, ao, lambda, tt, z, be, & 180 183 bbb, d1, e1, d2, e2, ccc, b1, c1, b2, c2, v1xysq, & … … 182 185 wrFIN, SSt, XX, YY, c1_gr, ao_gr, fao_gr, rc, grFIN, v1, t1, v2, t2, rrp 183 186 184 REAL, DIMENSION(pstart:pend) :: vST, tau, St 187 REAL, SAVE :: airvisc, airdens, anu, gravity, waterdens 188 189 REAL, DIMENSION(pstart:pend) :: St, tau 190 191 LOGICAL, SAVE :: first = .TRUE. 185 192 186 193 INTEGER :: i, j 187 194 188 airvisc = 0.1818 !dynamic viscosity in mg/cm*s 189 airdens = 1.2250 !air density in mg/cm**3 190 waterdens = 1000.0 !water density in mg/cm**3 191 gravity = 980.6650 !in cm/s**2 192 193 anu = airvisc/airdens ! kinetic viscosity in cm**2/s 194 195 Relamda = urms**2.0*sqrt(15.0/epsilon/anu) 196 197 Tl = urms**2.0/epsilon !in s 198 Lf = 0.5 * (urms**3.0)/epsilon !in cm 195 ! 196 !-- Initial assignment of constants 197 IF ( first ) THEN 198 199 first = .FALSE. 200 airvisc = 0.1818 !dynamic viscosity in mg/cm*s 201 airdens = 1.2250 !air density in mg/cm**3 202 waterdens = 1000.0 !water density in mg/cm**3 203 gravity = 980.6650 !in cm/s**2 204 anu = airvisc/airdens ! kinetic viscosity in cm**2/s 205 206 ENDIF 207 208 Relamda = urms**2*sqrt(15.0/epsilon/anu) 209 210 Tl = urms**2/epsilon !in s 211 Lf = 0.5 * (urms**3)/epsilon !in cm 199 212 200 213 tauk = (anu/epsilon)**0.5 !in s 201 eta = (anu**3 .0/epsilon)**0.25 !in cm214 eta = (anu**3/epsilon)**0.25 !in cm 202 215 vk = eta/tauk !in cm/s 203 216 … … 211 224 212 225 DO i = pstart, pend 213 vST(i) = winf(i) !in cm/s 214 tau(i) = vST(i)/gravity !in s 226 tau(i) = winf(i)/gravity !in s 215 227 St(i) = tau(i)/tauk 216 228 ENDDO … … 222 234 be = sqrt(2.0)*lambda/Lf 223 235 224 bbb = sqrt(1.0-2.0*be**2 .0)236 bbb = sqrt(1.0-2.0*be**2) 225 237 d1 = (1.+bbb)/2.0/bbb 226 238 e1 = Lf*(1.0+bbb)/2.0 !in cm … … 228 240 e2 = Lf*(1.0-bbb)/2.0 !in cm 229 241 230 ccc = sqrt(1.0-2.0*z**2 .0)242 ccc = sqrt(1.0-2.0*z**2) 231 243 b1 = (1.+ccc)/2./ccc 232 244 c1 = Tl*(1.+ccc)/2. !in s … … 235 247 236 248 DO i = pstart, pend 237 v1 = vST(i) !in cm/s249 v1 = winf(i) !in cm/s 238 250 t1 = tau(i) !in s 239 251 240 252 DO j = pstart,i 241 253 rrp = (particles(i)%radius + particles(j)%radius) * 100.0 !radius in cm 242 v2 = vST(j) !in cm/s254 v2 = winf(j) !in cm/s 243 255 t2 = tau(j) !in s 244 256 … … 247 259 - b2*d1*PHI(c2,e1,v1,t1) & 248 260 + b2*d2*PHI(c2,e2,v1,t1) 249 v1xysq = v1xysq * urms**2 .0/t1 !in cm**2/s**2261 v1xysq = v1xysq * urms**2/t1 !in cm**2/s**2 250 262 251 263 vrms1xy= sqrt(v1xysq) !in cm/s … … 255 267 - b2*d1*PHI(c2,e1,v2,t2) & 256 268 + b2*d2*PHI(c2,e2,v2,t2) 257 v2xysq = v2xysq * urms**2 .0/t2 !in cm**2/s**2269 v2xysq = v2xysq * urms**2/t2 !in cm**2/s**2 258 270 259 271 vrms2xy= sqrt(v2xysq) !in cm/s 260 272 261 IF( vST(i).ge.vST(j)) THEN262 v1 = vST(i)273 IF(winf(i).ge.winf(j)) THEN 274 v1 = winf(i) 263 275 t1 = tau(i) 264 v2 = vST(j)276 v2 = winf(j) 265 277 t2 = tau(j) 266 278 ELSE 267 v1 = vST(j)279 v1 = winf(j) 268 280 t1 = tau(j) 269 v2 = vST(i)281 v2 = winf(i) 270 282 t2 = tau(i) 271 283 ENDIF … … 276 288 + b2*d2*ZHI(c2,e2,v1,t1,v2,t2) 277 289 fR = d1 * exp(-rrp/e1) - d2 * exp(-rrp/e2) 278 v1v2xy = v1v2xy * fR * urms**2 .0/tau(i)/tau(j) !in cm**2/s**2279 280 wrtur2xy=vrms1xy**2 .0 + vrms2xy**2.0- 2.*v1v2xy !in cm**2/s**2290 v1v2xy = v1v2xy * fR * urms**2/tau(i)/tau(j) !in cm**2/s**2 291 292 wrtur2xy=vrms1xy**2 + vrms2xy**2 - 2.*v1v2xy !in cm**2/s**2 281 293 282 294 IF (wrtur2xy.le.0.0) wrtur2xy=0.0 283 295 284 wrgrav2=pi/8.*( vST(j)-vST(i))**2.0296 wrgrav2=pi/8.*(winf(j)-winf(i))**2 285 297 286 298 wrFIN = sqrt((2.0/pi)*(wrtur2xy+wrgrav2)) !in cm/s … … 295 307 ENDIF 296 308 297 XX = -0.1988*SSt**4 .0 + 1.5275*SSt**3.0&298 -4.2942*SSt**2 .0+ 5.3406*SSt309 XX = -0.1988*SSt**4 + 1.5275*SSt**3 & 310 -4.2942*SSt**2 + 5.3406*SSt 299 311 300 312 IF(XX.le.0.0) XX = 0.0 … … 304 316 c1_gr = XX/(gravity/(vk/tauk))**YY 305 317 306 ao_gr = ao + (pi/8.)*(gravity/(vk/tauk))**2 .0318 ao_gr = ao + (pi/8.)*(gravity/(vk/tauk))**2 307 319 fao_gr = 20.115 * (ao_gr/Relamda)**0.5 308 320 rc = sqrt( fao_gr * abs(St(j)-St(i)) ) * eta !in cm 309 321 310 grFIN = ((eta**2 .0+rc**2.0)/(rrp**2.0+rc**2.0))**(c1_gr/2.)322 grFIN = ((eta**2+rc**2)/(rrp**2+rc**2))**(c1_gr/2.) 311 323 IF (grFIN.lt.1.0) grFIN = 1.0 312 324 313 gck(i,j) = 2. * pi * rrp**2 .0* wrFIN * grFIN ! in cm**3/s325 gck(i,j) = 2. * pi * rrp**2 * wrFIN * grFIN ! in cm**3/s 314 326 gck(j,i) = gck(i,j) 315 327 … … 330 342 aa1 = 1./tau0 + 1./a + vsett/b 331 343 332 PHI = 1./aa1 - vsett/2.0/b/aa1**2 .0!in s344 PHI = 1./aa1 - vsett/2.0/b/aa1**2 !in s 333 345 334 346 END FUNCTION PHI … … 346 358 aa2 = vsett1/b + 1./tau1 + 1./a 347 359 aa3 = (vsett1-vsett2)/b + 1./tau1 + 1./tau2 348 aa4 = (vsett2/b)**2 .0 - (1./tau2 + 1./a)**2.0360 aa4 = (vsett2/b)**2 - (1./tau2 + 1./a)**2 349 361 aa5 = vsett2/b + 1./tau2 + 1./a 350 362 aa6 = 1./tau1 - 1./a + (1./tau2 + 1./a) * vsett1/vsett2 351 ZHI = (1./aa1 - 1./aa2) * (vsett1-vsett2)/2./b/aa3**2 .0&352 + (4./aa4 - 1./aa5**2 .0 - 1./aa1**2.0) * vsett2/2./b/aa6 &353 + (2.*b/aa2 - 2.*b/aa1 - vsett1/aa2**2 .0 + vsett2/aa1**2.0) &363 ZHI = (1./aa1 - 1./aa2) * (vsett1-vsett2)/2./b/aa3**2 & 364 + (4./aa4 - 1./aa5**2 - 1./aa1**2) * vsett2/2./b/aa6 & 365 + (2.*b/aa2 - 2.*b/aa1 - vsett1/aa2**2 + vsett2/aa1**2) & 354 366 * 1./2./b/aa3 ! in s**2 367 355 368 END FUNCTION ZHI 356 369 357 358 370 !------------------------------------------------------------------------------! 359 371 ! SUBROUTINE for calculation of terminal velocity winf 360 372 !------------------------------------------------------------------------------! 361 362 373 SUBROUTINE fallg 363 374 … … 371 382 INTEGER :: i, j 372 383 373 REAL :: eta, xlamb, rhoa, rhow, grav, cunh, t0, sigma, stok, stb, phy, py, & 374 x, y, xrey, bond 375 376 REAL, DIMENSION(1:7) :: b 377 REAL, DIMENSION(1:6) :: c 378 REAL, DIMENSION(1:20) :: rat 379 REAL, DIMENSION(1:15) :: r0 380 REAL, DIMENSION(1:15,1:20) :: ecoll 381 382 b = (/-0.318657e1,0.992696e0,-0.153193e-2,-0.987059e-3, & 384 LOGICAL, SAVE :: first = .TRUE. 385 386 REAL, SAVE :: eta, xlamb, rhoa, rhow, grav, cunh, t0, sigma, stok, stb, phy, py 387 388 REAL :: bond, x, xrey, y 389 390 REAL, DIMENSION(1:7), SAVE :: b 391 REAL, DIMENSION(1:6), SAVE :: c 392 393 ! 394 !-- Initial assignment of constants 395 IF ( first ) THEN 396 397 first = .FALSE. 398 b = (/-0.318657e1,0.992696e0,-0.153193e-2,-0.987059e-3, & 383 399 -0.578878e-3,0.855176e-4,-0.327815e-5/) 384 c = (/-0.500015e1,0.523778e1,-0.204914e1,0.475294e0, &400 c = (/-0.500015e1,0.523778e1,-0.204914e1,0.475294e0, & 385 401 -0.542819e-1, 0.238449e-2/) 386 402 387 eta = 1.818e-4 !in poise = g/(cm s) 388 xlamb = 6.62e-6 !in cm 389 390 rhow = 1.0 !in g/cm**3 391 rhoa = 1.225e-3 !in g/cm**3 392 393 grav = 980.665 !in cm/s**2 394 cunh = 1.257 * xlamb !in cm 395 t0 = 273.15 !in K 396 sigma= 76.1 - 0.155 * (293.15 - t0) !in N/m = g/s**2 397 stok = 2.0 * grav * (rhow - rhoa)/(9.0 * eta) ! in 1/(cm s) 398 stb = 32.0 * rhoa * (rhow-rhoa) * grav/(3.0 * eta * eta) ! in 1/cm**3 399 phy = (sigma**3.0) * (rhoa**2.0)/((eta**4.0) * grav * (rhow-rhoa)) 400 py = phy**(1.0/6.0) 403 eta = 1.818e-4 !in poise = g/(cm s) 404 xlamb = 6.62e-6 !in cm 405 406 rhow = 1.0 !in g/cm**3 407 rhoa = 1.225e-3 !in g/cm**3 408 409 grav = 980.665 !in cm/s**2 410 cunh = 1.257 * xlamb !in cm 411 t0 = 273.15 !in K 412 sigma= 76.1 - 0.155 * (293.15 - t0) !in N/m = g/s**2 413 stok = 2.0 * grav * (rhow - rhoa)/(9.0 * eta) ! in 1/(cm s) 414 stb = 32.0 * rhoa * (rhow-rhoa) * grav/(3.0 * eta * eta) ! in 1/cm**3 415 phy = (sigma**3) * (rhoa**2)/((eta**4) * grav * (rhow-rhoa)) 416 py = phy**(1.0/6.0) 417 418 ENDIF 401 419 402 420 !particle radius has to be in cm … … 405 423 IF (particles(j)%radius*100.0 .le. 1.e-3) THEN 406 424 407 winf(j)=stok*((particles(j)%radius*100.0)**2 .+cunh* particles(j)%radius*100.0) !in cm/s425 winf(j)=stok*((particles(j)%radius*100.0)**2+cunh* particles(j)%radius*100.0) !in cm/s 408 426 409 427 ELSEIF (particles(j)%radius*100.0.gt.1.e-3.and.particles(j)%radius*100.0.le.5.35e-2) THEN 410 428 411 x = log(stb*(particles(j)%radius*100.0)**3 .0)429 x = log(stb*(particles(j)%radius*100.0)**3) 412 430 y = 0.0 413 431 … … 421 439 ELSEIF (particles(j)%radius*100.0.gt.5.35e-2) THEN 422 440 423 bond = grav*(rhow-rhoa)*(particles(j)%radius*100.0)**2.0/sigma424 425 441 IF (particles(j)%radius*100.0.gt.0.35) THEN 426 442 bond = grav*(rhow-rhoa) * 0.35 * 0.35/sigma 443 ELSE 444 bond = grav*(rhow-rhoa)*(particles(j)%radius*100.0)**2/sigma 427 445 ENDIF 428 446 … … 435 453 436 454 xrey = py*exp(y) 437 winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0) !in cm/s438 455 439 456 IF (particles(j)%radius*100.0 .gt.0.35) THEN 440 457 winf(j) = xrey * eta/(2.0 * rhoa * 0.35) !in cm/s 458 ELSE 459 winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0) !in cm/s 441 460 ENDIF 442 461 443 462 ENDIF 444 463 ENDDO 445 464 RETURN 446 465 END SUBROUTINE fallg 447 466 … … 457 476 USE particle_attributes 458 477 478 !collision efficiencies of hall kernel 459 479 IMPLICIT NONE 460 480 … … 545 565 DO i = pstart, j 546 566 547 ! DO k = 2, 15 548 ! IF ((particles(j)%radius*1.0E06).le.r0(k).and.(particles(j)%radius*1.0E06).ge.r0(k-1)) THEN 549 ! ir = k 550 ! ELSEIF ((particles(j)%radius*1.0E06).gt.r0(15)) THEN 551 ! ir = 16 552 ! ELSEIF ((particles(j)%radius*1.0E06).lt.r0(1)) THEN 553 ! ir = 1 567 ir = ira(j) 568 569 rq = particles(i)%radius / particles(j)%radius 570 571 ! DO kk = 2, 21 572 ! IF ( rq <= rat(kk) ) THEN 573 ! iq = kk 574 ! EXIT 554 575 ! ENDIF 555 576 ! ENDDO 556 577 557 ir = ira(j)558 559 rq = particles(i)%radius / particles(j)%radius560 561 ! DO kk = 2, 21562 ! IF ( rq .le.rat(kk).and.rq.gt.rat(kk-1)) iq = kk563 ! ENDDO564 565 578 iq = INT( rq * 20 ) + 1 579 iq = MAX(iq , 2) 566 580 567 581 IF ( ir < 16 ) THEN … … 587 601 588 602 ec(i,j) = ec(j,i) 589 IF ( ec(i,j) < 1.0E-20 ) STOP 99603 IF ( ec(i,j) < 1.0E-20 ) ec(i,j) = 0.0 590 604 591 605 ENDDO … … 599 613 ! SUBROUTINE for calculation of enhancement factor collision efficencies 600 614 !------------------------------------------------------------------------------! 601 602 615 SUBROUTINE turb_enhance_eff 603 616 … … 607 620 USE arrays_3d 608 621 609 !collision efficiencies of hall kernel610 622 IMPLICIT NONE 611 623 612 624 INTEGER :: i, ik, ir, iq, j, k, kk 613 REAL :: rq, y1, pp, qq, y2, y3, x1, x2, x3 614 615 REAL, DIMENSION(1:11) :: rat 616 REAL, DIMENSION(1:7) :: r0 617 REAL, DIMENSION(1:7,1:11) :: ecoll_100, ecoll_400 618 619 r0 = (/10., 20., 30.,40., 50., 60.,100./) 620 rat = (/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0/) 625 626 INTEGER, DIMENSION(:), ALLOCATABLE :: ira 627 628 REAL :: rq, y1, particle_radius, pp, qq, y2, y3, x1, x2, x3 629 630 LOGICAL, SAVE :: first = .TRUE. 631 632 REAL, DIMENSION(1:11), SAVE :: rat 633 REAL, DIMENSION(1:7), SAVE :: r0 634 REAL, DIMENSION(1:7,1:11), SAVE :: ecoll_100, ecoll_400 635 636 ! 637 !-- Initial assignment of constants 638 IF ( first ) THEN 639 640 first = .FALSE. 641 642 r0 = (/10., 20., 30.,40., 50., 60.,100./) 643 rat = (/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0/) 621 644 622 645 ! 100 cm^2/s^3 623 ecoll_100(:,1) = (/1.74, 1.74, 1.773, 1.49, 1.207, 1.207, 1.0 /)624 ecoll_100(:,2) = (/1.46, 1.46, 1.421, 1.245, 1.069, 1.069, 1.0 /)625 ecoll_100(:,3) = (/1.32, 1.32, 1.245, 1.123, 1.000, 1.000, 1.0 /)626 ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025, 1.025, 1.0 /)627 ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056, 1.056, 1.0 /)628 ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028, 1.028, 1.0 /)629 ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046, 1.046, 1.0 /)630 ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029, 1.029, 1.0 /)631 ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021, 1.021, 1.0 /)632 ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088, 1.088, 1.0 /)633 ecoll_100(:,11)= (/20.3, 20.3, 14.6 , 8.61, 2.60, 2.60 , 1.0 /)646 ecoll_100(:,1) = (/1.74, 1.74, 1.773, 1.49, 1.207, 1.207, 1.0 /) 647 ecoll_100(:,2) = (/1.46, 1.46, 1.421, 1.245, 1.069, 1.069, 1.0 /) 648 ecoll_100(:,3) = (/1.32, 1.32, 1.245, 1.123, 1.000, 1.000, 1.0 /) 649 ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025, 1.025, 1.0 /) 650 ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056, 1.056, 1.0 /) 651 ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028, 1.028, 1.0 /) 652 ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046, 1.046, 1.0 /) 653 ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029, 1.029, 1.0 /) 654 ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021, 1.021, 1.0 /) 655 ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088, 1.088, 1.0 /) 656 ecoll_100(:,11)= (/20.3, 20.3, 14.6 , 8.61, 2.60, 2.60 , 1.0 /) 634 657 635 658 ! 400 cm^2/s^3 636 ecoll_400(:,1) = (/4.976, 4.976, 3.593, 2.519, 1.445, 1.445, 1.0 /) 637 ecoll_400(:,2) = (/2.984, 2.984, 2.181, 1.691, 1.201, 1.201, 1.0 /) 638 ecoll_400(:,3) = (/1.988, 1.988, 1.475, 1.313, 1.150, 1.150, 1.0 /) 639 ecoll_400(:,4) = (/1.490, 1.490, 1.187, 1.156, 1.126, 1.126, 1.0 /) 640 ecoll_400(:,5) = (/1.249, 1.249, 1.088, 1.090, 1.092, 1.092, 1.0 /) 641 ecoll_400(:,6) = (/1.139, 1.139, 1.130, 1.091, 1.051, 1.051, 1.0 /) 642 ecoll_400(:,7) = (/1.220, 1.220, 1.190, 1.138, 1.086, 1.086, 1.0 /) 643 ecoll_400(:,8) = (/1.325, 1.325, 1.267, 1.165, 1.063, 1.063, 1.0 /) 644 ecoll_400(:,9) = (/1.716, 1.716, 1.345, 1.223, 1.100, 1.100, 1.0 /) 645 ecoll_400(:,10)= (/3.788, 3.788, 1.501, 1.311, 1.120, 1.120, 1.0 /) 646 ecoll_400(:,11)= (/36.52, 36.52, 19.16, 22.80, 26.0, 26.0, 1.0 /) 659 ecoll_400(:,1) = (/4.976, 4.976, 3.593, 2.519, 1.445, 1.445, 1.0 /) 660 ecoll_400(:,2) = (/2.984, 2.984, 2.181, 1.691, 1.201, 1.201, 1.0 /) 661 ecoll_400(:,3) = (/1.988, 1.988, 1.475, 1.313, 1.150, 1.150, 1.0 /) 662 ecoll_400(:,4) = (/1.490, 1.490, 1.187, 1.156, 1.126, 1.126, 1.0 /) 663 ecoll_400(:,5) = (/1.249, 1.249, 1.088, 1.090, 1.092, 1.092, 1.0 /) 664 ecoll_400(:,6) = (/1.139, 1.139, 1.130, 1.091, 1.051, 1.051, 1.0 /) 665 ecoll_400(:,7) = (/1.220, 1.220, 1.190, 1.138, 1.086, 1.086, 1.0 /) 666 ecoll_400(:,8) = (/1.325, 1.325, 1.267, 1.165, 1.063, 1.063, 1.0 /) 667 ecoll_400(:,9) = (/1.716, 1.716, 1.345, 1.223, 1.100, 1.100, 1.0 /) 668 ecoll_400(:,10)= (/3.788, 3.788, 1.501, 1.311, 1.120, 1.120, 1.0 /) 669 ecoll_400(:,11)= (/36.52, 36.52, 19.16, 22.80, 26.0, 26.0, 1.0 /) 670 671 ENDIF 672 673 ! 674 !-- Calculate the radius class index of particles with respect to array r 675 ALLOCATE( ira(pstart:pend) ) 676 677 DO j = pstart, pend 678 particle_radius = particles(j)%radius * 1.0E6 679 DO k = 1, 7 680 IF ( particle_radius < r0(k) ) THEN 681 ira(j) = k 682 EXIT 683 ENDIF 684 ENDDO 685 IF ( particle_radius >= r0(7) ) ira(j) = 8 686 ENDDO 647 687 648 688 ! two-dimensional linear interpolation of the collision efficiency … … 650 690 DO i = pstart, j 651 691 652 DO k = 2, 7 653 IF ((particles(j)%radius*1.0E06).le.r0(k).and.(particles(j)%radius*1.0E06).ge.r0(k-1)) THEN 654 ir = k 655 ELSEIF ((particles(j)%radius*1.0E06).gt.r0(7)) THEN 656 ir = 8 657 ELSEIF ((particles(j)%radius*1.0E06).lt.r0(1)) THEN 658 ir = 1 692 ir = ira(j) 693 694 rq = particles(i)%radius/particles(j)%radius 695 696 DO kk = 2, 11 697 IF ( rq <= rat(kk) ) THEN 698 iq = kk 699 EXIT 659 700 ENDIF 660 ENDDO661 662 rq = particles(i)%radius*1.0E06/(particles(j)%radius*1.0E06)663 664 DO kk = 2, 11665 IF (rq.le.rat(kk).and.rq.gt.rat(kk-1)) iq = kk666 701 ENDDO 667 702
Note: See TracChangeset
for help on using the changeset viewer.
