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

Last change on this file since 791 was 791, checked in by raasch, 13 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 22.4 KB
Line 
1MODULE wang_kernel_mod
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: wang_kernel.f90 791 2011-11-29 03:33:42Z raasch $
11!
12! 790 2011-11-29 03:11:20Z raasch
13! initial revision
14!
15! Description:
16! ------------
17! This routine calculates the effect of (SGS) turbulence on the collision
18! efficiency of droplets.
19! It is based on the original kernel developed by Wang (...)
20!------------------------------------------------------------------------------!
21
22    USE arrays_3d
23    USE cloud_parameters
24    USE constants
25    USE particle_attributes
26
27    IMPLICIT NONE
28
29    PRIVATE
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
38
39!
40!-- 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
52
53    INTERFACE turb_enhance_eff
54       MODULE PROCEDURE turb_enhance_eff
55    END INTERFACE turb_enhance_eff
56
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
68
69    CONTAINS
70
71!------------------------------------------------------------------------------!
72! SUBROUTINE for calculation of collision kernel
73!------------------------------------------------------------------------------!
74    SUBROUTINE colker(i1,j1,k1,kernel)
75
76       USE arrays_3d
77       USE cloud_parameters
78       USE constants
79       USE indices
80       USE particle_attributes
81
82       IMPLICIT NONE
83
84       INTEGER :: i,i1,j,j1,k1
85
86       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
87
88       ip = i1
89       jp = j1
90       kp = k1
91
92       pstart = prt_start_index(kp,jp,ip)
93       pend   = prt_start_index(kp,jp,ip) + prt_count(kp,jp,ip) - 1
94       psum   = prt_count(kp,jp,ip)
95
96       ALLOCATE(winf(pstart:pend), ec(pstart:pend,pstart:pend))
97
98       IF ( turbulence_effects_on_collision ) THEN
99
100          ALLOCATE(gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend))
101
102          epsilon = diss(kp,jp,ip) * 10000.0  !dissipation rate in cm**2/s**-3
103          urms     = 202.0 * (epsilon/400.0)**(1.0/3.0)
104
105          IF (epsilon <= 0.001) THEN
106
107             CALL fallg
108             CALL effic
109
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))
113                ENDDO
114             ENDDO
115
116          ELSE
117
118             CALL TurbSD
119             CALL turb_enhance_eff
120             CALL effic
121
122             DO j = pstart, pend, 1
123                DO i =  pstart, pend, 1
124                   kernel(i,j) = ec(i,j) * gck(i,j) * ecf(i,j)
125                ENDDO
126             ENDDO
127          ENDIF
128
129          DEALLOCATE(gck, ecf)
130
131       ELSE
132
133          CALL fallg
134          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))
139             ENDDO
140          ENDDO
141
142       ENDIF
143
144       DEALLOCATE(winf, ec)
145
146       RETURN
147
148    END SUBROUTINE colker
149
150!------------------------------------------------------------------------------!
151! SUBROUTINE for calculation of w, g and gck
152!------------------------------------------------------------------------------!
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
155       USE constants
156       USE cloud_parameters
157       USE particle_attributes
158       USE arrays_3d
159
160       IMPLICIT NONE
161
162       REAL :: airvisc, airdens, waterdens, gravity, anu, Relamda, &
163               Tl, Lf, tauk, eta, vk, ao, lambda, tt, z, be,       &
164               bbb, d1, e1, d2, e2, ccc, b1, c1, b2, c2, v1xysq,  &
165               vrms1xy, v2xysq, vrms2xy, v1v2xy, fR, wrtur2xy, wrgrav2, &
166               wrFIN, SSt, XX, YY, c1_gr, ao_gr, fao_gr, rc, grFIN, v1, t1, v2, t2, rrp
167
168       REAL, DIMENSION(pstart:pend) :: vST, tau, St
169
170       INTEGER :: i, j
171
172       airvisc   = 0.1818   !dynamic viscosity in mg/cm*s
173       airdens   = 1.2250   !air density in mg/cm**3
174       waterdens = 1000.0   !water density in mg/cm**3
175       gravity   = 980.6650 !in cm/s**2
176
177       anu = airvisc/airdens  ! kinetic viscosity in cm**2/s
178
179       Relamda = urms**2.0*sqrt(15.0/epsilon/anu)
180
181       Tl = urms**2.0/epsilon       !in s
182       Lf = 0.5 * (urms**3.0)/epsilon !in cm
183
184       tauk = (anu/epsilon)**0.5       !in s
185       eta  = (anu**3.0/epsilon)**0.25  !in cm
186       vk   = eta/tauk                   !in cm/s
187
188       ao = (11.+7.*Relamda)/(205.+Relamda)
189
190       lambda = urms * sqrt(15.*anu/epsilon)     !in cm
191
192       tt = sqrt(2.*Relamda/(15.**0.5)/ao) * tauk !in s
193
194       CALL fallg !gives winf in cm/s
195
196       DO i = pstart, pend
197          vST(i) = winf(i)           !in cm/s
198          tau(i) = vST(i)/gravity    !in s
199          St(i)  = tau(i)/tauk
200       ENDDO
201
202!***** TO CALCULATE wr ********************************
203!from Aayala 2008b, page 38f
204
205       z  = tt/Tl
206       be = sqrt(2.0)*lambda/Lf
207
208       bbb = sqrt(1.0-2.0*be**2.0)
209       d1 = (1.+bbb)/2.0/bbb
210       e1 = Lf*(1.0+bbb)/2.0   !in cm
211       d2 = (1.0-bbb)/2.0/bbb
212       e2 = Lf*(1.0-bbb)/2.0   !in cm
213
214       ccc = sqrt(1.0-2.0*z**2.0)
215       b1 = (1.+ccc)/2./ccc
216       c1 = Tl*(1.+ccc)/2.   !in s
217       b2 = (1.-ccc)/2./ccc
218       c2 = Tl*(1.-ccc)/2.   !in s
219
220       DO i = pstart, pend
221          v1 = vST(i)         !in cm/s
222          t1 = tau(i)         !in s
223
224          DO j = pstart,i
225             rrp = (particles(i)%radius + particles(j)%radius) * 100.0 !radius in cm
226             v2  = vST(j)         !in cm/s
227             t2  = tau(j)         !in s
228
229             v1xysq =  b1*d1*PHI(c1,e1,v1,t1) &
230                     - b1*d2*PHI(c1,e2,v1,t1) &
231                     - b2*d1*PHI(c2,e1,v1,t1) &
232                     + b2*d2*PHI(c2,e2,v1,t1)
233             v1xysq = v1xysq * urms**2.0/t1    !in cm**2/s**2
234
235             vrms1xy= sqrt(v1xysq)             !in cm/s
236
237             v2xysq =  b1*d1*PHI(c1,e1,v2,t2) &
238                     - b1*d2*PHI(c1,e2,v2,t2) &
239                     - b2*d1*PHI(c2,e1,v2,t2) &
240                     + b2*d2*PHI(c2,e2,v2,t2)
241             v2xysq = v2xysq * urms**2.0/t2    !in cm**2/s**2
242
243              vrms2xy= sqrt(v2xysq)             !in cm/s
244
245             IF(vST(i).ge.vST(j)) THEN
246                v1 = vST(i)
247                t1 = tau(i)
248                v2 = vST(j)
249                t2 = tau(j)
250             ELSE
251                v1 = vST(j)
252                t1 = tau(j)
253                v2 = vST(i)
254                t2 = tau(i)
255             ENDIF
256
257             v1v2xy =  b1*d1*ZHI(c1,e1,v1,t1,v2,t2) &
258                     - b1*d2*ZHI(c1,e2,v1,t1,v2,t2) &
259                     - b2*d1*ZHI(c2,e1,v1,t1,v2,t2) &
260                     + b2*d2*ZHI(c2,e2,v1,t1,v2,t2)
261             fR = d1 * exp(-rrp/e1) - d2 * exp(-rrp/e2)
262             v1v2xy = v1v2xy * fR * urms**2.0/tau(i)/tau(j)    !in cm**2/s**2
263
264             wrtur2xy=vrms1xy**2.0 + vrms2xy**2.0 - 2.*v1v2xy  !in cm**2/s**2
265
266              IF (wrtur2xy.le.0.0) wrtur2xy=0.0
267
268              wrgrav2=pi/8.*(vST(j)-vST(i))**2.0
269
270              wrFIN = sqrt((2.0/pi)*(wrtur2xy+wrgrav2))    !in cm/s
271
272
273!***** TO CALCULATE gr ********************************
274
275             IF(St(j).gt.St(i)) THEN
276                SSt = St(j)
277             ELSE
278                SSt = St(i)
279             ENDIF
280
281             XX = -0.1988*SSt**4.0 + 1.5275*SSt**3.0 &
282                  -4.2942*SSt**2.0 + 5.3406*SSt
283
284             IF(XX.le.0.0) XX = 0.0
285
286             YY = 0.1886*exp(20.306/Relamda)
287
288             c1_gr = XX/(gravity/(vk/tauk))**YY
289
290             ao_gr  = ao + (pi/8.)*(gravity/(vk/tauk))**2.0
291             fao_gr = 20.115 * (ao_gr/Relamda)**0.5
292             rc     = sqrt( fao_gr * abs(St(j)-St(i)) ) * eta   !in cm
293
294             grFIN = ((eta**2.0+rc**2.0)/(rrp**2.0+rc**2.0))**(c1_gr/2.)
295             IF (grFIN.lt.1.0) grFIN = 1.0
296
297             gck(i,j) = 2. * pi * rrp**2.0 * wrFIN * grFIN       ! in cm**3/s
298             gck(j,i) = gck(i,j)
299
300          ENDDO
301       ENDDO
302
303
304       RETURN
305    END SUBROUTINE TurbSD
306
307!------------------------------------------------------------------------------!
308! PHI as a function
309!------------------------------------------------------------------------------!
310  REAL FUNCTION PHI(a,b,vsett,tau0)
311
312       IMPLICIT NONE
313
314       REAL :: a, aa1, b, vsett, tau0
315
316       aa1 = 1./tau0 + 1./a + vsett/b
317
318       PHI = 1./aa1 - vsett/2.0/b/aa1**2.0  !in s
319       RETURN
320    END FUNCTION PHI
321
322!------------------------------------------------------------------------------!
323! ZETA as a function
324!------------------------------------------------------------------------------!
325  REAL FUNCTION ZHI(a,b,vsett1,tau1,vsett2,tau2)
326
327       IMPLICIT NONE
328
329       REAL :: a, aa1, aa2, aa3, aa4, aa5, aa6, b, vsett1, tau1, vsett2, tau2
330
331        aa1 = vsett2/b - 1./tau2 - 1./a
332        aa2 = vsett1/b + 1./tau1 + 1./a
333        aa3 = (vsett1-vsett2)/b + 1./tau1 + 1./tau2
334        aa4 = (vsett2/b)**2.0 - (1./tau2 + 1./a)**2.0
335        aa5 = vsett2/b + 1./tau2 + 1./a
336        aa6 = 1./tau1 - 1./a + (1./tau2 + 1./a) * vsett1/vsett2
337        ZHI =  (1./aa1 - 1./aa2) * (vsett1-vsett2)/2./b/aa3**2.0          &
338             + (4./aa4 - 1./aa5**2.0 - 1./aa1**2.0) * vsett2/2./b/aa6      &
339             + (2.*b/aa2 - 2.*b/aa1 - vsett1/aa2**2.0 + vsett2/aa1**2.0)   &
340                                                           * 1./2./b/aa3             ! in s**2
341        RETURN
342    END FUNCTION ZHI
343
344
345!------------------------------------------------------------------------------!
346! SUBROUTINE for calculation of terminal velocity winf
347!------------------------------------------------------------------------------!
348
349   SUBROUTINE fallg
350
351       USE constants
352       USE cloud_parameters
353       USE particle_attributes
354       USE arrays_3d
355
356      IMPLICIT NONE
357
358      INTEGER :: i, j
359
360      REAL :: eta, xlamb, rhoa, rhow, grav, cunh, t0, sigma, stok, stb, phy, py,  &
361              x, y, xrey, bond
362
363      REAL, DIMENSION(1:7)  :: b
364      REAL, DIMENSION(1:6)  :: c
365      REAL, DIMENSION(1:20) :: rat
366      REAL, DIMENSION(1:15) :: r0
367      REAL, DIMENSION(1:15,1:20) :: ecoll
368
369      b = (/-0.318657e1,0.992696e0,-0.153193e-2,-0.987059e-3, &
370             -0.578878e-3,0.855176e-4,-0.327815e-5/)
371      c = (/-0.500015e1,0.523778e1,-0.204914e1,0.475294e0,    &
372             -0.542819e-1, 0.238449e-2/)
373
374      eta   = 1.818e-4     !in poise = g/(cm s)
375      xlamb = 6.62e-6    !in cm
376
377      rhow = 1.0       !in g/cm**3
378      rhoa = 1.225e-3    !in g/cm**3
379
380      grav = 980.665     !in cm/s**2
381      cunh = 1.257 * xlamb !in cm
382      t0   = 273.15        !in K
383      sigma= 76.1 - 0.155 * (293.15 - t0) !in N/m = g/s**2
384      stok = 2.0 * grav * (rhow - rhoa)/(9.0 * eta) ! in 1/(cm s)
385      stb  = 32.0 * rhoa * (rhow-rhoa) * grav/(3.0 * eta * eta) ! in 1/cm**3
386      phy  = (sigma**3.0) * (rhoa**2.0)/((eta**4.0) * grav * (rhow-rhoa))
387      py   = phy**(1.0/6.0)
388
389!particle radius has to be in cm
390      DO j = pstart, pend
391
392         IF (particles(j)%radius*100.0 .le. 1.e-3) THEN
393
394            winf(j)=stok*((particles(j)%radius*100.0)**2.+cunh* particles(j)%radius*100.0) !in cm/s
395
396         ELSEIF (particles(j)%radius*100.0.gt.1.e-3.and.particles(j)%radius*100.0.le.5.35e-2) THEN
397
398            x = log(stb*(particles(j)%radius*100.0)**3.0)
399            y = 0.0
400
401            DO i = 1, 7
402               y = y + b(i) * (x**(i-1))
403            ENDDO
404
405            xrey = (1.0 + cunh/(particles(j)%radius*100.0)) * exp(y)
406            winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0)  !in cm/s
407
408         ELSEIF (particles(j)%radius*100.0.gt.5.35e-2) THEN
409
410            bond = grav*(rhow-rhoa)*(particles(j)%radius*100.0)**2.0/sigma
411
412            IF (particles(j)%radius*100.0.gt.0.35)  THEN
413               bond = grav*(rhow-rhoa) * 0.35 * 0.35/sigma
414            ENDIF
415
416            x = log(16.0*bond*py/3.0)
417            y = 0.0
418
419            DO i = 1, 6
420               y = y + c(i) * (x**(i-1))
421            ENDDO
422
423            xrey = py*exp(y)
424            winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0)  !in cm/s
425
426            IF (particles(j)%radius*100.0 .gt.0.35)  THEN
427              winf(j) = xrey * eta/(2.0 * rhoa * 0.35) !in cm/s
428            ENDIF
429
430         ENDIF
431      ENDDO
432      RETURN
433   END SUBROUTINE fallg
434
435!------------------------------------------------------------------------------!
436! SUBROUTINE for calculation of collision efficencies
437!------------------------------------------------------------------------------!
438
439   SUBROUTINE effic
440
441       USE constants
442       USE cloud_parameters
443       USE particle_attributes
444       USE arrays_3d
445
446!collision efficiencies of hall kernel
447      IMPLICIT NONE
448
449      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)
513               ELSE
514                  qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
515                  ec(j,i) = (1.-qq) * ecoll(1,iq-1) + qq * ecoll(1,iq)
516               ENDIF
517            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)
521            ENDIF
522            ec(i,j) = ec(j,i)
523            IF (ec(i,j).lt.1.e-20) stop 99
524         ENDDO
525      ENDDO
526      RETURN
527   END SUBROUTINE effic
528
529!------------------------------------------------------------------------------!
530! SUBROUTINE for calculation of enhancement factor collision efficencies
531!------------------------------------------------------------------------------!
532
533   SUBROUTINE turb_enhance_eff
534
535       USE constants
536       USE cloud_parameters
537       USE particle_attributes
538       USE arrays_3d
539
540!collision efficiencies of hall kernel
541      IMPLICIT NONE
542
543      INTEGER :: i, ik, ir, iq, j, k, kk
544      REAL    :: rq, y1, pp, qq, y2, y3, x1, x2, x3
545
546      REAL, DIMENSION(1:11) :: rat
547      REAL, DIMENSION(1:7)  :: r0
548      REAL, DIMENSION(1:7,1:11) :: ecoll_100, ecoll_400
549
550      r0  = (/10., 20., 30.,40., 50., 60.,100./)
551      rat = (/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0/)
552
553! 100 cm^2/s^3
554      ecoll_100(:,1) = (/1.74,  1.74,  1.773, 1.49,  1.207,  1.207,  1.0 /)
555      ecoll_100(:,2) = (/1.46,  1.46,  1.421, 1.245, 1.069,  1.069,  1.0 /)
556      ecoll_100(:,3) = (/1.32,  1.32,  1.245, 1.123, 1.000,  1.000,  1.0 /)
557      ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025,  1.025,  1.0 /)
558      ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056,  1.056,  1.0 /)
559      ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028,  1.028,  1.0 /)
560      ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046,  1.046,  1.0 /)
561      ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029,  1.029,  1.0 /)
562      ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021,  1.021,  1.0 /)
563      ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088,  1.088,  1.0 /)
564      ecoll_100(:,11)= (/20.3,  20.3,  14.6 , 8.61,  2.60,   2.60 ,  1.0 /)
565
566! 400 cm^2/s^3
567      ecoll_400(:,1) = (/4.976, 4.976,  3.593, 2.519, 1.445,  1.445,  1.0 /)
568      ecoll_400(:,2) = (/2.984, 2.984,  2.181, 1.691, 1.201,  1.201,  1.0 /)
569      ecoll_400(:,3) = (/1.988, 1.988,  1.475, 1.313, 1.150,  1.150,  1.0 /)
570      ecoll_400(:,4) = (/1.490, 1.490,  1.187, 1.156, 1.126,  1.126,  1.0 /)
571      ecoll_400(:,5) = (/1.249, 1.249,  1.088, 1.090, 1.092,  1.092,  1.0 /)
572      ecoll_400(:,6) = (/1.139, 1.139,  1.130, 1.091, 1.051,  1.051,  1.0 /)
573      ecoll_400(:,7) = (/1.220, 1.220,  1.190, 1.138, 1.086,  1.086,  1.0 /)
574      ecoll_400(:,8) = (/1.325, 1.325,  1.267, 1.165, 1.063,  1.063,  1.0 /)
575      ecoll_400(:,9) = (/1.716, 1.716,  1.345, 1.223, 1.100,  1.100,  1.0 /)
576      ecoll_400(:,10)= (/3.788, 3.788,  1.501, 1.311, 1.120,  1.120,  1.0 /)
577      ecoll_400(:,11)= (/36.52, 36.52,  19.16, 22.80,  26.0,   26.0,  1.0 /)
578
579! two-dimensional linear interpolation of the collision efficiency
580      DO j =  pstart, pend
581         DO i = pstart, j
582
583            DO k = 2, 7
584               IF ((particles(j)%radius*1.0E06).le.r0(k).and.(particles(j)%radius*1.0E06).ge.r0(k-1)) THEN
585                  ir = k
586               ELSEIF ((particles(j)%radius*1.0E06).gt.r0(7)) THEN
587                  ir = 8
588               ELSEIF ((particles(j)%radius*1.0E06).lt.r0(1)) THEN
589                  ir = 1
590               ENDIF
591            ENDDO
592
593            rq = particles(i)%radius*1.0E06/(particles(j)%radius*1.0E06)
594
595            DO kk = 2, 11
596               IF (rq.le.rat(kk).and.rq.gt.rat(kk-1)) iq = kk
597            ENDDO
598
599! 0  cm2/s3
600            y1 = 1.0
601! 100 cm2/s3, 400 cm2/s3
602            IF (ir.lt.8) THEN
603               IF (ir.ge.2) THEN
604                  pp = ((particles(j)%radius*1.0E06)-r0(ir-1))/(r0(ir)-r0(ir-1))
605                  qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
606                  y2= (1.-pp)*(1.-qq)*ecoll_100(ir-1,iq-1)+  &
607                            pp*(1.-qq)*ecoll_100(ir,iq-1)+    &
608                            qq*(1.-pp)*ecoll_100(ir-1,iq)+    &
609                            pp*qq*ecoll_100(ir,iq)
610                  y3= (1.-pp)*(1.-qq)*ecoll_400(ir-1,iq-1)+  &
611                            pp*(1.-qq)*ecoll_400(ir,iq-1)+    &
612                            qq*(1.-pp)*ecoll_400(ir-1,iq)+    &
613                            pp*qq*ecoll_400(ir,iq)
614               ELSE
615                  qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
616                  y2= (1.-qq)*ecoll_100(1,iq-1)+qq*ecoll_100(1,iq)
617                  y3= (1.-qq)*ecoll_400(1,iq-1)+qq*ecoll_400(1,iq)
618               ENDIF
619            ELSE
620            qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
621            y2 = (1.-qq) * ecoll_100(7,iq-1) + qq * ecoll_100(7,iq)
622            y3 = (1.-qq) * ecoll_400(7,iq-1) + qq * ecoll_400(7,iq)
623         ENDIF
624! linear interpolation
625! dissipation rate in cm2/s3
626         x1 = 0.0
627         x2 = 100.0
628         x3 = 400.0
629
630         IF (epsilon.le.100.) THEN
631            ecf(j,i) = (epsilon-100.)/(0.-100.) * y1 &
632                    +  (epsilon-0.)/(100.-0.) * y2
633         ELSE IF(epsilon.le.600.)THEN
634            ecf(j,i) = (epsilon-400.)/(100.-400.) * y2 &
635                    +  (epsilon-100.)/(400.-100.) * y3
636
637         ELSE
638            ecf(j,i) = (600.-400.)/(100.-400.) * y2 &
639                    +  (600.-100.)/(400.-100.) * y3
640         ENDIF
641
642         IF (ecf(j,i).lt.1.0) ecf(j,i) = 1.0
643
644         ecf(i,j)=ecf(j,i)
645      ENDDO
646   ENDDO
647
648   RETURN
649   END SUBROUTINE turb_enhance_eff
650
651 END MODULE wang_kernel_mod
Note: See TracBrowser for help on using the repository browser.