source: palm/trunk/SOURCE/wang_kernel.f90 @ 821

Last change on this file since 821 was 800, checked in by franke, 13 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 24.8 KB
Line 
1MODULE wang_kernel_mod
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: wang_kernel.f90 800 2011-12-21 18:13:51Z maronga $
11!
12! 799 2011-12-21 17:48:03Z franke
13! speed optimizations and formatting
14! Bugfix: iq=1 is not allowed (routine effic)
15! Bugfix: replaced stop by ec=0.0 in case of very small ec (routine effic)
16!
17! 790 2011-11-29 03:11:20Z raasch
18! initial revision
19!
20! Description:
21! ------------
22! This routine calculates the effect of (SGS) turbulence on the collision
23! efficiency of droplets.
24! It is based on the original kernel developed by Wang (...)
25!------------------------------------------------------------------------------!
26
27    USE arrays_3d
28    USE cloud_parameters
29    USE constants
30    USE particle_attributes
31
32    IMPLICIT NONE
33
34    PRIVATE
35
36    PUBLIC  colker, effic, fallg, phi, turbsd, turb_enhance_eff, zhi 
37
38    INTEGER, SAVE ::  ip, jp, kp, pend, pstart, psum
39
40    REAL, SAVE :: epsilon, eps2, urms, urms2
41
42    REAL, DIMENSION(:),   ALLOCATABLE, SAVE ::  winf
43    REAL, DIMENSION(:,:), ALLOCATABLE, SAVE ::  ec, ecf, gck
44
45!
46!-- Public interfaces
47    INTERFACE colker
48       MODULE PROCEDURE colker
49    END INTERFACE colker
50
51    INTERFACE effic
52       MODULE PROCEDURE effic
53    END INTERFACE effic
54
55    INTERFACE fallg
56       MODULE PROCEDURE fallg
57    END INTERFACE fallg
58
59    INTERFACE phi
60       MODULE PROCEDURE phi
61    END INTERFACE phi
62
63    INTERFACE turbsd
64       MODULE PROCEDURE turbsd
65    END INTERFACE turbsd
66
67    INTERFACE turb_enhance_eff
68       MODULE PROCEDURE turb_enhance_eff
69    END INTERFACE turb_enhance_eff
70
71    INTERFACE zhi
72       MODULE PROCEDURE zhi
73    END INTERFACE zhi
74
75    CONTAINS
76
77!------------------------------------------------------------------------------!
78! SUBROUTINE for calculation of collision kernel
79!------------------------------------------------------------------------------!
80    SUBROUTINE colker( i1, j1, k1, kernel )
81
82       USE arrays_3d
83       USE cloud_parameters
84       USE constants
85       USE cpulog
86       USE indices
87       USE interfaces
88       USE particle_attributes
89
90       IMPLICIT NONE
91
92       INTEGER ::  i, i1, j, j1, k1
93
94       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
95
96!       CALL cpu_log( log_point_s(46), 'colker', 'start' )
97
98       ip = i1
99       jp = j1
100       kp = k1
101
102       pstart = prt_start_index(kp,jp,ip)
103       pend   = prt_start_index(kp,jp,ip) + prt_count(kp,jp,ip) - 1
104       psum   = prt_count(kp,jp,ip)
105
106       ALLOCATE( ec(pstart:pend,pstart:pend), winf(pstart:pend) )
107
108       IF ( turbulence_effects_on_collision )  THEN
109
110          ALLOCATE( gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend) )
111
112          epsilon = diss(kp,jp,ip) * 1.0E5  !dissipation rate in cm**2/s**-3
113          urms     = 202.0 * ( epsilon/ 400.0 )**( 1.0 / 3.0 )
114
115          IF ( epsilon <= 0.001 )  THEN
116
117             CALL fallg
118             CALL effic
119
120             DO  j = pstart, pend
121                DO  i =  pstart, pend
122                   kernel(i,j) = pi * ( particles(j)%radius * 100.0 +    &
123                                        particles(i)%radius * 100.0 )**2 &
124                                    * ec(i,j) * ABS( winf(j) - winf(i) )
125                ENDDO
126             ENDDO
127
128          ELSE
129
130             CALL turbsd
131             CALL turb_enhance_eff
132             CALL effic
133
134             DO  j = pstart, pend, 1
135                DO  i =  pstart, pend, 1
136                   kernel(i,j) = ec(i,j) * gck(i,j) * ecf(i,j)
137                ENDDO
138             ENDDO
139
140          ENDIF
141
142          DEALLOCATE(gck, ecf)
143
144       ELSE
145
146!          CALL cpu_log( log_point_s(50), 'colker_fallg', 'start' )
147          CALL fallg
148!          CALL cpu_log( log_point_s(50), 'colker_fallg', 'stop' )
149!          CALL cpu_log( log_point_s(51), 'colker_effic', 'start' )
150          CALL effic
151!          CALL cpu_log( log_point_s(51), 'colker_effic', 'stop' )
152
153          DO  j = pstart, pend
154             DO  i =  pstart, pend
155                kernel(i,j) = pi * ( particles(j)%radius * 100.0 +    &
156                                     particles(i)%radius * 100.0 )**2 &
157                                 * ec(i,j) * ABS( winf(j) - winf(i) )
158             ENDDO
159          ENDDO
160
161       ENDIF
162
163       DEALLOCATE( ec, winf )
164
165!       CALL cpu_log( log_point_s(46), 'colker', 'stop' )
166
167    END SUBROUTINE colker
168
169!------------------------------------------------------------------------------!
170! SUBROUTINE for calculation of w, g and gck
171!------------------------------------------------------------------------------!
172    SUBROUTINE turbsd
173! from Aayala 2008b, page 37ff, necessary input parameter water density, radii
174! of droplets, air density, air viscosity, turbulent dissipation rate,
175! taylor microscale reynolds number, gravitational acceleration
176
177       USE constants
178       USE cloud_parameters
179       USE particle_attributes
180       USE arrays_3d
181
182       IMPLICIT NONE
183
184       REAL :: Relamda, &
185               Tl, Lf, tauk, eta, vk, ao, lambda, tt, z, be,       &
186               bbb, d1, e1, d2, e2, ccc, b1, c1, b2, c2, v1xysq,  &
187               vrms1xy, v2xysq, vrms2xy, v1v2xy, fR, wrtur2xy, wrgrav2, &
188               wrFIN, SSt, XX, YY, c1_gr, ao_gr, fao_gr, rc, grFIN, v1, t1, v2, t2, rrp
189
190       REAL, SAVE :: airvisc, airdens, anu, gravity, waterdens
191
192       REAL, DIMENSION(pstart:pend) :: St, tau
193
194       LOGICAL, SAVE ::  first = .TRUE.
195
196       INTEGER :: i, j
197
198!
199!--   Initial assignment of constants
200       IF ( first )  THEN
201
202          first = .FALSE.
203          airvisc   = 0.1818   !dynamic viscosity in mg/cm*s
204          airdens   = 1.2250   !air density in mg/cm**3
205          waterdens = 1000.0   !water density in mg/cm**3
206          gravity   = 980.6650 !in cm/s**2
207          anu = airvisc/airdens  ! kinetic viscosity in cm**2/s
208
209       ENDIF   
210
211       Relamda = urms**2*sqrt(15.0/epsilon/anu)
212
213       Tl = urms**2/epsilon       !in s
214       Lf = 0.5 * (urms**3)/epsilon !in cm
215
216       tauk = (anu/epsilon)**0.5       !in s
217       eta  = (anu**3/epsilon)**0.25  !in cm
218       vk   = eta/tauk                   !in cm/s
219
220       ao = (11.+7.*Relamda)/(205.+Relamda)
221
222       lambda = urms * sqrt(15.*anu/epsilon)     !in cm
223
224       tt = sqrt(2.*Relamda/(15.**0.5)/ao) * tauk !in s
225
226       CALL fallg !gives winf in cm/s
227
228       DO i = pstart, pend
229          tau(i) = winf(i)/gravity    !in s
230          St(i)  = tau(i)/tauk
231       ENDDO
232
233!***** TO CALCULATE wr ********************************
234!from Aayala 2008b, page 38f
235
236       z  = tt/Tl
237       be = sqrt(2.0)*lambda/Lf
238
239       bbb = sqrt(1.0-2.0*be**2)
240       d1 = (1.+bbb)/2.0/bbb
241       e1 = Lf*(1.0+bbb)/2.0   !in cm
242       d2 = (1.0-bbb)/2.0/bbb
243       e2 = Lf*(1.0-bbb)/2.0   !in cm
244
245       ccc = sqrt(1.0-2.0*z**2)
246       b1 = (1.+ccc)/2./ccc
247       c1 = Tl*(1.+ccc)/2.   !in s
248       b2 = (1.-ccc)/2./ccc
249       c2 = Tl*(1.-ccc)/2.   !in s
250
251       DO i = pstart, pend
252          v1 = winf(i)         !in cm/s
253          t1 = tau(i)         !in s
254
255          DO j = pstart,i
256             rrp = (particles(i)%radius + particles(j)%radius) * 100.0 !radius in cm
257             v2  = winf(j)         !in cm/s
258             t2  = tau(j)         !in s
259
260             v1xysq =  b1*d1*PHI(c1,e1,v1,t1) &
261                     - b1*d2*PHI(c1,e2,v1,t1) &
262                     - b2*d1*PHI(c2,e1,v1,t1) &
263                     + b2*d2*PHI(c2,e2,v1,t1)
264             v1xysq = v1xysq * urms**2/t1    !in cm**2/s**2
265
266             vrms1xy= sqrt(v1xysq)             !in cm/s
267
268             v2xysq =  b1*d1*PHI(c1,e1,v2,t2) &
269                     - b1*d2*PHI(c1,e2,v2,t2) &
270                     - b2*d1*PHI(c2,e1,v2,t2) &
271                     + b2*d2*PHI(c2,e2,v2,t2)
272             v2xysq = v2xysq * urms**2/t2    !in cm**2/s**2
273
274              vrms2xy= sqrt(v2xysq)             !in cm/s
275
276             IF(winf(i).ge.winf(j)) THEN
277                v1 = winf(i)
278                t1 = tau(i)
279                v2 = winf(j)
280                t2 = tau(j)
281             ELSE
282                v1 = winf(j)
283                t1 = tau(j)
284                v2 = winf(i)
285                t2 = tau(i)
286             ENDIF
287
288             v1v2xy =  b1*d1*ZHI(c1,e1,v1,t1,v2,t2) &
289                     - b1*d2*ZHI(c1,e2,v1,t1,v2,t2) &
290                     - b2*d1*ZHI(c2,e1,v1,t1,v2,t2) &
291                     + b2*d2*ZHI(c2,e2,v1,t1,v2,t2)
292             fR = d1 * exp(-rrp/e1) - d2 * exp(-rrp/e2)
293             v1v2xy = v1v2xy * fR * urms**2/tau(i)/tau(j)    !in cm**2/s**2
294
295             wrtur2xy=vrms1xy**2 + vrms2xy**2 - 2.*v1v2xy  !in cm**2/s**2
296
297              IF (wrtur2xy.le.0.0) wrtur2xy=0.0
298
299              wrgrav2=pi/8.*(winf(j)-winf(i))**2
300
301              wrFIN = sqrt((2.0/pi)*(wrtur2xy+wrgrav2))    !in cm/s
302
303
304!***** TO CALCULATE gr ********************************
305
306             IF(St(j).gt.St(i)) THEN
307                SSt = St(j)
308             ELSE
309                SSt = St(i)
310             ENDIF
311
312             XX = -0.1988*SSt**4 + 1.5275*SSt**3 &
313                  -4.2942*SSt**2 + 5.3406*SSt
314
315             IF(XX.le.0.0) XX = 0.0
316
317             YY = 0.1886*exp(20.306/Relamda)
318
319             c1_gr = XX/(gravity/(vk/tauk))**YY
320
321             ao_gr  = ao + (pi/8.)*(gravity/(vk/tauk))**2
322             fao_gr = 20.115 * (ao_gr/Relamda)**0.5
323             rc     = sqrt( fao_gr * abs(St(j)-St(i)) ) * eta   !in cm
324
325             grFIN = ((eta**2+rc**2)/(rrp**2+rc**2))**(c1_gr/2.)
326             IF (grFIN.lt.1.0) grFIN = 1.0
327
328             gck(i,j) = 2. * pi * rrp**2 * wrFIN * grFIN       ! in cm**3/s
329             gck(j,i) = gck(i,j)
330
331          ENDDO
332       ENDDO
333
334    END SUBROUTINE TurbSD
335
336!------------------------------------------------------------------------------!
337! PHI as a function
338!------------------------------------------------------------------------------!
339  REAL FUNCTION PHI(a,b,vsett,tau0)
340
341       IMPLICIT NONE
342
343       REAL :: a, aa1, b, vsett, tau0
344
345       aa1 = 1./tau0 + 1./a + vsett/b
346
347       PHI = 1./aa1 - vsett/2.0/b/aa1**2  !in s
348
349    END FUNCTION PHI
350
351!------------------------------------------------------------------------------!
352! ZETA as a function
353!------------------------------------------------------------------------------!
354  REAL FUNCTION ZHI(a,b,vsett1,tau1,vsett2,tau2)
355
356       IMPLICIT NONE
357
358       REAL :: a, aa1, aa2, aa3, aa4, aa5, aa6, b, vsett1, tau1, vsett2, tau2
359
360        aa1 = vsett2/b - 1./tau2 - 1./a
361        aa2 = vsett1/b + 1./tau1 + 1./a
362        aa3 = (vsett1-vsett2)/b + 1./tau1 + 1./tau2
363        aa4 = (vsett2/b)**2 - (1./tau2 + 1./a)**2
364        aa5 = vsett2/b + 1./tau2 + 1./a
365        aa6 = 1./tau1 - 1./a + (1./tau2 + 1./a) * vsett1/vsett2
366        ZHI =  (1./aa1 - 1./aa2) * (vsett1-vsett2)/2./b/aa3**2          &
367             + (4./aa4 - 1./aa5**2 - 1./aa1**2) * vsett2/2./b/aa6      &
368             + (2.*b/aa2 - 2.*b/aa1 - vsett1/aa2**2 + vsett2/aa1**2)   &
369                                                           * 1./2./b/aa3             ! in s**2
370
371    END FUNCTION ZHI
372
373!------------------------------------------------------------------------------!
374! SUBROUTINE for calculation of terminal velocity winf
375!------------------------------------------------------------------------------!
376   SUBROUTINE fallg
377
378       USE constants
379       USE cloud_parameters
380       USE particle_attributes
381       USE arrays_3d
382
383      IMPLICIT NONE
384
385      INTEGER :: i, j
386
387      LOGICAL, SAVE ::  first = .TRUE.
388
389      REAL, SAVE :: eta, xlamb, rhoa, rhow, grav, cunh, t0, sigma, stok, stb, phy, py
390
391      REAL :: bond, x, xrey, y
392
393      REAL, DIMENSION(1:7), SAVE  :: b
394      REAL, DIMENSION(1:6), SAVE  :: c
395
396!
397!--   Initial assignment of constants
398      IF ( first )  THEN
399
400         first = .FALSE.
401         b = (/-0.318657e1,0.992696e0,-0.153193e-2,-0.987059e-3, &
402             -0.578878e-3,0.855176e-4,-0.327815e-5/)
403         c = (/-0.500015e1,0.523778e1,-0.204914e1,0.475294e0,    &
404             -0.542819e-1, 0.238449e-2/)
405
406         eta   = 1.818e-4     !in poise = g/(cm s)
407         xlamb = 6.62e-6    !in cm
408
409         rhow = 1.0       !in g/cm**3
410         rhoa = 1.225e-3    !in g/cm**3
411
412         grav = 980.665     !in cm/s**2
413         cunh = 1.257 * xlamb !in cm
414         t0   = 273.15        !in K
415         sigma= 76.1 - 0.155 * (293.15 - t0) !in N/m = g/s**2
416         stok = 2.0 * grav * (rhow - rhoa)/(9.0 * eta) ! in 1/(cm s)
417         stb  = 32.0 * rhoa * (rhow-rhoa) * grav/(3.0 * eta * eta) ! in 1/cm**3
418         phy  = (sigma**3) * (rhoa**2)/((eta**4) * grav * (rhow-rhoa))
419         py   = phy**(1.0/6.0)
420
421      ENDIF
422
423!particle radius has to be in cm
424      DO j = pstart, pend
425
426         IF (particles(j)%radius*100.0 .le. 1.e-3) THEN
427
428            winf(j)=stok*((particles(j)%radius*100.0)**2+cunh* particles(j)%radius*100.0) !in cm/s
429
430         ELSEIF (particles(j)%radius*100.0.gt.1.e-3.and.particles(j)%radius*100.0.le.5.35e-2) THEN
431
432            x = log(stb*(particles(j)%radius*100.0)**3)
433            y = 0.0
434
435            DO i = 1, 7
436               y = y + b(i) * (x**(i-1))
437            ENDDO
438
439            xrey = (1.0 + cunh/(particles(j)%radius*100.0)) * exp(y)
440            winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0)  !in cm/s
441
442         ELSEIF (particles(j)%radius*100.0.gt.5.35e-2) THEN
443
444            IF (particles(j)%radius*100.0.gt.0.35)  THEN
445               bond = grav*(rhow-rhoa) * 0.35 * 0.35/sigma
446            ELSE
447               bond = grav*(rhow-rhoa)*(particles(j)%radius*100.0)**2/sigma
448            ENDIF
449
450            x = log(16.0*bond*py/3.0)
451            y = 0.0
452
453            DO i = 1, 6
454               y = y + c(i) * (x**(i-1))
455            ENDDO
456
457            xrey = py*exp(y)
458
459            IF (particles(j)%radius*100.0 .gt.0.35)  THEN
460              winf(j) = xrey * eta/(2.0 * rhoa * 0.35) !in cm/s
461            ELSE
462               winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0)  !in cm/s
463            ENDIF
464
465         ENDIF
466      ENDDO
467      RETURN
468   END SUBROUTINE fallg
469
470!------------------------------------------------------------------------------!
471! SUBROUTINE for calculation of collision efficencies
472!------------------------------------------------------------------------------!
473
474   SUBROUTINE effic
475
476      USE arrays_3d
477      USE constants
478      USE cloud_parameters
479      USE particle_attributes
480
481!collision efficiencies of hall kernel
482      IMPLICIT NONE
483
484      INTEGER :: i, ir, iq, j, k, kk
485
486      INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
487
488      LOGICAL, SAVE ::  first = .TRUE.
489
490      REAL ::  ek, particle_radius, pp, qq, rq
491
492      REAL, DIMENSION(1:21), SAVE :: rat
493      REAL, DIMENSION(1:15), SAVE :: r0
494      REAL, DIMENSION(1:15,1:21), SAVE :: ecoll
495
496!
497!--   Initial assignment of constants
498      IF ( first )  THEN
499
500         first = .FALSE.
501         r0  = (/6.,8.,10.,15.,20.,25.,30.,40.,50., 60.,70.,100.,150.,200., &
502                 300./)
503         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,  &
504                 0.65, 0.7,0.75,0.8,0.85,0.9,0.95,1.0/)
505
506         ecoll(:,1) = (/0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001, &
507                        0.001,0.001,0.001,0.001,0.001,0.001/)
508         ecoll(:,2) = (/0.003,0.003,0.003,0.004,0.005,0.005,0.005,0.010,0.100, &
509                        0.050,0.200,0.500,0.770,0.870,0.970/)
510         ecoll(:,3) = (/0.007,0.007,0.007,0.008,0.009,0.010,0.010,0.070,0.400, &
511                        0.430,0.580,0.790,0.930,0.960,1.000/)
512         ecoll(:,4) = (/0.009,0.009,0.009,0.012,0.015,0.010,0.020,0.280,0.600, &
513                        0.640,0.750,0.910,0.970,0.980,1.000/)
514         ecoll(:,5) = (/0.014,0.014,0.014,0.015,0.016,0.030,0.060,0.500,0.700, &
515                        0.770,0.840,0.950,0.970,1.000,1.000/)
516         ecoll(:,6) = (/0.017,0.017,0.017,0.020,0.022,0.060,0.100,0.620,0.780, &
517                        0.840,0.880,0.950,1.000,1.000,1.000/)
518         ecoll(:,7) = (/0.030,0.030,0.024,0.022,0.032,0.062,0.200,0.680,0.830, &
519                        0.870,0.900,0.950,1.000,1.000,1.000/)
520         ecoll(:,8) = (/0.025,0.025,0.025,0.036,0.043,0.130,0.270,0.740,0.860, &
521                        0.890,0.920,1.000,1.000,1.000,1.000/)
522         ecoll(:,9) = (/0.027,0.027,0.027,0.040,0.052,0.200,0.400,0.780,0.880, &
523                        0.900,0.940,1.000,1.000,1.000,1.000/)
524         ecoll(:,10)= (/0.030,0.030,0.030,0.047,0.064,0.250,0.500,0.800,0.900, &
525                        0.910,0.950,1.000,1.000,1.000,1.000/)
526         ecoll(:,11)= (/0.040,0.040,0.033,0.037,0.068,0.240,0.550,0.800,0.900, &
527                        0.910,0.950,1.000,1.000,1.000,1.000/)
528         ecoll(:,12)= (/0.035,0.035,0.035,0.055,0.079,0.290,0.580,0.800,0.900, &
529                        0.910,0.950,1.000,1.000,1.000,1.000/)
530         ecoll(:,13)= (/0.037,0.037,0.037,0.062,0.082,0.290,0.590,0.780,0.900, &
531                        0.910,0.950,1.000,1.000,1.000,1.000/)
532         ecoll(:,14)= (/0.037,0.037,0.037,0.060,0.080,0.290,0.580,0.770,0.890, &
533                        0.910,0.950,1.000,1.000,1.000,1.000/)
534         ecoll(:,15)= (/0.037,0.037,0.037,0.041,0.075,0.250,0.540,0.760,0.880, &
535                        0.920,0.950,1.000,1.000,1.000,1.000/)
536         ecoll(:,16)= (/0.037,0.037,0.037,0.052,0.067,0.250,0.510,0.770,0.880, &
537                        0.930,0.970,1.000,1.000,1.000,1.000/)
538         ecoll(:,17)= (/0.037,0.037,0.037,0.047,0.057,0.250,0.490,0.770,0.890, &
539                        0.950,1.000,1.000,1.000,1.000,1.000/)
540         ecoll(:,18)= (/0.036,0.036,0.036,0.042,0.048,0.230,0.470,0.780,0.920, &
541                        1.000,1.020,1.020,1.020,1.020,1.020/)
542         ecoll(:,19)= (/0.040,0.040,0.035,0.033,0.040,0.112,0.450,0.790,1.010, &
543                        1.030,1.040,1.040,1.040,1.040,1.040/)
544         ecoll(:,20)= (/0.033,0.033,0.033,0.033,0.033,0.119,0.470,0.950,1.300, &
545                        1.700,2.300,2.300,2.300,2.300,2.300/)
546         ecoll(:,21)= (/0.027,0.027,0.027,0.027,0.027,0.125,0.520,1.400,2.300, &
547                        3.000,4.000,4.000,4.000,4.000,4.000/)
548      ENDIF
549
550!
551!--   Calculate the radius class index of particles with respect to array r
552      ALLOCATE( ira(pstart:pend) )
553      DO  j = pstart, pend
554         particle_radius = particles(j)%radius * 1.0E6
555         DO  k = 1, 15
556            IF ( particle_radius < r0(k) )  THEN
557               ira(j) = k
558               EXIT
559            ENDIF
560         ENDDO
561         IF ( particle_radius >= r0(15) )  ira(j) = 16
562      ENDDO
563
564!
565!--   Two-dimensional linear interpolation of the collision efficiency.
566!--   Radius has to be in µm
567      DO  j = pstart, pend
568         DO  i = pstart, j
569
570            ir = ira(j)
571
572            rq = particles(i)%radius / particles(j)%radius
573
574!            DO kk = 2, 21
575!               IF ( rq <= rat(kk) )  THEN
576!                  iq = kk
577!                  EXIT
578!               ENDIF
579!            ENDDO
580
581            iq = INT( rq * 20 ) + 1
582            iq = MAX(iq , 2)
583
584            IF ( ir < 16 )  THEN
585
586               IF ( ir >= 2 )  THEN
587                  pp = ( ( particles(j)%radius * 1.0E06 ) - r0(ir-1) ) / &
588                       ( r0(ir) - r0(ir-1) )
589                  qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
590                  ec(j,i) = ( 1.0-pp ) * ( 1.0-qq ) * ecoll(ir-1,iq-1)   &
591                            + pp * ( 1.0-qq ) * ecoll(ir,iq-1)           &
592                            + qq * ( 1.0-pp ) * ecoll(ir-1,iq)           &
593                            + pp * qq * ecoll(ir,iq)
594               ELSE
595                  qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
596                  ec(j,i) = (1.-qq) * ecoll(1,iq-1) + qq * ecoll(1,iq)
597               ENDIF
598
599            ELSE
600               qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
601               ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq)
602               ec(j,i) = MIN( ek, 1.0 )
603            ENDIF
604
605            ec(i,j) = ec(j,i)
606            IF ( ec(i,j) < 1.0E-20 )  ec(i,j) = 0.0
607
608         ENDDO
609      ENDDO
610
611      DEALLOCATE( ira )
612
613   END SUBROUTINE effic
614
615!------------------------------------------------------------------------------!
616! SUBROUTINE for calculation of enhancement factor collision efficencies
617!------------------------------------------------------------------------------!
618   SUBROUTINE turb_enhance_eff
619
620       USE constants
621       USE cloud_parameters
622       USE particle_attributes
623       USE arrays_3d
624
625      IMPLICIT NONE
626
627      INTEGER :: i, ik, ir, iq, j, k, kk
628
629      INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
630
631      REAL    :: rq, y1, particle_radius, pp, qq, y2, y3, x1, x2, x3
632
633      LOGICAL, SAVE ::  first = .TRUE.
634
635      REAL, DIMENSION(1:11), SAVE :: rat
636      REAL, DIMENSION(1:7), SAVE  :: r0
637      REAL, DIMENSION(1:7,1:11), SAVE :: ecoll_100, ecoll_400
638
639!
640!--   Initial assignment of constants
641      IF ( first )  THEN
642
643         first = .FALSE.
644
645         r0  = (/10., 20., 30.,40., 50., 60.,100./)
646         rat = (/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0/)
647
648! 100 cm^2/s^3
649         ecoll_100(:,1) = (/1.74,  1.74,  1.773, 1.49,  1.207,  1.207,  1.0 /)
650         ecoll_100(:,2) = (/1.46,  1.46,  1.421, 1.245, 1.069,  1.069,  1.0 /)
651         ecoll_100(:,3) = (/1.32,  1.32,  1.245, 1.123, 1.000,  1.000,  1.0 /)
652         ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025,  1.025,  1.0 /)
653         ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056,  1.056,  1.0 /)
654         ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028,  1.028,  1.0 /)
655         ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046,  1.046,  1.0 /)
656         ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029,  1.029,  1.0 /)
657         ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021,  1.021,  1.0 /)
658         ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088,  1.088,  1.0 /)
659         ecoll_100(:,11)= (/20.3,  20.3,  14.6 , 8.61,  2.60,   2.60 ,  1.0 /)
660
661! 400 cm^2/s^3
662         ecoll_400(:,1) = (/4.976, 4.976,  3.593, 2.519, 1.445,  1.445,  1.0 /)
663         ecoll_400(:,2) = (/2.984, 2.984,  2.181, 1.691, 1.201,  1.201,  1.0 /)
664         ecoll_400(:,3) = (/1.988, 1.988,  1.475, 1.313, 1.150,  1.150,  1.0 /)
665         ecoll_400(:,4) = (/1.490, 1.490,  1.187, 1.156, 1.126,  1.126,  1.0 /)
666         ecoll_400(:,5) = (/1.249, 1.249,  1.088, 1.090, 1.092,  1.092,  1.0 /)
667         ecoll_400(:,6) = (/1.139, 1.139,  1.130, 1.091, 1.051,  1.051,  1.0 /)
668         ecoll_400(:,7) = (/1.220, 1.220,  1.190, 1.138, 1.086,  1.086,  1.0 /)
669         ecoll_400(:,8) = (/1.325, 1.325,  1.267, 1.165, 1.063,  1.063,  1.0 /)
670         ecoll_400(:,9) = (/1.716, 1.716,  1.345, 1.223, 1.100,  1.100,  1.0 /)
671         ecoll_400(:,10)= (/3.788, 3.788,  1.501, 1.311, 1.120,  1.120,  1.0 /)
672         ecoll_400(:,11)= (/36.52, 36.52,  19.16, 22.80,  26.0,   26.0,  1.0 /)
673
674      ENDIF
675
676!
677!--   Calculate the radius class index of particles with respect to array r
678      ALLOCATE( ira(pstart:pend) )
679
680      DO  j = pstart, pend
681         particle_radius = particles(j)%radius * 1.0E6
682         DO  k = 1, 7
683            IF ( particle_radius < r0(k) )  THEN
684               ira(j) = k
685               EXIT
686            ENDIF
687         ENDDO
688         IF ( particle_radius >= r0(7) )  ira(j) = 8
689      ENDDO
690
691! two-dimensional linear interpolation of the collision efficiency
692      DO j =  pstart, pend
693         DO i = pstart, j
694
695            ir = ira(j)
696
697            rq = particles(i)%radius/particles(j)%radius
698
699            DO kk = 2, 11
700               IF ( rq <= rat(kk) )  THEN
701                  iq = kk
702                  EXIT
703               ENDIF
704            ENDDO
705
706! 0  cm2/s3
707            y1 = 1.0
708! 100 cm2/s3, 400 cm2/s3
709            IF (ir.lt.8) THEN
710               IF (ir.ge.2) THEN
711                  pp = ((particles(j)%radius*1.0E06)-r0(ir-1))/(r0(ir)-r0(ir-1))
712                  qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
713                  y2= (1.-pp)*(1.-qq)*ecoll_100(ir-1,iq-1)+  &
714                            pp*(1.-qq)*ecoll_100(ir,iq-1)+    &
715                            qq*(1.-pp)*ecoll_100(ir-1,iq)+    &
716                            pp*qq*ecoll_100(ir,iq)
717                  y3= (1.-pp)*(1.-qq)*ecoll_400(ir-1,iq-1)+  &
718                            pp*(1.-qq)*ecoll_400(ir,iq-1)+    &
719                            qq*(1.-pp)*ecoll_400(ir-1,iq)+    &
720                            pp*qq*ecoll_400(ir,iq)
721               ELSE
722                  qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
723                  y2= (1.-qq)*ecoll_100(1,iq-1)+qq*ecoll_100(1,iq)
724                  y3= (1.-qq)*ecoll_400(1,iq-1)+qq*ecoll_400(1,iq)
725               ENDIF
726            ELSE
727            qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
728            y2 = (1.-qq) * ecoll_100(7,iq-1) + qq * ecoll_100(7,iq)
729            y3 = (1.-qq) * ecoll_400(7,iq-1) + qq * ecoll_400(7,iq)
730         ENDIF
731! linear interpolation
732! dissipation rate in cm2/s3
733         x1 = 0.0
734         x2 = 100.0
735         x3 = 400.0
736
737         IF (epsilon.le.100.) THEN
738            ecf(j,i) = (epsilon-100.)/(0.-100.) * y1 &
739                    +  (epsilon-0.)/(100.-0.) * y2
740         ELSE IF(epsilon.le.600.)THEN
741            ecf(j,i) = (epsilon-400.)/(100.-400.) * y2 &
742                    +  (epsilon-100.)/(400.-100.) * y3
743
744         ELSE
745            ecf(j,i) = (600.-400.)/(100.-400.) * y2 &
746                    +  (600.-100.)/(400.-100.) * y3
747         ENDIF
748
749         IF (ecf(j,i).lt.1.0) ecf(j,i) = 1.0
750
751         ecf(i,j)=ecf(j,i)
752      ENDDO
753   ENDDO
754
755   RETURN
756   END SUBROUTINE turb_enhance_eff
757
758 END MODULE wang_kernel_mod
Note: See TracBrowser for help on using the repository browser.