Home

Context Navigation

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

Last change on this file since 936 was 850, checked in by raasch, 13 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 35.4 KB

Line
1	MODULE lpm_collision_kernels_mod
2
3	!------------------------------------------------------------------------------!
4	! Current revisions:
5	! -----------------
6	!
7	!
8	! Former revisions:
9	! -----------------
10	! $Id: lpm_collision_kernels.f90 850 2012-03-15 12:09:25Z raasch $
11	!
12	! 849 2012-03-15 10:35:09Z raasch
13	! routine collision_efficiency_rogers added (moved from former advec_particles
14	! to here)
15	!
16	! 835 2012-02-22 11:21:19Z raasch $
17	! Bugfix: array diss can be used only in case of Wang kernel
18	!
19	! 828 2012-02-21 12:00:36Z raasch
20	! code has been completely reformatted, routine colker renamed
21	! recalculate_kernel,
22	! routine init_kernels added, radius is now communicated to the collision
23	! routines by array radclass
24	!
25	! Bugfix: transformation factor for dissipation changed from 1E5 to 1E4
26	!
27	! 825 2012-02-19 03:03:44Z raasch
28	! routine renamed from wang_kernel to lpm_collision_kernels,
29	! turbulence_effects on collision replaced by wang_kernel
30	!
31	! 799 2011-12-21 17:48:03Z franke
32	! speed optimizations and formatting
33	! Bugfix: iq=1 is not allowed (routine effic)
34	! Bugfix: replaced stop by ec=0.0 in case of very small ec (routine effic)
35	!
36	! 790 2011-11-29 03:11:20Z raasch
37	! initial revision
38	!
39	! Description:
40	! ------------
41	! This module calculates collision efficiencies either due to pure gravitational
42	! effects (Hall kernel, see Hall, 1980: J. Atmos. Sci., 2486-2507) or
43	! including the effects of (SGS) turbulence (Wang kernel, see Wang and
44	! Grabowski, 2009: Atmos. Sci. Lett., 10, 1-8). The original code has been
45	! provided by L.-P. Wang but is substantially reformatted and speed optimized
46	! here.
47	!
48	! ATTENTION:
49	! Physical quantities (like g, densities, etc.) used in this module still
50	! have to be adjusted to those values used in the main PALM code.
51	! Also, quantities in CGS-units should be converted to SI-units eventually.
52	!------------------------------------------------------------------------------!
53
54	USE arrays_3d
55	USE cloud_parameters
56	USE constants
57	USE particle_attributes
58	USE pegrid
59
60
61	IMPLICIT NONE
62
63	PRIVATE
64
65	PUBLIC ckernel, collision_efficiency_rogers, init_kernels, &
66	rclass_lbound, rclass_ubound, recalculate_kernel
67
68	REAL :: epsilon, eps2, rclass_lbound, rclass_ubound, urms, urms2
69
70	REAL, DIMENSION(:), ALLOCATABLE :: epsclass, radclass, winf
71	REAL, DIMENSION(:,:), ALLOCATABLE :: ec, ecf, gck, hkernel, hwratio
72	REAL, DIMENSION(:,:,:), ALLOCATABLE :: ckernel
73
74	SAVE
75
76	!
77	!-- Public interfaces
78	INTERFACE collision_efficiency_rogers
79	MODULE PROCEDURE collision_efficiency_rogers
80	END INTERFACE collision_efficiency_rogers
81
82	INTERFACE init_kernels
83	MODULE PROCEDURE init_kernels
84	END INTERFACE init_kernels
85
86	INTERFACE recalculate_kernel
87	MODULE PROCEDURE recalculate_kernel
88	END INTERFACE recalculate_kernel
89
90
91	CONTAINS
92
93
94	SUBROUTINE init_kernels
95	!------------------------------------------------------------------------------!
96	! Initialization of the collision efficiency matrix with fixed radius and
97	! dissipation classes, calculated at simulation start only.
98	!------------------------------------------------------------------------------!
99
100	IMPLICIT NONE
101
102	INTEGER :: i, j, k
103
104
105	!
106	!-- Calculate collision efficiencies for fixed radius- and dissipation
107	!-- classes
108	IF ( collision_kernel(6:9) == 'fast' ) THEN
109
110	ALLOCATE( ckernel(1:radius_classes,1:radius_classes, &
111	0:dissipation_classes), epsclass(1:dissipation_classes), &
112	radclass(1:radius_classes) )
113
114	!
115	!-- Calculate the radius class bounds with logarithmic distances
116	!-- in the interval [1.0E-6, 2.0E-4] m
117	rclass_lbound = LOG( 1.0E-6 )
118	rclass_ubound = LOG( 2.0E-4 )
119	radclass(1) = 1.0E-6
120	DO i = 2, radius_classes
121	radclass(i) = EXP( rclass_lbound + &
122	( rclass_ubound - rclass_lbound ) * ( i-1.0 ) /&
123	( radius_classes - 1.0 ) )
124	! IF ( myid == 0 ) THEN
125	! PRINT, 'i=', i, ' r = ', radclass(i)1.0E6
126	! ENDIF
127	ENDDO
128	!
129	!-- Collision routines expect radius to be in cm
130	radclass = radclass * 100.0
131
132	!
133	!-- Set the class bounds for dissipation in interval [0, 1000] cm2/s3
134	DO i = 1, dissipation_classes
135	epsclass(i) = 1000.0 * REAL( i ) / dissipation_classes
136	! IF ( myid == 0 ) THEN
137	! PRINT*, 'i=', i, ' eps = ', epsclass(i)
138	! ENDIF
139	ENDDO
140	!
141	!-- Calculate collision efficiencies of the Wang/ayala kernel
142	ALLOCATE( ec(1:radius_classes,1:radius_classes), &
143	ecf(1:radius_classes,1:radius_classes), &
144	gck(1:radius_classes,1:radius_classes), &
145	winf(1:radius_classes) )
146
147	DO k = 1, dissipation_classes
148
149	epsilon = epsclass(k)
150	urms = 202.0 * ( epsilon / 400.0 )**( 1.0 / 3.0 )
151
152	CALL turbsd
153	CALL turb_enhance_eff
154	CALL effic
155
156	DO j = 1, radius_classes
157	DO i = 1, radius_classes
158	ckernel(i,j,k) = ec(i,j) * gck(i,j) * ecf(i,j)
159	ENDDO
160	ENDDO
161
162	ENDDO
163
164	!
165	!-- Calculate collision efficiencies of the Hall kernel
166	ALLOCATE( hkernel(1:radius_classes,1:radius_classes), &
167	hwratio(1:radius_classes,1:radius_classes) )
168
169	CALL fallg
170	CALL effic
171
172	DO j = 1, radius_classes
173	DO i = 1, radius_classes
174	hkernel(i,j) = pi * ( radclass(j) + radclass(i) )**2 &
175	* ec(i,j) * ABS( winf(j) - winf(i) )
176	ckernel(i,j,0) = hkernel(i,j) ! hall kernel stored on index 0
177	ENDDO
178	ENDDO
179
180	!
181	!-- Test output of efficiencies
182	IF ( j == -1 ) THEN
183
184	PRINT, '** Hall kernel'
185	WRITE ( ,'(5X,20(F4.0,1X))' ) ( radclass(i)1.0E4, i = 1,radius_classes )
186	DO j = 1, radius_classes
187	WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j), ( hkernel(i,j), i = 1,radius_classes )
188	ENDDO
189
190	DO k = 1, dissipation_classes
191	DO i = 1, radius_classes
192	DO j = 1, radius_classes
193	IF ( hkernel(i,j) == 0.0 ) THEN
194	hwratio(i,j) = 9999999.9
195	ELSE
196	hwratio(i,j) = ckernel(i,j,k) / hkernel(i,j)
197	ENDIF
198	ENDDO
199	ENDDO
200
201	PRINT, '** epsilon = ', epsclass(k)
202	WRITE ( ,'(5X,20(F4.0,1X))' ) ( radclass(i)1.0E4, i = 1,radius_classes )
203	DO j = 1, radius_classes
204	! WRITE ( ,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)1.0E4, ( ckernel(i,j,k), i = 1,radius_classes )
205	WRITE ( ,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)1.0E4, ( hwratio(i,j), i = 1,radius_classes )
206	ENDDO
207	ENDDO
208
209	ENDIF
210
211	DEALLOCATE( ec, ecf, epsclass, gck, hkernel, winf )
212
213	ckernel = ckernel * 1.0E-6 ! kernel is needed in m**3/s
214
215	ELSEIF( collision_kernel == 'hall' .OR. collision_kernel == 'wang' ) &
216	THEN
217	!
218	!-- Initial settings for Hall- and Wang-Kernel
219	!-- To be done: move here parts from turbsd, fallg, ecoll, etc.
220	ENDIF
221
222	END SUBROUTINE init_kernels
223
224
225	!------------------------------------------------------------------------------!
226	! Calculation of collision kernels during each timestep and for each grid box
227	!------------------------------------------------------------------------------!
228	SUBROUTINE recalculate_kernel( i1, j1, k1 )
229
230	USE arrays_3d
231	USE cloud_parameters
232	USE constants
233	USE cpulog
234	USE indices
235	USE interfaces
236	USE particle_attributes
237
238	IMPLICIT NONE
239
240	INTEGER :: i, i1, j, j1, k1, pend, pstart
241
242
243	pstart = prt_start_index(k1,j1,i1)
244	pend = prt_start_index(k1,j1,i1) + prt_count(k1,j1,i1) - 1
245	radius_classes = prt_count(k1,j1,i1)
246
247	ALLOCATE( ec(1:radius_classes,1:radius_classes), &
248	radclass(1:radius_classes), winf(1:radius_classes) )
249
250	!
251	!-- Store particle radii on the radclass array. Collision routines
252	!-- expect radii to be in cm.
253	radclass(1:radius_classes) = particles(pstart:pend)%radius * 100.0
254
255	IF ( wang_kernel ) THEN
256	epsilon = diss(k1,j1,i1) * 1.0E4 ! dissipation rate in cm2/s-3
257	ELSE
258	epsilon = 0.0
259	ENDIF
260	urms = 202.0 * ( epsilon / 400.0 )**( 0.33333333333 )
261
262	IF ( wang_kernel .AND. epsilon > 0.001 ) THEN
263	!
264	!-- Call routines to calculate efficiencies for the Wang kernel
265	ALLOCATE( gck(1:radius_classes,1:radius_classes), &
266	ecf(1:radius_classes,1:radius_classes) )
267
268	CALL turbsd
269	CALL turb_enhance_eff
270	CALL effic
271
272	DO j = 1, radius_classes
273	DO i = 1, radius_classes
274	ckernel(pstart+i-1,pstart+j-1,1) = ec(i,j) * gck(i,j) * ecf(i,j)
275	ENDDO
276	ENDDO
277
278	DEALLOCATE( gck, ecf )
279
280	ELSE
281	!
282	!-- Call routines to calculate efficiencies for the Hall kernel
283	CALL fallg
284	CALL effic
285
286	DO j = 1, radius_classes
287	DO i = 1, radius_classes
288	ckernel(pstart+i-1,pstart+j-1,1) = pi * &
289	( radclass(j) + radclass(i) )**2 &
290	* ec(i,j) * ABS( winf(j) - winf(i) )
291	ENDDO
292	ENDDO
293
294	ENDIF
295
296	ckernel = ckernel * 1.0E-6 ! kernel is needed in m**3/s
297
298	DEALLOCATE( ec, radclass, winf )
299
300	END SUBROUTINE recalculate_kernel
301
302
303	!------------------------------------------------------------------------------!
304	! Calculation of gck
305	! This is from Aayala 2008b, page 37ff.
306	! Necessary input parameters: water density, radii of droplets, air density,
307	! air viscosity, turbulent dissipation rate, taylor microscale reynolds number,
308	! gravitational acceleration --> to be replaced by PALM parameters
309	!------------------------------------------------------------------------------!
310	SUBROUTINE turbsd
311
312	USE constants
313	USE cloud_parameters
314	USE particle_attributes
315	USE arrays_3d
316
317	IMPLICIT NONE
318
319	INTEGER :: i, j
320
321	LOGICAL, SAVE :: first = .TRUE.
322
323	REAL :: ao, ao_gr, bbb, be, b1, b2, ccc, c1, c1_gr, c2, d1, d2, eta, &
324	e1, e2, fao_gr, fr, grfin, lambda, lambda_re, lf, rc, rrp, &
325	sst, tauk, tl, t2, tt, t1, vk, vrms1xy, vrms2xy, v1, v1v2xy, &
326	v1xysq, v2, v2xysq, wrfin, wrgrav2, wrtur2xy, xx, yy, z
327
328	REAL, SAVE :: airdens, airvisc, anu, gravity, waterdens
329
330	REAL, DIMENSION(1:radius_classes) :: st, tau
331
332
333	!
334	!-- Initial assignment of constants
335	IF ( first ) THEN
336
337	first = .FALSE.
338	airvisc = 0.1818 ! dynamic viscosity in mg/cm*s
339	airdens = 1.2250 ! air density in mg/cm**3
340	waterdens = 1000.0 ! water density in mg/cm**3
341	gravity = 980.6650 ! in cm/s**2
342	anu = airvisc/airdens ! kinetic viscosity in cm**2/s
343
344	ENDIF
345
346	lambda = urms * SQRT( 15.0 * anu / epsilon ) ! in cm
347	lambda_re = urms*2 SQRT( 15.0 / epsilon / anu )
348	tl = urms**2 / epsilon ! in s
349	lf = 0.5 * urms**3 / epsilon ! in cm
350	tauk = SQRT( anu / epsilon ) ! in s
351	eta = ( anu3 / epsilon )0.25 ! in cm
352	vk = eta / tauk ! in cm/s
353
354	ao = ( 11.0 + 7.0 * lambda_re ) / ( 205.0 + lambda_re )
355	tt = SQRT( 2.0 * lambda_re / ( SQRT( 15.0 ) * ao ) ) * tauk ! in s
356
357	CALL fallg ! gives winf in cm/s
358
359	DO i = 1, radius_classes
360	tau(i) = winf(i) / gravity ! in s
361	st(i) = tau(i) / tauk
362	ENDDO
363
364	!
365	!-- Calculate wr (from Aayala 2008b, page 38f)
366	z = tt / tl
367	be = SQRT( 2.0 ) * lambda / lf
368	bbb = SQRT( 1.0 - 2.0 * be**2 )
369	d1 = ( 1.0 + bbb ) / ( 2.0 * bbb )
370	e1 = lf * ( 1.0 + bbb ) * 0.5 ! in cm
371	d2 = ( 1.0 - bbb ) * 0.5 / bbb
372	e2 = lf * ( 1.0 - bbb ) * 0.5 ! in cm
373	ccc = SQRT( 1.0 - 2.0 * z**2 )
374	b1 = ( 1.0 + ccc ) * 0.5 / ccc
375	c1 = tl * ( 1.0 + ccc ) * 0.5 ! in s
376	b2 = ( 1.0 - ccc ) * 0.5 / ccc
377	c2 = tl * ( 1.0 - ccc ) * 0.5 ! in s
378
379	DO i = 1, radius_classes
380
381	v1 = winf(i) ! in cm/s
382	t1 = tau(i) ! in s
383
384	DO j = 1, i
385	rrp = radclass(i) + radclass(j) ! radius in cm
386	v2 = winf(j) ! in cm/s
387	t2 = tau(j) ! in s
388
389	v1xysq = b1 * d1 * phi(c1,e1,v1,t1) - b1 * d2 * phi(c1,e2,v1,t1) &
390	- b2 * d1 * phi(c2,e1,v1,t1) + b2 * d2 * phi(c2,e2,v1,t1)
391	v1xysq = v1xysq * urms2 / t1 ! in cm2/s**2
392	vrms1xy = SQRT( v1xysq ) ! in cm/s
393
394	v2xysq = b1 * d1 * phi(c1,e1,v2,t2) - b1 * d2 * phi(c1,e2,v2,t2) &
395	- b2 * d1 * phi(c2,e1,v2,t2) + b2 * d2 * phi(c2,e2,v2,t2)
396	v2xysq = v2xysq * urms2 / t2 ! in cm2/s**2
397	vrms2xy = SQRT( v2xysq ) ! in cm/s
398
399	IF ( winf(i) >= winf(j) ) THEN
400	v1 = winf(i)
401	t1 = tau(i)
402	v2 = winf(j)
403	t2 = tau(j)
404	ELSE
405	v1 = winf(j)
406	t1 = tau(j)
407	v2 = winf(i)
408	t2 = tau(i)
409	ENDIF
410
411	v1v2xy = b1 * d1 * zhi(c1,e1,v1,t1,v2,t2) - &
412	b1 * d2 * zhi(c1,e2,v1,t1,v2,t2) - &
413	b2 * d1 * zhi(c2,e1,v1,t1,v2,t2) + &
414	b2 * d2* zhi(c2,e2,v1,t1,v2,t2)
415	fr = d1 * EXP( -rrp / e1 ) - d2 * EXP( -rrp / e2 )
416	v1v2xy = v1v2xy * fr * urms2 / tau(i) / tau(j) ! in cm2/s**2
417	wrtur2xy = vrms1xy2 + vrms2xy2 - 2.0 * v1v2xy ! in cm2/s2
418	IF ( wrtur2xy < 0.0 ) wrtur2xy = 0.0
419	wrgrav2 = pi / 8.0 * ( winf(j) - winf(i) )**2
420	wrfin = SQRT( ( 2.0 / pi ) * ( wrtur2xy + wrgrav2) ) ! in cm/s
421
422	!
423	!-- Calculate gr
424	IF ( st(j) > st(i) ) THEN
425	sst = st(j)
426	ELSE
427	sst = st(i)
428	ENDIF
429
430	xx = -0.1988 * sst*4 + 1.5275 sst*3 - 4.2942 sst**2 + &
431	5.3406 * sst
432	IF ( xx < 0.0 ) xx = 0.0
433	yy = 0.1886 * EXP( 20.306 / lambda_re )
434
435	c1_gr = xx / ( gravity / vk * tauk )**yy
436
437	ao_gr = ao + ( pi / 8.0) * ( gravity / vk * tauk )**2
438	fao_gr = 20.115 * SQRT( ao_gr / lambda_re )
439	rc = SQRT( fao_gr * ABS( st(j) - st(i) ) ) * eta ! in cm
440
441	grfin = ( ( eta2 + rc2 ) / ( rrp2 + rc2) )*( c1_gr0.5 )
442	IF ( grfin < 1.0 ) grfin = 1.0
443
444	gck(i,j) = 2.0 * pi * rrp*2 wrfin * grfin ! in cm**3/s
445	gck(j,i) = gck(i,j)
446
447	ENDDO
448	ENDDO
449
450	END SUBROUTINE turbsd
451
452
453	!------------------------------------------------------------------------------!
454	! phi as a function
455	!------------------------------------------------------------------------------!
456	REAL FUNCTION phi( a, b, vsett, tau0 )
457
458	IMPLICIT NONE
459
460	REAL :: a, aa1, b, tau0, vsett
461
462	aa1 = 1.0 / tau0 + 1.0 / a + vsett / b
463	phi = 1.0 / aa1 - 0.5 * vsett / b / aa1**2 ! in s
464
465	END FUNCTION phi
466
467
468	!------------------------------------------------------------------------------!
469	! zeta as a function
470	!------------------------------------------------------------------------------!
471	REAL FUNCTION zhi( a, b, vsett1, tau1, vsett2, tau2 )
472
473	IMPLICIT NONE
474
475	REAL :: a, aa1, aa2, aa3, aa4, aa5, aa6, b, tau1, tau2, vsett1, vsett2
476
477	aa1 = vsett2 / b - 1.0 / tau2 - 1.0 / a
478	aa2 = vsett1 / b + 1.0 / tau1 + 1.0 / a
479	aa3 = ( vsett1 - vsett2 ) / b + 1.0 / tau1 + 1.0 / tau2
480	aa4 = ( vsett2 / b )2 - ( 1.0 / tau2 + 1.0 / a )2
481	aa5 = vsett2 / b + 1.0 / tau2 + 1.0 / a
482	aa6 = 1.0 / tau1 - 1.0 / a + ( 1.0 / tau2 + 1.0 / a) * vsett1 / vsett2
483	zhi = (1.0 / aa1 - 1.0 / aa2 ) * ( vsett1 - vsett2 ) * 0.5 / b / aa3**2 &
484	+ (4.0 / aa4 - 1.0 / aa52 - 1.0 / aa12 ) * vsett2 * 0.5 / b /aa6&
485	+ (2.0 * ( b / aa2 - b / aa1 ) - vsett1 / aa22 + vsett2 / aa12 )&
486	* 0.5 / b / aa3 ! in s**2
487
488	END FUNCTION zhi
489
490
491	!------------------------------------------------------------------------------!
492	! Calculation of terminal velocity winf
493	!------------------------------------------------------------------------------!
494	SUBROUTINE fallg
495
496	USE constants
497	USE cloud_parameters
498	USE particle_attributes
499	USE arrays_3d
500
501	IMPLICIT NONE
502
503	INTEGER :: i, j
504
505	LOGICAL, SAVE :: first = .TRUE.
506
507	REAL, SAVE :: cunh, eta, grav, phy, py, rhoa, rhow, sigma, stb, stok, &
508	t0, xlamb
509
510	REAL :: bond, x, xrey, y
511
512	REAL, DIMENSION(1:7), SAVE :: b
513	REAL, DIMENSION(1:6), SAVE :: c
514
515	!
516	!-- Initial assignment of constants
517	IF ( first ) THEN
518
519	first = .FALSE.
520	b = (/ -0.318657E1, 0.992696E0, -0.153193E-2, -0.987059E-3, &
521	-0.578878E-3, 0.855176E-4, -0.327815E-5 /)
522	c = (/ -0.500015E1, 0.523778E1, -0.204914E1, 0.475294E0, &
523	-0.542819E-1, 0.238449E-2 /)
524
525	eta = 1.818E-4 ! in poise = g/(cm s)
526	xlamb = 6.62E-6 ! in cm
527	rhow = 1.0 ! in g/cm**3
528	rhoa = 1.225E-3 ! in g/cm**3
529	grav = 980.665 ! in cm/s**2
530	cunh = 1.257 * xlamb ! in cm
531	t0 = 273.15 ! in K
532	sigma = 76.1 - 0.155 * ( 293.15 - t0 ) ! in N/m = g/s**2
533	stok = 2.0 * grav * ( rhow - rhoa ) / ( 9.0 * eta ) ! in 1/(cm s)
534	stb = 32.0 * rhoa * ( rhow - rhoa) * grav / (3.0 * eta * eta)
535	! in 1/cm**3
536	phy = sigma*3 rhoa2 / ( eta4 * grav * ( rhow - rhoa ) )
537	py = phy**( 1.0 / 6.0 )
538
539	ENDIF
540
541	DO j = 1, radius_classes
542
543	IF ( radclass(j) <= 1.0E-3 ) THEN
544
545	winf(j) = stok * ( radclass(j)*2 + cunh radclass(j) ) ! in cm/s
546
547	ELSEIF ( radclass(j) > 1.0E-3 .AND. radclass(j) <= 5.35E-2 ) THEN
548
549	x = LOG( stb * radclass(j)**3 )
550	y = 0.0
551
552	DO i = 1, 7
553	y = y + b(i) * x**(i-1)
554	ENDDO
555
556	xrey = ( 1.0 + cunh / radclass(j) ) * EXP( y )
557	winf(j) = xrey * eta / ( 2.0 * rhoa * radclass(j) ) ! in cm/s
558
559	ELSEIF ( radclass(j) > 5.35E-2 ) THEN
560
561	IF ( radclass(j) > 0.35 ) THEN
562	bond = grav * ( rhow - rhoa ) * 0.35**2 / sigma
563	ELSE
564	bond = grav * ( rhow - rhoa ) * radclass(j)**2 / sigma
565	ENDIF
566
567	x = LOG( 16.0 * bond * py / 3.0 )
568	y = 0.0
569
570	DO i = 1, 6
571	y = y + c(i) * x**(i-1)
572	ENDDO
573
574	xrey = py * EXP( y )
575
576	IF ( radclass(j) > 0.35 ) THEN
577	winf(j) = xrey * eta / ( 2.0 * rhoa * 0.35 ) ! in cm/s
578	ELSE
579	winf(j) = xrey * eta / ( 2.0 * rhoa * radclass(j) ) ! in cm/s
580	ENDIF
581
582	ENDIF
583
584	ENDDO
585
586	END SUBROUTINE fallg
587
588
589	!------------------------------------------------------------------------------!
590	! Calculation of collision efficencies for the Hall kernel
591	!------------------------------------------------------------------------------!
592	SUBROUTINE effic
593
594	USE arrays_3d
595	USE cloud_parameters
596	USE constants
597	USE particle_attributes
598
599	IMPLICIT NONE
600
601	INTEGER :: i, iq, ir, j, k, kk
602
603	INTEGER, DIMENSION(:), ALLOCATABLE :: ira
604
605	LOGICAL, SAVE :: first = .TRUE.
606
607	REAL :: ek, particle_radius, pp, qq, rq
608
609	REAL, DIMENSION(1:21), SAVE :: rat
610	REAL, DIMENSION(1:15), SAVE :: r0
611	REAL, DIMENSION(1:15,1:21), SAVE :: ecoll
612
613	!
614	!-- Initial assignment of constants
615	IF ( first ) THEN
616
617	first = .FALSE.
618	r0 = (/ 6.0, 8.0, 10.0, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0, 60., &
619	70.0, 100.0, 150.0, 200.0, 300.0 /)
620	rat = (/ 0.00, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, &
621	0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, &
622	1.00 /)
623
624	ecoll(:,1) = (/0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, &
625	0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001/)
626	ecoll(:,2) = (/0.003, 0.003, 0.003, 0.004, 0.005, 0.005, 0.005, &
627	0.010, 0.100, 0.050, 0.200, 0.500, 0.770, 0.870, 0.970/)
628	ecoll(:,3) = (/0.007, 0.007, 0.007, 0.008, 0.009, 0.010, 0.010, &
629	0.070, 0.400, 0.430, 0.580, 0.790, 0.930, 0.960, 1.000/)
630	ecoll(:,4) = (/0.009, 0.009, 0.009, 0.012, 0.015, 0.010, 0.020, &
631	0.280, 0.600, 0.640, 0.750, 0.910, 0.970, 0.980, 1.000/)
632	ecoll(:,5) = (/0.014, 0.014, 0.014, 0.015, 0.016, 0.030, 0.060, &
633	0.500, 0.700, 0.770, 0.840, 0.950, 0.970, 1.000, 1.000/)
634	ecoll(:,6) = (/0.017, 0.017, 0.017, 0.020, 0.022, 0.060, 0.100, &
635	0.620, 0.780, 0.840, 0.880, 0.950, 1.000, 1.000, 1.000/)
636	ecoll(:,7) = (/0.030, 0.030, 0.024, 0.022, 0.032, 0.062, 0.200, &
637	0.680, 0.830, 0.870, 0.900, 0.950, 1.000, 1.000, 1.000/)
638	ecoll(:,8) = (/0.025, 0.025, 0.025, 0.036, 0.043, 0.130, 0.270, &
639	0.740, 0.860, 0.890, 0.920, 1.000, 1.000, 1.000, 1.000/)
640	ecoll(:,9) = (/0.027, 0.027, 0.027, 0.040, 0.052, 0.200, 0.400, &
641	0.780, 0.880, 0.900, 0.940, 1.000, 1.000, 1.000, 1.000/)
642	ecoll(:,10)= (/0.030, 0.030, 0.030, 0.047, 0.064, 0.250, 0.500, &
643	0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
644	ecoll(:,11)= (/0.040, 0.040, 0.033, 0.037, 0.068, 0.240, 0.550, &
645	0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
646	ecoll(:,12)= (/0.035, 0.035, 0.035, 0.055, 0.079, 0.290, 0.580, &
647	0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
648	ecoll(:,13)= (/0.037, 0.037, 0.037, 0.062, 0.082, 0.290, 0.590, &
649	0.780, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
650	ecoll(:,14)= (/0.037, 0.037, 0.037, 0.060, 0.080, 0.290, 0.580, &
651	0.770, 0.890, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
652	ecoll(:,15)= (/0.037, 0.037, 0.037, 0.041, 0.075, 0.250, 0.540, &
653	0.760, 0.880, 0.920, 0.950, 1.000, 1.000, 1.000, 1.000/)
654	ecoll(:,16)= (/0.037, 0.037, 0.037, 0.052, 0.067, 0.250, 0.510, &
655	0.770, 0.880, 0.930, 0.970, 1.000, 1.000, 1.000, 1.000/)
656	ecoll(:,17)= (/0.037, 0.037, 0.037, 0.047, 0.057, 0.250, 0.490, &
657	0.770, 0.890, 0.950, 1.000, 1.000, 1.000, 1.000, 1.000/)
658	ecoll(:,18)= (/0.036, 0.036, 0.036, 0.042, 0.048, 0.230, 0.470, &
659	0.780, 0.920, 1.000, 1.020, 1.020, 1.020, 1.020, 1.020/)
660	ecoll(:,19)= (/0.040, 0.040, 0.035, 0.033, 0.040, 0.112, 0.450, &
661	0.790, 1.010, 1.030, 1.040, 1.040, 1.040, 1.040, 1.040/)
662	ecoll(:,20)= (/0.033, 0.033, 0.033, 0.033, 0.033, 0.119, 0.470, &
663	0.950, 1.300, 1.700, 2.300, 2.300, 2.300, 2.300, 2.300/)
664	ecoll(:,21)= (/0.027, 0.027, 0.027, 0.027, 0.027, 0.125, 0.520, &
665	1.400, 2.300, 3.000, 4.000, 4.000, 4.000, 4.000, 4.000/)
666	ENDIF
667
668	!
669	!-- Calculate the radius class index of particles with respect to array r
670	ALLOCATE( ira(1:radius_classes) )
671	DO j = 1, radius_classes
672	particle_radius = radclass(j) * 1.0E4 ! in microm
673	DO k = 1, 15
674	IF ( particle_radius < r0(k) ) THEN
675	ira(j) = k
676	EXIT
677	ENDIF
678	ENDDO
679	IF ( particle_radius >= r0(15) ) ira(j) = 16
680	ENDDO
681
682	!
683	!-- Two-dimensional linear interpolation of the collision efficiency.
684	!-- Radius has to be in Âµm
685	DO j = 1, radius_classes
686	DO i = 1, j
687
688	ir = ira(j)
689	rq = radclass(i) / radclass(j)
690	iq = INT( rq * 20 ) + 1
691	iq = MAX( iq , 2)
692
693	IF ( ir < 16 ) THEN
694	IF ( ir >= 2 ) THEN
695	pp = ( ( radclass(j) * 1.0E04 ) - r0(ir-1) ) / &
696	( r0(ir) - r0(ir-1) )
697	qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
698	ec(j,i) = ( 1.0-pp ) * ( 1.0-qq ) * ecoll(ir-1,iq-1) &
699	+ pp * ( 1.0-qq ) * ecoll(ir,iq-1) &
700	+ qq * ( 1.0-pp ) * ecoll(ir-1,iq) &
701	+ pp * qq * ecoll(ir,iq)
702	ELSE
703	qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
704	ec(j,i) = (1.0-qq) * ecoll(1,iq-1) + qq * ecoll(1,iq)
705	ENDIF
706	ELSE
707	qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
708	ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq)
709	ec(j,i) = MIN( ek, 1.0 )
710	ENDIF
711
712	ec(i,j) = ec(j,i)
713	IF ( ec(i,j) < 1.0E-20 ) ec(i,j) = 0.0
714
715	ENDDO
716	ENDDO
717
718	DEALLOCATE( ira )
719
720	END SUBROUTINE effic
721
722
723	!------------------------------------------------------------------------------!
724	! Calculation of enhancement factor for collision efficencies due to turbulence
725	!------------------------------------------------------------------------------!
726	SUBROUTINE turb_enhance_eff
727
728	USE constants
729	USE cloud_parameters
730	USE particle_attributes
731	USE arrays_3d
732
733	IMPLICIT NONE
734
735	INTEGER :: i, ik, iq, ir, j, k, kk
736
737	INTEGER, DIMENSION(:), ALLOCATABLE :: ira
738
739	REAL :: particle_radius, pp, qq, rq, x1, x2, x3, y1, y2, y3
740
741	LOGICAL, SAVE :: first = .TRUE.
742
743	REAL, DIMENSION(1:11), SAVE :: rat
744	REAL, DIMENSION(1:7), SAVE :: r0
745	REAL, DIMENSION(1:7,1:11), SAVE :: ecoll_100, ecoll_400
746
747	!
748	!-- Initial assignment of constants
749	IF ( first ) THEN
750
751	first = .FALSE.
752
753	r0 = (/ 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 100.0 /)
754	rat = (/ 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 /)
755	!
756	!-- In 100 cm^2/s^3
757	ecoll_100(:,1) = (/1.74, 1.74, 1.773, 1.49, 1.207, 1.207, 1.0 /)
758	ecoll_100(:,2) = (/1.46, 1.46, 1.421, 1.245, 1.069, 1.069, 1.0 /)
759	ecoll_100(:,3) = (/1.32, 1.32, 1.245, 1.123, 1.000, 1.000, 1.0 /)
760	ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025, 1.025, 1.0 /)
761	ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056, 1.056, 1.0 /)
762	ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028, 1.028, 1.0 /)
763	ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046, 1.046, 1.0 /)
764	ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029, 1.029, 1.0 /)
765	ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021, 1.021, 1.0 /)
766	ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088, 1.088, 1.0 /)
767	ecoll_100(:,11)= (/20.3, 20.3, 14.6 , 8.61, 2.60, 2.60 , 1.0 /)
768	!
769	!-- In 400 cm^2/s^3
770	ecoll_400(:,1) = (/4.976, 4.976, 3.593, 2.519, 1.445, 1.445, 1.0 /)
771	ecoll_400(:,2) = (/2.984, 2.984, 2.181, 1.691, 1.201, 1.201, 1.0 /)
772	ecoll_400(:,3) = (/1.988, 1.988, 1.475, 1.313, 1.150, 1.150, 1.0 /)
773	ecoll_400(:,4) = (/1.490, 1.490, 1.187, 1.156, 1.126, 1.126, 1.0 /)
774	ecoll_400(:,5) = (/1.249, 1.249, 1.088, 1.090, 1.092, 1.092, 1.0 /)
775	ecoll_400(:,6) = (/1.139, 1.139, 1.130, 1.091, 1.051, 1.051, 1.0 /)
776	ecoll_400(:,7) = (/1.220, 1.220, 1.190, 1.138, 1.086, 1.086, 1.0 /)
777	ecoll_400(:,8) = (/1.325, 1.325, 1.267, 1.165, 1.063, 1.063, 1.0 /)
778	ecoll_400(:,9) = (/1.716, 1.716, 1.345, 1.223, 1.100, 1.100, 1.0 /)
779	ecoll_400(:,10)= (/3.788, 3.788, 1.501, 1.311, 1.120, 1.120, 1.0 /)
780	ecoll_400(:,11)= (/36.52, 36.52, 19.16, 22.80, 26.0, 26.0, 1.0 /)
781
782	ENDIF
783
784	!
785	!-- Calculate the radius class index of particles with respect to array r0
786	ALLOCATE( ira(1:radius_classes) )
787
788	DO j = 1, radius_classes
789	particle_radius = radclass(j) * 1.0E4 ! in microm
790	DO k = 1, 7
791	IF ( particle_radius < r0(k) ) THEN
792	ira(j) = k
793	EXIT
794	ENDIF
795	ENDDO
796	IF ( particle_radius >= r0(7) ) ira(j) = 8
797	ENDDO
798
799	!
800	!-- Two-dimensional linear interpolation of the collision efficiencies
801	DO j = 1, radius_classes
802	DO i = 1, j
803
804	ir = ira(j)
805	rq = radclass(i) / radclass(j)
806
807	DO kk = 2, 11
808	IF ( rq <= rat(kk) ) THEN
809	iq = kk
810	EXIT
811	ENDIF
812	ENDDO
813
814	y1 = 1.0 ! 0 cm2/s3
815	!
816	!-- 100 cm2/s3, 400 cm2/s3
817	IF ( ir < 8 ) THEN
818	IF ( ir >= 2 ) THEN
819	pp = ( radclass(j)*1.0E4 - r0(ir-1) ) / ( r0(ir) - r0(ir-1) )
820	qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
821	y2 = ( 1.0-pp ) * ( 1.0-qq ) * ecoll_100(ir-1,iq-1) + &
822	pp * ( 1.0-qq ) * ecoll_100(ir,iq-1) + &
823	qq * ( 10.-pp ) * ecoll_100(ir-1,iq) + &
824	pp * qq * ecoll_100(ir,iq)
825	y3 = ( 1.0-pp ) * ( 1.0-qq ) * ecoll_400(ir-1,iq-1) + &
826	pp * ( 1.0-qq ) * ecoll_400(ir,iq-1) + &
827	qq * ( 1.0-pp ) * ecoll_400(ir-1,iq) + &
828	pp * qq * ecoll_400(ir,iq)
829	ELSE
830	qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
831	y2 = ( 1.0-qq ) * ecoll_100(1,iq-1) + qq * ecoll_100(1,iq)
832	y3 = ( 1.0-qq ) * ecoll_400(1,iq-1) + qq * ecoll_400(1,iq)
833	ENDIF
834	ELSE
835	qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
836	y2 = ( 1.0-qq ) * ecoll_100(7,iq-1) + qq * ecoll_100(7,iq)
837	y3 = ( 1.0-qq ) * ecoll_400(7,iq-1) + qq * ecoll_400(7,iq)
838	ENDIF
839	!
840	!-- Linear interpolation of dissipation rate in cm2/s3
841	IF ( epsilon <= 100.0 ) THEN
842	ecf(j,i) = ( epsilon - 100.0 ) / ( 0.0 - 100.0 ) * y1 &
843	+ ( epsilon - 0.0 ) / ( 100.0 - 0.0 ) * y2
844	ELSEIF ( epsilon <= 600.0 ) THEN
845	ecf(j,i) = ( epsilon - 400.0 ) / ( 100.0 - 400.0 ) * y2 &
846	+ ( epsilon - 100.0 ) / ( 400.0 - 100.0 ) * y3
847	ELSE
848	ecf(j,i) = ( 600.0 - 400.0 ) / ( 100.0 - 400.0 ) * y2 &
849	+ ( 600.0 - 100.0 ) / ( 400.0 - 100.0 ) * y3
850	ENDIF
851
852	IF ( ecf(j,i) < 1.0 ) ecf(j,i) = 1.0
853
854	ecf(i,j) = ecf(j,i)
855
856	ENDDO
857	ENDDO
858
859	END SUBROUTINE turb_enhance_eff
860
861
862
863	SUBROUTINE collision_efficiency_rogers( mean_r, r, e)
864	!------------------------------------------------------------------------------!
865	! Collision efficiencies from table 8.2 in Rogers and Yau (1989, 3rd edition).
866	! Values are calculated from table by bilinear interpolation.
867	!------------------------------------------------------------------------------!
868
869	IMPLICIT NONE
870
871	INTEGER :: i, j, k
872
873	LOGICAL, SAVE :: first = .TRUE.
874
875	REAL :: aa, bb, cc, dd, dx, dy, e, gg, mean_r, mean_rm, r, &
876	rm, x, y
877
878	REAL, DIMENSION(1:9), SAVE :: collected_r = 0.0
879	REAL, DIMENSION(1:19), SAVE :: collector_r = 0.0
880	REAL, DIMENSION(1:9,1:19), SAVE :: ef = 0.0
881
882	mean_rm = mean_r * 1.0E06
883	rm = r * 1.0E06
884
885	IF ( first ) THEN
886
887	collected_r = (/ 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0, 20.0, 25.0 /)
888	collector_r = (/ 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 80.0, 100.0, &
889	150.0, 200.0, 300.0, 400.0, 500.0, 600.0, 1000.0, &
890	1400.0, 1800.0, 2400.0, 3000.0 /)
891
892	ef(:,1) = (/0.017, 0.027, 0.037, 0.052, 0.052, 0.052, 0.052, 0.0, &
893	0.0 /)
894	ef(:,2) = (/0.001, 0.016, 0.027, 0.060, 0.12, 0.17, 0.17, 0.17, 0.0 /)
895	ef(:,3) = (/0.001, 0.001, 0.02, 0.13, 0.28, 0.37, 0.54, 0.55, 0.47/)
896	ef(:,4) = (/0.001, 0.001, 0.02, 0.23, 0.4, 0.55, 0.7, 0.75, 0.75/)
897	ef(:,5) = (/0.01, 0.01, 0.03, 0.3, 0.4, 0.58, 0.73, 0.75, 0.79/)
898	ef(:,6) = (/0.01, 0.01, 0.13, 0.38, 0.57, 0.68, 0.80, 0.86, 0.91/)
899	ef(:,7) = (/0.01, 0.085, 0.23, 0.52, 0.68, 0.76, 0.86, 0.92, 0.95/)
900	ef(:,8) = (/0.01, 0.14, 0.32, 0.60, 0.73, 0.81, 0.90, 0.94, 0.96/)
901	ef(:,9) = (/0.025, 0.25, 0.43, 0.66, 0.78, 0.83, 0.92, 0.95, 0.96/)
902	ef(:,10)= (/0.039, 0.3, 0.46, 0.69, 0.81, 0.87, 0.93, 0.95, 0.96/)
903	ef(:,11)= (/0.095, 0.33, 0.51, 0.72, 0.82, 0.87, 0.93, 0.96, 0.97/)
904	ef(:,12)= (/0.098, 0.36, 0.51, 0.73, 0.83, 0.88, 0.93, 0.96, 0.97/)
905	ef(:,13)= (/0.1, 0.36, 0.52, 0.74, 0.83, 0.88, 0.93, 0.96, 0.97/)
906	ef(:,14)= (/0.17, 0.4, 0.54, 0.72, 0.83, 0.88, 0.94, 0.98, 1.0 /)
907	ef(:,15)= (/0.15, 0.37, 0.52, 0.74, 0.82, 0.88, 0.94, 0.98, 1.0 /)
908	ef(:,16)= (/0.11, 0.34, 0.49, 0.71, 0.83, 0.88, 0.94, 0.95, 1.0 /)
909	ef(:,17)= (/0.08, 0.29, 0.45, 0.68, 0.8, 0.86, 0.96, 0.94, 1.0 /)
910	ef(:,18)= (/0.04, 0.22, 0.39, 0.62, 0.75, 0.83, 0.92, 0.96, 1.0 /)
911	ef(:,19)= (/0.02, 0.16, 0.33, 0.55, 0.71, 0.81, 0.90, 0.94, 1.0 /)
912
913	ENDIF
914
915	DO k = 1, 8
916	IF ( collected_r(k) <= mean_rm ) i = k
917	ENDDO
918
919	DO k = 1, 18
920	IF ( collector_r(k) <= rm ) j = k
921	ENDDO
922
923	IF ( rm < 10.0 ) THEN
924	e = 0.0
925	ELSEIF ( mean_rm < 2.0 ) THEN
926	e = 0.001
927	ELSEIF ( mean_rm >= 25.0 ) THEN
928	IF( j <= 2 ) e = 0.0
929	IF( j == 3 ) e = 0.47
930	IF( j == 4 ) e = 0.8
931	IF( j == 5 ) e = 0.9
932	IF( j >=6 ) e = 1.0
933	ELSEIF ( rm >= 3000.0 ) THEN
934	IF( i == 1 ) e = 0.02
935	IF( i == 2 ) e = 0.16
936	IF( i == 3 ) e = 0.33
937	IF( i == 4 ) e = 0.55
938	IF( i == 5 ) e = 0.71
939	IF( i == 6 ) e = 0.81
940	IF( i == 7 ) e = 0.90
941	IF( i >= 8 ) e = 0.94
942	ELSE
943	x = mean_rm - collected_r(i)
944	y = rm - collector_r(j)
945	dx = collected_r(i+1) - collected_r(i)
946	dy = collector_r(j+1) - collector_r(j)
947	aa = x2 + y2
948	bb = ( dx - x )2 + y2
949	cc = x2 + ( dy - y )2
950	dd = ( dx - x )2 + ( dy - y )2
951	gg = aa + bb + cc + dd
952
953	e = ( (gg-aa)ef(i,j) + (gg-bb)ef(i+1,j) + (gg-cc)*ef(i,j+1) + &
954	(gg-dd)ef(i+1,j+1) ) / (3.0gg)
955	ENDIF
956
957	END SUBROUTINE collision_efficiency_rogers
958
959	END MODULE lpm_collision_kernels_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |