Home

Context Navigation

source: palm/trunk/SOURCE/lpm_droplet_collision.f90 @ 1928

Last change on this file since 1928 was 1885, checked in by hoffmann, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 22.7 KB

Line
1	!> @file lpm_droplet_collision.f90
2	!--------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the terms
6	! of the GNU General Public License as published by the Free Software Foundation,
7	! either version 3 of the License, or (at your option) any later version.
8	!
9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12	!
13	! You should have received a copy of the GNU General Public License along with
14	! PALM. If not, see <http://www.gnu.org/licenses/>.
15	!
16	! Copyright 1997-2016 Leibniz Universitaet Hannover
17	!--------------------------------------------------------------------------------!
18	!
19	! Current revisions:
20	! ------------------
21	!
22	!
23	! Former revisions:
24	! -----------------
25	! $Id: lpm_droplet_collision.f90 1885 2016-04-21 11:12:26Z hellstea $
26	!
27	! 1884 2016-04-21 11:11:40Z hoffmann
28	! Conservation of mass should only be checked if collisions took place.
29	!
30	! 1860 2016-04-13 13:21:28Z hoffmann
31	! Interpolation of dissipation rate adjusted to more reasonable values.
32	!
33	! 1822 2016-04-07 07:49:42Z hoffmann
34	! Integration of a new collision algortithm based on Shima et al. (2009) and
35	! Soelch and Kaercher (2010) called all_or_nothing. The previous implemented
36	! collision algorithm is called average_impact. Moreover, both algorithms are
37	! now positive definit due to their construction, i.e., no negative weighting
38	! factors should occur.
39	!
40	! 1682 2015-10-07 23:56:08Z knoop
41	! Code annotations made doxygen readable
42	!
43	! 1359 2014-04-11 17:15:14Z hoffmann
44	! New particle structure integrated.
45	! Kind definition added to all floating point numbers.
46	!
47	! 1322 2014-03-20 16:38:49Z raasch
48	! REAL constants defined as wp_kind
49	!
50	! 1320 2014-03-20 08:40:49Z raasch
51	! ONLY-attribute added to USE-statements,
52	! kind-parameters added to all INTEGER and REAL declaration statements,
53	! kinds are defined in new module kinds,
54	! revision history before 2012 removed,
55	! comment fields (!:) to be used for variable explanations added to
56	! all variable declaration statements
57	!
58	! 1092 2013-02-02 11:24:22Z raasch
59	! unused variables removed
60	!
61	! 1071 2012-11-29 16:54:55Z franke
62	! Calculation of Hall and Wang kernel now uses collision-coalescence formulation
63	! proposed by Wang instead of the continuous collection equation (for more
64	! information about new method see PALM documentation)
65	! Bugfix: message identifiers added
66	!
67	! 1036 2012-10-22 13:43:42Z raasch
68	! code put under GPL (PALM 3.9)
69	!
70	! 849 2012-03-15 10:35:09Z raasch
71	! initial revision (former part of advec_particles)
72	!
73	!
74	! Description:
75	! ------------
76	!> Calculates change in droplet radius by collision. Droplet collision is
77	!> calculated for each grid box seperately. Collision is parameterized by
78	!> using collision kernels. Two different kernels are available:
79	!> Hall kernel: Kernel from Hall (1980, J. Atmos. Sci., 2486-2507), which
80	!> considers collision due to pure gravitational effects.
81	!> Wang kernel: Beside gravitational effects (treated with the Hall-kernel) also
82	!> the effects of turbulence on the collision are considered using
83	!> parameterizations of Ayala et al. (2008, New J. Phys., 10,
84	!> 075015) and Wang and Grabowski (2009, Atmos. Sci. Lett., 10,
85	!> 1-8). This kernel includes three possible effects of turbulence:
86	!> the modification of the relative velocity between the droplets,
87	!> the effect of preferential concentration, and the enhancement of
88	!> collision efficiencies.
89	!------------------------------------------------------------------------------!
90	SUBROUTINE lpm_droplet_collision (i,j,k)
91
92
93
94	USE arrays_3d, &
95	ONLY: diss, ql_v, ql_vp
96
97	USE cloud_parameters, &
98	ONLY: rho_l
99
100	USE constants, &
101	ONLY: pi
102
103	USE control_parameters, &
104	ONLY: dt_3d, message_string, dz
105
106	USE cpulog, &
107	ONLY: cpu_log, log_point_s
108
109	USE grid_variables, &
110	ONLY: dx, dy
111
112	USE kinds
113
114	USE lpm_collision_kernels_mod, &
115	ONLY: ckernel, recalculate_kernel
116
117	USE particle_attributes, &
118	ONLY: all_or_nothing, average_impact, dissipation_classes, &
119	hall_kernel, iran_part, number_of_particles, particles, &
120	particle_type, prt_count, use_kernel_tables, wang_kernel
121
122	USE random_function_mod, &
123	ONLY: random_function
124
125	USE pegrid
126
127	IMPLICIT NONE
128
129	INTEGER(iwp) :: eclass !<
130	INTEGER(iwp) :: i !<
131	INTEGER(iwp) :: j !<
132	INTEGER(iwp) :: k !<
133	INTEGER(iwp) :: n !<
134	INTEGER(iwp) :: m !<
135	INTEGER(iwp) :: rclass_l !<
136	INTEGER(iwp) :: rclass_s !<
137
138	REAL(wp) :: collection_probability !< probability for collection
139	REAL(wp) :: ddV !< inverse grid box volume
140	REAL(wp) :: epsilon !< dissipation rate
141	REAL(wp) :: factor_volume_to_mass !< 4.0 / 3.0 * pi * rho_l
142	REAL(wp) :: xm !< mean mass of droplet m
143	REAL(wp) :: xn !< mean mass of droplet n
144
145	REAL(wp), DIMENSION(:), ALLOCATABLE :: weight !< weighting factor
146	REAL(wp), DIMENSION(:), ALLOCATABLE :: mass !< total mass of super droplet
147
148	CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'start' )
149
150	number_of_particles = prt_count(k,j,i)
151	factor_volume_to_mass = 4.0_wp / 3.0_wp * pi * rho_l
152	ddV = 1 / ( dx * dy * dz )
153	!
154	!-- Collision requires at least one super droplet inside the box
155	IF ( number_of_particles > 0 ) THEN
156
157	!
158	!-- Now apply the different kernels
159	IF ( use_kernel_tables ) THEN
160	!
161	!-- Fast method with pre-calculated collection kernels for
162	!-- discrete radius- and dissipation-classes.
163	!--
164	!-- Determine dissipation class index of this gridbox
165	IF ( wang_kernel ) THEN
166	eclass = INT( diss(k,j,i) * 1.0E4_wp / 600.0_wp * &
167	dissipation_classes ) + 1
168	epsilon = diss(k,j,i)
169	ELSE
170	epsilon = 0.0_wp
171	ENDIF
172	IF ( hall_kernel .OR. epsilon * 1.0E4_wp < 0.001_wp ) THEN
173	eclass = 0 ! Hall kernel is used
174	ELSE
175	eclass = MIN( dissipation_classes, eclass )
176	ENDIF
177
178	!
179	!-- Droplet collision are calculated using collision-coalescence
180	!-- formulation proposed by Wang (see PALM documentation)
181	!-- Temporary fields for total mass of super-droplet and weighting factors
182	!-- are allocated.
183	ALLOCATE(mass(1:number_of_particles), weight(1:number_of_particles))
184
185	mass(1:number_of_particles) = particles(1:number_of_particles)%weight_factor * &
186	particles(1:number_of_particles)%radius*3 &
187	factor_volume_to_mass
188	weight(1:number_of_particles) = particles(1:number_of_particles)%weight_factor
189
190	IF ( average_impact ) THEN ! select collision algorithm
191
192	DO n = 1, number_of_particles
193
194	rclass_l = particles(n)%class
195	xn = mass(n) / weight(n)
196
197	DO m = n, number_of_particles
198
199	rclass_s = particles(m)%class
200	xm = mass(m) / weight(m)
201
202	IF ( xm .LT. xn ) THEN
203
204	!
205	!-- Particle n collects smaller particle m
206	collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
207	weight(n) * ddV * dt_3d
208
209	mass(n) = mass(n) + mass(m) * collection_probability
210	weight(m) = weight(m) - weight(m) * collection_probability
211	mass(m) = mass(m) - mass(m) * collection_probability
212	ELSEIF ( xm .GT. xn ) THEN
213	!
214	!-- Particle m collects smaller particle n
215	collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
216	weight(m) * ddV * dt_3d
217
218	mass(m) = mass(m) + mass(n) * collection_probability
219	weight(n) = weight(n) - weight(n) * collection_probability
220	mass(n) = mass(n) - mass(n) * collection_probability
221	ELSE
222	!
223	!-- Same-size collections. If n = m, weight is reduced by the
224	!-- number of possible same-size collections; the total mass
225	!-- is not changed during same-size collection.
226	!-- Same-size collections of different
227	!-- particles ( n /= m ) are treated as same-size collections
228	!-- of ONE partilce with weight = weight(n) + weight(m) and
229	!-- mass = mass(n) + mass(m).
230	!-- Accordingly, each particle loses the same number of
231	!-- droplets to the other particle, but this has no effect on
232	!-- total mass mass, since the exchanged droplets have the
233	!-- same radius.
234
235	!-- Note: For m = n this equation is an approximation only
236	!-- valid for weight >> 1 (which is usually the case). The
237	!-- approximation is weight(n)-1 = weight(n).
238	weight(n) = weight(n) - 0.5_wp * weight(n) * &
239	ckernel(rclass_l,rclass_s,eclass) * &
240	weight(m) * ddV * dt_3d
241	IF ( n .NE. m ) THEN
242	weight(m) = weight(m) - 0.5_wp * weight(m) * &
243	ckernel(rclass_l,rclass_s,eclass) * &
244	weight(n) * ddV * dt_3d
245	ENDIF
246	ENDIF
247
248	ENDDO
249
250	ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
251
252	ENDDO
253
254	ELSEIF ( all_or_nothing ) THEN ! select collision algorithm
255
256	DO n = 1, number_of_particles
257
258	rclass_l = particles(n)%class
259	xn = mass(n) / weight(n) ! mean mass of droplet n
260
261	DO m = n, number_of_particles
262
263	rclass_s = particles(m)%class
264	xm = mass(m) / weight(m) ! mean mass of droplet m
265
266	IF ( weight(n) .LT. weight(m) ) THEN
267	!
268	!-- Particle n collects weight(n) droplets of particle m
269	collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
270	weight(m) * ddV * dt_3d
271
272	IF ( collection_probability .GT. random_function( iran_part ) ) THEN
273	mass(n) = mass(n) + weight(n) * xm
274	weight(m) = weight(m) - weight(n)
275	mass(m) = mass(m) - weight(n) * xm
276	ENDIF
277
278	ELSEIF ( weight(m) .LT. weight(n) ) THEN
279	!
280	!-- Particle m collects weight(m) droplets of particle n
281	collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
282	weight(n) * ddV * dt_3d
283
284	IF ( collection_probability .GT. random_function( iran_part ) ) THEN
285	mass(m) = mass(m) + weight(m) * xn
286	weight(n) = weight(n) - weight(m)
287	mass(n) = mass(n) - weight(m) * xn
288	ENDIF
289	ELSE
290	!
291	!-- Collisions of particles of the same weighting factor.
292	!-- Particle n collects 1/2 weight(n) droplets of particle m,
293	!-- particle m collects 1/2 weight(m) droplets of particle n.
294	!-- The total mass mass changes accordingly.
295	!-- If n = m, the first half of the droplets coalesces with the
296	!-- second half of the droplets; mass is unchanged because
297	!-- xm = xn for n = m.
298
299	!-- Note: For m = n this equation is an approximation only
300	!-- valid for weight >> 1 (which is usually the case). The
301	!-- approximation is weight(n)-1 = weight(n).
302	collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
303	weight(n) * ddV * dt_3d
304
305	IF ( collection_probability .GT. random_function( iran_part ) ) THEN
306	mass(n) = mass(n) + 0.5_wp * weight(n) * ( xm - xn )
307	mass(m) = mass(m) + 0.5_wp * weight(m) * ( xn - xm )
308	weight(n) = weight(n) - 0.5_wp * weight(m)
309	weight(m) = weight(n)
310	ENDIF
311	ENDIF
312
313	ENDDO
314
315	ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
316
317	ENDDO
318
319	ENDIF
320
321
322
323
324	IF ( ANY(weight < 0.0_wp) ) THEN
325	WRITE( message_string, * ) 'negative weighting'
326	CALL message( 'lpm_droplet_collision', 'PA0028', &
327	2, 2, -1, 6, 1 )
328	ENDIF
329
330	particles(1:number_of_particles)%radius = ( mass(1:number_of_particles) / &
331	( weight(1:number_of_particles) &
332	* factor_volume_to_mass &
333	) &
334	)**0.33333333333333_wp
335
336	particles(1:number_of_particles)%weight_factor = weight(1:number_of_particles)
337
338	DEALLOCATE(weight, mass)
339
340	ELSEIF ( .NOT. use_kernel_tables ) THEN
341	!
342	!-- Collection kernels are calculated for every new
343	!-- grid box. First, allocate memory for kernel table.
344	!-- Third dimension is 1, because table is re-calculated for
345	!-- every new dissipation value.
346	ALLOCATE( ckernel(1:number_of_particles,1:number_of_particles,1:1) )
347	!
348	!-- Now calculate collection kernel for this box. Note that
349	!-- the kernel is based on the previous time step
350	CALL recalculate_kernel( i, j, k )
351	!
352	!-- Droplet collision are calculated using collision-coalescence
353	!-- formulation proposed by Wang (see PALM documentation)
354	!-- Temporary fields for total mass of super-droplet and weighting factors
355	!-- are allocated.
356	ALLOCATE(mass(1:number_of_particles), weight(1:number_of_particles))
357
358	mass(1:number_of_particles) = particles(1:number_of_particles)%weight_factor * &
359	particles(1:number_of_particles)%radius*3 &
360	factor_volume_to_mass
361
362	weight(1:number_of_particles) = particles(1:number_of_particles)%weight_factor
363
364	IF ( average_impact ) THEN ! select collision algorithm
365
366	DO n = 1, number_of_particles
367
368	xn = mass(n) / weight(n) ! mean mass of droplet n
369
370	DO m = n, number_of_particles
371
372	xm = mass(m) / weight(m) !mean mass of droplet m
373
374	IF ( xm .LT. xn ) THEN
375	!
376	!-- Particle n collects smaller particle m
377	collection_probability = ckernel(n,m,1) * weight(n) * &
378	ddV * dt_3d
379
380	mass(n) = mass(n) + mass(m) * collection_probability
381	weight(m) = weight(m) - weight(m) * collection_probability
382	mass(m) = mass(m) - mass(m) * collection_probability
383	ELSEIF ( xm .GT. xn ) THEN
384	!
385	!-- Particle m collects smaller particle n
386	collection_probability = ckernel(n,m,1) * weight(m) * &
387	ddV * dt_3d
388
389	mass(m) = mass(m) + mass(n) * collection_probability
390	weight(n) = weight(n) - weight(n) * collection_probability
391	mass(n) = mass(n) - mass(n) * collection_probability
392	ELSE
393	!
394	!-- Same-size collections. If n = m, weight is reduced by the
395	!-- number of possible same-size collections; the total mass
396	!-- mass is not changed during same-size collection.
397	!-- Same-size collections of different
398	!-- particles ( n /= m ) are treated as same-size collections
399	!-- of ONE partilce with weight = weight(n) + weight(m) and
400	!-- mass = mass(n) + mass(m).
401	!-- Accordingly, each particle loses the same number of
402	!-- droplets to the other particle, but this has no effect on
403	!-- total mass mass, since the exchanged droplets have the
404	!-- same radius.
405	!--
406	!-- Note: For m = n this equation is an approximation only
407	!-- valid for weight >> 1 (which is usually the case). The
408	!-- approximation is weight(n)-1 = weight(n).
409	weight(n) = weight(n) - 0.5_wp * weight(n) * &
410	ckernel(n,m,1) * weight(m) * &
411	ddV * dt_3d
412	IF ( n .NE. m ) THEN
413	weight(m) = weight(m) - 0.5_wp * weight(m) * &
414	ckernel(n,m,1) * weight(n) * &
415	ddV * dt_3d
416	ENDIF
417	ENDIF
418
419
420	ENDDO
421
422	ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
423
424	ENDDO
425
426	ELSEIF ( all_or_nothing ) THEN ! select collision algorithm
427
428	DO n = 1, number_of_particles
429
430	xn = mass(n) / weight(n) ! mean mass of droplet n
431
432	DO m = n, number_of_particles
433
434	xm = mass(m) / weight(m) !mean mass of droplet m
435
436	IF ( weight(n) .LT. weight(m) ) THEN
437	!
438	!-- Particle n collects smaller particle m
439	collection_probability = ckernel(n,m,1) * weight(m) * &
440	ddV * dt_3d
441
442	IF ( collection_probability .GT. random_function( iran_part ) ) THEN
443	mass(n) = mass(n) + weight(n) * xm
444	weight(m) = weight(m) - weight(n)
445	mass(m) = mass(m) - weight(n) * xm
446	ENDIF
447
448	ELSEIF ( weight(m) .LT. weight(n) ) THEN
449	!
450	!-- Particle m collects smaller particle n
451	collection_probability = ckernel(n,m,1) * weight(n) * &
452	ddV * dt_3d
453
454	IF ( collection_probability .GT. random_function( iran_part ) ) THEN
455	mass(m) = mass(m) + weight(m) * xn
456	weight(n) = weight(n) - weight(m)
457	mass(n) = mass(n) - weight(m) * xn
458	ENDIF
459	ELSE
460	!
461	!-- Collisions of particles of the same weighting factor.
462	!-- Particle n collects 1/2 weight(n) droplets of particle m,
463	!-- particle m collects 1/2 weight(m) droplets of particle n.
464	!-- The total mass mass changes accordingly.
465	!-- If n = m, the first half of the droplets coalesces with the
466	!-- second half of the droplets; mass is unchanged because
467	!-- xm = xn for n = m.
468	!--
469	!-- Note: For m = n this equation is an approximation only
470	!-- valid for weight >> 1 (which is usually the case). The
471	!-- approximation is weight(n)-1 = weight(n).
472	collection_probability = ckernel(n,m,1) * weight(n) * &
473	ddV * dt_3d
474
475	IF ( collection_probability .GT. random_function( iran_part ) ) THEN
476	mass(n) = mass(n) + 0.5_wp * weight(n) * ( xm - xn )
477	mass(m) = mass(m) + 0.5_wp * weight(m) * ( xn - xm )
478	weight(n) = weight(n) - 0.5_wp * weight(m)
479	weight(m) = weight(n)
480	ENDIF
481	ENDIF
482
483
484	ENDDO
485
486	ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
487
488	ENDDO
489
490	ENDIF
491
492	IF ( ANY(weight < 0.0_wp) ) THEN
493	WRITE( message_string, * ) 'negative weighting'
494	CALL message( 'lpm_droplet_collision', 'PA0028', &
495	2, 2, -1, 6, 1 )
496	ENDIF
497
498	particles(1:number_of_particles)%radius = ( mass(1:number_of_particles) / &
499	( weight(1:number_of_particles) &
500	* factor_volume_to_mass &
501	) &
502	)**0.33333333333333_wp
503
504	particles(1:number_of_particles)%weight_factor = weight(1:number_of_particles)
505
506	DEALLOCATE( weight, mass, ckernel )
507
508	ENDIF
509
510	!
511	!-- Check if LWC is conserved during collision process
512	IF ( ql_v(k,j,i) /= 0.0_wp ) THEN
513	IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001_wp .OR. &
514	ql_vp(k,j,i) / ql_v(k,j,i) <= 0.9999_wp ) THEN
515	WRITE( message_string, * ) ' LWC is not conserved during', &
516	' collision! ', &
517	' LWC after condensation: ', ql_v(k,j,i), &
518	' LWC after collision: ', ql_vp(k,j,i)
519	CALL message( 'lpm_droplet_collision', 'PA0040', 2, 2, -1, 6, 1 )
520	ENDIF
521	ENDIF
522
523	ENDIF
524
525	CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'stop' )
526
527	END SUBROUTINE lpm_droplet_collision

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |