Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |