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

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

further adjustments for speedup of particle code

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