Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |