Changeset 792 for palm/trunk/SOURCE/wang_kernel.f90
- Timestamp:
- Dec 1, 2011 12:23:23 AM (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/wang_kernel.f90
r791 r792 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! speed optimizations + formatting 7 7 ! 8 8 ! Former revisions: … … 29 29 PRIVATE 30 30 31 PUBLIC TurbSD, turb_enhance_eff, effic, fallg, PHI, ZHI, colker 32 33 INTEGER, SAVE :: ip, jp, kp, pstart, pend, psum 34 REAL, SAVE :: epsilon, urms 35 REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: gck, ec, ecf 36 REAL, DIMENSION(:), ALLOCATABLE, SAVE :: winf 37 REAL, SAVE :: eps2, urms2 31 PUBLIC colker, effic, fallg, phi, turbsd, turb_enhance_eff, zhi 32 33 INTEGER, SAVE :: ip, jp, kp, pend, pstart, psum 34 35 REAL, SAVE :: epsilon, eps2, urms, urms2 36 37 REAL, DIMENSION(:), ALLOCATABLE, SAVE :: winf 38 REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: ec, ecf, gck 38 39 39 40 ! 40 41 !-- Public interfaces 41 INTERFACE TurbSD 42 MODULE PROCEDURE TurbSD 43 END INTERFACE TurbSD 44 45 INTERFACE PHI 46 MODULE PROCEDURE PHI 47 END INTERFACE PHI 48 49 INTERFACE ZHI 50 MODULE PROCEDURE ZHI 51 END INTERFACE ZHI 42 INTERFACE colker 43 MODULE PROCEDURE colker 44 END INTERFACE colker 45 46 INTERFACE effic 47 MODULE PROCEDURE effic 48 END INTERFACE effic 49 50 INTERFACE fallg 51 MODULE PROCEDURE fallg 52 END INTERFACE fallg 53 54 INTERFACE phi 55 MODULE PROCEDURE phi 56 END INTERFACE phi 57 58 INTERFACE turbsd 59 MODULE PROCEDURE turbsd 60 END INTERFACE turbsd 52 61 53 62 INTERFACE turb_enhance_eff … … 55 64 END INTERFACE turb_enhance_eff 56 65 57 INTERFACE effic 58 MODULE PROCEDURE effic 59 END INTERFACE effic 60 61 INTERFACE fallg 62 MODULE PROCEDURE fallg 63 END INTERFACE fallg 64 65 INTERFACE colker 66 MODULE PROCEDURE colker 67 END INTERFACE colker 66 INTERFACE zhi 67 MODULE PROCEDURE zhi 68 END INTERFACE zhi 68 69 69 70 CONTAINS … … 72 73 ! SUBROUTINE for calculation of collision kernel 73 74 !------------------------------------------------------------------------------! 74 SUBROUTINE colker( i1,j1,k1,kernel)75 SUBROUTINE colker( i1, j1, k1, kernel ) 75 76 76 77 USE arrays_3d 77 78 USE cloud_parameters 78 79 USE constants 80 USE cpulog 79 81 USE indices 82 USE interfaces 80 83 USE particle_attributes 81 84 82 85 IMPLICIT NONE 83 86 84 INTEGER :: i,i1,j,j1,k187 INTEGER :: i, i1, j, j1, k1 85 88 86 89 REAL, DIMENSION(prt_start_index(k1,j1,i1):prt_start_index(k1,j1,i1)+prt_count(k1,j1,i1)-1,prt_start_index(k1,j1,i1):prt_start_index(k1,j1,i1)+prt_count(k1,j1,i1)-1) :: kernel 90 91 ! CALL cpu_log( log_point_s(46), 'colker', 'start' ) 87 92 88 93 ip = i1 … … 94 99 psum = prt_count(kp,jp,ip) 95 100 96 ALLOCATE( winf(pstart:pend), ec(pstart:pend,pstart:pend))97 98 IF ( turbulence_effects_on_collision ) THEN99 100 ALLOCATE( gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend))101 102 epsilon = diss(kp,jp,ip) * 1 0000.0!dissipation rate in cm**2/s**-3103 urms = 202.0 * ( epsilon/400.0)**(1.0/3.0)104 105 IF ( epsilon <= 0.001)THEN101 ALLOCATE( ec(pstart:pend,pstart:pend), winf(pstart:pend) ) 102 103 IF ( turbulence_effects_on_collision ) THEN 104 105 ALLOCATE( gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend) ) 106 107 epsilon = diss(kp,jp,ip) * 1.0E5 !dissipation rate in cm**2/s**-3 108 urms = 202.0 * ( epsilon/ 400.0 )**( 1.0 / 3.0 ) 109 110 IF ( epsilon <= 0.001 ) THEN 106 111 107 112 CALL fallg 108 113 CALL effic 109 114 110 DO j = pstart, pend, 1 111 DO i = pstart, pend, 1 112 kernel(i,j) = pi * (particles(j)%radius*100.0 + particles(i)%radius*100.0)**2.0 * ec(i,j) * abs(winf(j)-winf(i)) 115 DO j = pstart, pend 116 DO i = pstart, pend 117 kernel(i,j) = pi * ( particles(j)%radius * 100.0 + & 118 particles(i)%radius * 100.0 )**2 & 119 * ec(i,j) * ABS( winf(j) - winf(i) ) 113 120 ENDDO 114 121 ENDDO … … 116 123 ELSE 117 124 118 CALL TurbSD125 CALL turbsd 119 126 CALL turb_enhance_eff 120 127 CALL effic 121 128 122 DO j = pstart, pend, 1123 DO i = pstart, pend, 1129 DO j = pstart, pend, 1 130 DO i = pstart, pend, 1 124 131 kernel(i,j) = ec(i,j) * gck(i,j) * ecf(i,j) 125 132 ENDDO 126 133 ENDDO 134 127 135 ENDIF 128 136 … … 131 139 ELSE 132 140 141 ! CALL cpu_log( log_point_s(50), 'colker_fallg', 'start' ) 133 142 CALL fallg 143 ! CALL cpu_log( log_point_s(50), 'colker_fallg', 'stop' ) 144 ! CALL cpu_log( log_point_s(51), 'colker_effic', 'start' ) 134 145 CALL effic 135 136 DO j = pstart, pend, 1 137 DO i = pstart, pend, 1 138 kernel(i,j) = pi * (particles(j)%radius*100.0 + particles(i)%radius*100.0)**2.0 * ec(i,j) * abs(winf(j)-winf(i)) 146 ! CALL cpu_log( log_point_s(51), 'colker_effic', 'stop' ) 147 148 DO j = pstart, pend 149 DO i = pstart, pend 150 kernel(i,j) = pi * ( particles(j)%radius * 100.0 + & 151 particles(i)%radius * 100.0 )**2 & 152 * ec(i,j) * ABS( winf(j) - winf(i) ) 139 153 ENDDO 140 154 ENDDO … … 142 156 ENDIF 143 157 144 DEALLOCATE( winf, ec)145 146 RETURN 158 DEALLOCATE( ec, winf ) 159 160 ! CALL cpu_log( log_point_s(46), 'colker', 'stop' ) 147 161 148 162 END SUBROUTINE colker … … 151 165 ! SUBROUTINE for calculation of w, g and gck 152 166 !------------------------------------------------------------------------------! 153 SUBROUTINE TurbSD 154 !from Aayala 2008b, page 37ff, necessary input parameter water density, radii of droplets, air density, air viscosity, turbulent dissipation rate, taylor microscale reynolds number, gravitational acceleration 167 SUBROUTINE turbsd 168 ! from Aayala 2008b, page 37ff, necessary input parameter water density, radii 169 ! of droplets, air density, air viscosity, turbulent dissipation rate, 170 ! taylor microscale reynolds number, gravitational acceleration 155 171 USE constants 156 172 USE cloud_parameters … … 301 317 ENDDO 302 318 303 304 RETURN305 319 END SUBROUTINE TurbSD 306 320 … … 317 331 318 332 PHI = 1./aa1 - vsett/2.0/b/aa1**2.0 !in s 319 RETURN 333 320 334 END FUNCTION PHI 321 335 … … 339 353 + (2.*b/aa2 - 2.*b/aa1 - vsett1/aa2**2.0 + vsett2/aa1**2.0) & 340 354 * 1./2./b/aa3 ! in s**2 341 RETURN342 355 END FUNCTION ZHI 343 356 … … 430 443 ENDIF 431 444 ENDDO 432 RETURN 445 433 446 END SUBROUTINE fallg 434 447 … … 439 452 SUBROUTINE effic 440 453 441 USE constants 442 USE cloud_parameters 443 USE particle_attributes 444 USE arrays_3d 445 446 !collision efficiencies of hall kernel 454 USE arrays_3d 455 USE constants 456 USE cloud_parameters 457 USE particle_attributes 458 447 459 IMPLICIT NONE 448 460 449 461 INTEGER :: i, ir, iq, j, k, kk 450 REAL :: ek, rq, pp, qq 451 452 REAL, DIMENSION(1:21) :: rat 453 REAL, DIMENSION(1:15) :: r0 454 REAL, DIMENSION(1:15,1:21) :: ecoll 455 456 r0 = (/6.,8.,10.,15.,20.,25.,30.,40.,50., & 457 60.,70.,100.,150.,200.,300./) 458 rat = (/0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5, & 459 0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0/) 460 461 ecoll(:,1) = (/0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001/) 462 ecoll(:,2) = (/0.003,0.003,0.003,0.004,0.005,0.005,0.005,0.010,0.100,0.050,0.200,0.500,0.770,0.870,0.970/) 463 ecoll(:,3) = (/0.007,0.007,0.007,0.008,0.009,0.010,0.010,0.070,0.400,0.430,0.580,0.790,0.930,0.960,1.000/) 464 ecoll(:,4) = (/0.009,0.009,0.009,0.012,0.015,0.010,0.020,0.280,0.600,0.640,0.750,0.910,0.970,0.980,1.000/) 465 ecoll(:,5) = (/0.014,0.014,0.014,0.015,0.016,0.030,0.060,0.500,0.700,0.770,0.840,0.950,0.970,1.000,1.000/) 466 ecoll(:,6) = (/0.017,0.017,0.017,0.020,0.022,0.060,0.100,0.620,0.780,0.840,0.880,0.950,1.000,1.000,1.000/) 467 ecoll(:,7) = (/0.030,0.030,0.024,0.022,0.032,0.062,0.200,0.680,0.830,0.870,0.900,0.950,1.000,1.000,1.000/) 468 ecoll(:,8) = (/0.025,0.025,0.025,0.036,0.043,0.130,0.270,0.740,0.860,0.890,0.920,1.000,1.000,1.000,1.000/) 469 ecoll(:,9) = (/0.027,0.027,0.027,0.040,0.052,0.200,0.400,0.780,0.880,0.900,0.940,1.000,1.000,1.000,1.000/) 470 ecoll(:,10)= (/0.030,0.030,0.030,0.047,0.064,0.250,0.500,0.800,0.900,0.910,0.950,1.000,1.000,1.000,1.000/) 471 ecoll(:,11)= (/0.040,0.040,0.033,0.037,0.068,0.240,0.550,0.800,0.900,0.910,0.950,1.000,1.000,1.000,1.000/) 472 ecoll(:,12)= (/0.035,0.035,0.035,0.055,0.079,0.290,0.580,0.800,0.900,0.910,0.950,1.000,1.000,1.000,1.000/) 473 ecoll(:,13)= (/0.037,0.037,0.037,0.062,0.082,0.290,0.590,0.780,0.900,0.910,0.950,1.000,1.000,1.000,1.000/) 474 ecoll(:,14)= (/0.037,0.037,0.037,0.060,0.080,0.290,0.580,0.770,0.890,0.910,0.950,1.000,1.000,1.000,1.000/) 475 ecoll(:,15)= (/0.037,0.037,0.037,0.041,0.075,0.250,0.540,0.760,0.880,0.920,0.950,1.000,1.000,1.000,1.000/) 476 ecoll(:,16)= (/0.037,0.037,0.037,0.052,0.067,0.250,0.510,0.770,0.880,0.930,0.970,1.000,1.000,1.000,1.000/) 477 ecoll(:,17)= (/0.037,0.037,0.037,0.047,0.057,0.250,0.490,0.770,0.890,0.950,1.000,1.000,1.000,1.000,1.000/) 478 ecoll(:,18)= (/0.036,0.036,0.036,0.042,0.048,0.230,0.470,0.780,0.920,1.000,1.020,1.020,1.020,1.020,1.020/) 479 ecoll(:,19)= (/0.040,0.040,0.035,0.033,0.040,0.112,0.450,0.790,1.010,1.030,1.040,1.040,1.040,1.040,1.040/) 480 ecoll(:,20)= (/0.033,0.033,0.033,0.033,0.033,0.119,0.470,0.950,1.300,1.700,2.300,2.300,2.300,2.300,2.300/) 481 ecoll(:,21)= (/0.027,0.027,0.027,0.027,0.027,0.125,0.520,1.400,2.300,3.000,4.000,4.000,4.000,4.000,4.000/) 482 483 ! two-dimensional linear interpolation of the collision efficiency 484 ! radius has to be in µm 485 486 DO j = pstart, pend 487 DO i = pstart, j 488 489 DO k = 2, 15 490 IF ((particles(j)%radius*1.0E06).le.r0(k).and.(particles(j)%radius*1.0E06).ge.r0(k-1)) THEN 491 ir = k 492 ELSEIF ((particles(j)%radius*1.0E06).gt.r0(15)) THEN 493 ir = 16 494 ELSEIF ((particles(j)%radius*1.0E06).lt.r0(1)) THEN 495 ir = 1 496 ENDIF 497 ENDDO 498 499 rq = particles(i)%radius*1.0E06/(particles(j)%radius*1.0E06) 500 501 DO kk = 2, 21 502 IF (rq.le.rat(kk).and.rq.gt.rat(kk-1)) iq = kk 503 ENDDO 504 505 IF (ir.lt.16) THEN 506 IF (ir.ge.2) THEN 507 pp =((particles(j)%radius * 1.0E06) - r0(ir-1)) / (r0(ir)-r0(ir-1)) 508 qq =(rq-rat(iq-1))/(rat(iq)-rat(iq-1)) 509 ec(j,i) = (1.-pp) * (1.-qq) * ecoll(ir-1,iq-1) + & 510 pp * (1.-qq) * ecoll(ir,iq-1) + & 511 qq * (1.-pp) * ecoll(ir-1,iq) + & 512 pp * qq * ecoll(ir,iq) 462 463 INTEGER, DIMENSION(:), ALLOCATABLE :: ira 464 465 LOGICAL, SAVE :: first = .TRUE. 466 467 REAL :: ek, particle_radius, pp, qq, rq 468 469 REAL, DIMENSION(1:21), SAVE :: rat 470 REAL, DIMENSION(1:15), SAVE :: r0 471 REAL, DIMENSION(1:15,1:21), SAVE :: ecoll 472 473 ! 474 !-- Initial assignment of constants 475 IF ( first ) THEN 476 477 first = .FALSE. 478 r0 = (/6.,8.,10.,15.,20.,25.,30.,40.,50., 60.,70.,100.,150.,200., & 479 300./) 480 rat = (/0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6, & 481 0.65, 0.7,0.75,0.8,0.85,0.9,0.95,1.0/) 482 483 ecoll(:,1) = (/0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001, & 484 0.001,0.001,0.001,0.001,0.001,0.001/) 485 ecoll(:,2) = (/0.003,0.003,0.003,0.004,0.005,0.005,0.005,0.010,0.100, & 486 0.050,0.200,0.500,0.770,0.870,0.970/) 487 ecoll(:,3) = (/0.007,0.007,0.007,0.008,0.009,0.010,0.010,0.070,0.400, & 488 0.430,0.580,0.790,0.930,0.960,1.000/) 489 ecoll(:,4) = (/0.009,0.009,0.009,0.012,0.015,0.010,0.020,0.280,0.600, & 490 0.640,0.750,0.910,0.970,0.980,1.000/) 491 ecoll(:,5) = (/0.014,0.014,0.014,0.015,0.016,0.030,0.060,0.500,0.700, & 492 0.770,0.840,0.950,0.970,1.000,1.000/) 493 ecoll(:,6) = (/0.017,0.017,0.017,0.020,0.022,0.060,0.100,0.620,0.780, & 494 0.840,0.880,0.950,1.000,1.000,1.000/) 495 ecoll(:,7) = (/0.030,0.030,0.024,0.022,0.032,0.062,0.200,0.680,0.830, & 496 0.870,0.900,0.950,1.000,1.000,1.000/) 497 ecoll(:,8) = (/0.025,0.025,0.025,0.036,0.043,0.130,0.270,0.740,0.860, & 498 0.890,0.920,1.000,1.000,1.000,1.000/) 499 ecoll(:,9) = (/0.027,0.027,0.027,0.040,0.052,0.200,0.400,0.780,0.880, & 500 0.900,0.940,1.000,1.000,1.000,1.000/) 501 ecoll(:,10)= (/0.030,0.030,0.030,0.047,0.064,0.250,0.500,0.800,0.900, & 502 0.910,0.950,1.000,1.000,1.000,1.000/) 503 ecoll(:,11)= (/0.040,0.040,0.033,0.037,0.068,0.240,0.550,0.800,0.900, & 504 0.910,0.950,1.000,1.000,1.000,1.000/) 505 ecoll(:,12)= (/0.035,0.035,0.035,0.055,0.079,0.290,0.580,0.800,0.900, & 506 0.910,0.950,1.000,1.000,1.000,1.000/) 507 ecoll(:,13)= (/0.037,0.037,0.037,0.062,0.082,0.290,0.590,0.780,0.900, & 508 0.910,0.950,1.000,1.000,1.000,1.000/) 509 ecoll(:,14)= (/0.037,0.037,0.037,0.060,0.080,0.290,0.580,0.770,0.890, & 510 0.910,0.950,1.000,1.000,1.000,1.000/) 511 ecoll(:,15)= (/0.037,0.037,0.037,0.041,0.075,0.250,0.540,0.760,0.880, & 512 0.920,0.950,1.000,1.000,1.000,1.000/) 513 ecoll(:,16)= (/0.037,0.037,0.037,0.052,0.067,0.250,0.510,0.770,0.880, & 514 0.930,0.970,1.000,1.000,1.000,1.000/) 515 ecoll(:,17)= (/0.037,0.037,0.037,0.047,0.057,0.250,0.490,0.770,0.890, & 516 0.950,1.000,1.000,1.000,1.000,1.000/) 517 ecoll(:,18)= (/0.036,0.036,0.036,0.042,0.048,0.230,0.470,0.780,0.920, & 518 1.000,1.020,1.020,1.020,1.020,1.020/) 519 ecoll(:,19)= (/0.040,0.040,0.035,0.033,0.040,0.112,0.450,0.790,1.010, & 520 1.030,1.040,1.040,1.040,1.040,1.040/) 521 ecoll(:,20)= (/0.033,0.033,0.033,0.033,0.033,0.119,0.470,0.950,1.300, & 522 1.700,2.300,2.300,2.300,2.300,2.300/) 523 ecoll(:,21)= (/0.027,0.027,0.027,0.027,0.027,0.125,0.520,1.400,2.300, & 524 3.000,4.000,4.000,4.000,4.000,4.000/) 525 ENDIF 526 527 ! 528 !-- Calculate the radius class index of particles with respect to array r 529 ALLOCATE( ira(pstart:pend) ) 530 DO j = pstart, pend 531 particle_radius = particles(j)%radius * 1.0E6 532 DO k = 1, 15 533 IF ( particle_radius < r0(k) ) THEN 534 ira(j) = k 535 EXIT 536 ENDIF 537 ENDDO 538 IF ( particle_radius >= r0(15) ) ira(j) = 16 539 ENDDO 540 541 ! 542 !-- Two-dimensional linear interpolation of the collision efficiency. 543 !-- Radius has to be in µm 544 DO j = pstart, pend 545 DO i = pstart, j 546 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 554 ! ENDIF 555 ! ENDDO 556 557 ir = ira(j) 558 559 rq = particles(i)%radius / particles(j)%radius 560 561 ! DO kk = 2, 21 562 ! IF ( rq .le.rat(kk).and.rq.gt.rat(kk-1)) iq = kk 563 ! ENDDO 564 565 iq = INT( rq * 20 ) + 1 566 567 IF ( ir < 16 ) THEN 568 569 IF ( ir >= 2 ) THEN 570 pp = ( ( particles(j)%radius * 1.0E06 ) - r0(ir-1) ) / & 571 ( r0(ir) - r0(ir-1) ) 572 qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) ) 573 ec(j,i) = ( 1.0-pp ) * ( 1.0-qq ) * ecoll(ir-1,iq-1) & 574 + pp * ( 1.0-qq ) * ecoll(ir,iq-1) & 575 + qq * ( 1.0-pp ) * ecoll(ir-1,iq) & 576 + pp * qq * ecoll(ir,iq) 513 577 ELSE 514 578 qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1)) 515 579 ec(j,i) = (1.-qq) * ecoll(1,iq-1) + qq * ecoll(1,iq) 516 580 ENDIF 581 517 582 ELSE 518 qq = ( rq - rat(iq-1))/(rat(iq)-rat(iq-1))519 ek = ( 1.-qq)* ecoll(15,iq-1) + qq * ecoll(15,iq)520 ec(j,i) = min(ek,1.0)583 qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) ) 584 ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq) 585 ec(j,i) = MIN( ek, 1.0 ) 521 586 ENDIF 587 522 588 ec(i,j) = ec(j,i) 523 IF (ec(i,j).lt.1.e-20) stop 99 589 IF ( ec(i,j) < 1.0E-20 ) STOP 99 590 524 591 ENDDO 525 592 ENDDO 526 RETURN 593 594 DEALLOCATE( ira ) 595 527 596 END SUBROUTINE effic 528 597
Note: See TracChangeset
for help on using the changeset viewer.