source: palm/trunk/SOURCE/lpm_collision_kernels.f90 @ 826

Last change on this file since 826 was 826, checked in by raasch, 10 years ago

last commit documented

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