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

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