Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |