Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |