Home

Context Navigation

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

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