Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |