1 | MODULE 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 |
---|