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

Last change on this file since 1861 was 1861, checked in by hoffmann, 8 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 hoffmann $
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.