Changeset 799 for palm/trunk/SOURCE/wang_kernel.f90
- Timestamp:
- Dec 21, 2011 5:48:03 PM (13 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/wang_kernel.f90 (modified) (25 diffs)
Legend:
- Unmodified
- Added
- Removed
-
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.
