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

Last change on this file since 799 was 799, checked in by franke, 10 years ago

Implementation of Wang collision kernel and bugfix for calculation of vpt, pt_p, and ec in case of cloud droplets

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