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

Last change on this file since 1682 was 1682, checked in by knoop, 6 years ago

Code annotations made doxygen readable

  • Property svn:keywords set to Id
File size: 27.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-2014 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! ------------------
21! Code annotations made doxygen readable
22!
23! Former revisions:
24! -----------------
25! $Id: lpm_droplet_collision.f90 1682 2015-10-07 23:56:08Z knoop $
26!
27! 1359 2014-04-11 17:15:14Z hoffmann
28! New particle structure integrated.
29! Kind definition added to all floating point numbers.
30!
31! 1322 2014-03-20 16:38:49Z raasch
32! REAL constants defined as wp_kind
33!
34! 1320 2014-03-20 08:40:49Z raasch
35! ONLY-attribute added to USE-statements,
36! kind-parameters added to all INTEGER and REAL declaration statements,
37! kinds are defined in new module kinds,
38! revision history before 2012 removed,
39! comment fields (!:) to be used for variable explanations added to
40! all variable declaration statements
41!
42! 1092 2013-02-02 11:24:22Z raasch
43! unused variables removed
44!
45! 1071 2012-11-29 16:54:55Z franke
46! Calculation of Hall and Wang kernel now uses collision-coalescence formulation
47! proposed by Wang instead of the continuous collection equation (for more
48! information about new method see PALM documentation)
49! Bugfix: message identifiers added
50!
51! 1036 2012-10-22 13:43:42Z raasch
52! code put under GPL (PALM 3.9)
53!
54! 849 2012-03-15 10:35:09Z raasch
55! initial revision (former part of advec_particles)
56!
57!
58! Description:
59! ------------
60!> Calculates change in droplet radius by collision. Droplet collision is
61!> calculated for each grid box seperately. Collision is parameterized by
62!> using collision kernels. Three different kernels are available:
63!> PALM kernel: Kernel is approximated using a method from Rogers and
64!>              Yau (1989, A Short Course in Cloud Physics, Pergamon Press).
65!>              All droplets smaller than the treated one are represented by
66!>              one droplet with mean features. Collision efficiencies are taken
67!>              from the respective table in Rogers and Yau.
68!> Hall kernel: Kernel from Hall (1980, J. Atmos. Sci., 2486-2507), which
69!>              considers collision due to pure gravitational effects.
70!> Wang kernel: Beside gravitational effects (treated with the Hall-kernel) also
71!>              the effects of turbulence on the collision are considered using
72!>              parameterizations of Ayala et al. (2008, New J. Phys., 10,
73!>              075015) and Wang and Grabowski (2009, Atmos. Sci. Lett., 10,
74!>              1-8). This kernel includes three possible effects of turbulence:
75!>              the modification of the relative velocity between the droplets,
76!>              the effect of preferential concentration, and the enhancement of
77!>              collision efficiencies.
78!------------------------------------------------------------------------------!
79 SUBROUTINE lpm_droplet_collision (i,j,k)
80 
81
82
83    USE arrays_3d,                                                             &
84        ONLY:  diss, ql, ql_v, ql_vp, u, v, w, zu, zw
85
86    USE cloud_parameters,                                                      &
87        ONLY:  effective_coll_efficiency
88
89    USE constants,                                                             &
90        ONLY:  pi
91
92    USE control_parameters,                                                    &
93        ONLY:  dt_3d, message_string, u_gtrans, v_gtrans, dz
94
95    USE cpulog,                                                                &
96        ONLY:  cpu_log, log_point_s
97
98    USE grid_variables,                                                        &
99        ONLY:  ddx, dx, ddy, dy
100
101    USE indices,                                                               &
102        ONLY:  nxl, nxr, nyn, nys, nzb, nzt
103
104    USE kinds
105
106    USE lpm_collision_kernels_mod,                                             &
107        ONLY:  ckernel, collision_efficiency_rogers, recalculate_kernel
108
109    USE particle_attributes,                                                   &
110        ONLY:  deleted_particles, dissipation_classes, hall_kernel,            &
111               palm_kernel, particles, particle_type,                          &
112               prt_count, use_kernel_tables, wang_kernel
113
114    USE pegrid
115
116    IMPLICIT NONE
117
118    INTEGER(iwp) ::  eclass   !<
119    INTEGER(iwp) ::  i        !<
120    INTEGER(iwp) ::  ii       !<
121    INTEGER(iwp) ::  inc      !<
122    INTEGER(iwp) ::  is       !<
123    INTEGER(iwp) ::  j        !<
124    INTEGER(iwp) ::  jj       !<
125    INTEGER(iwp) ::  js       !<
126    INTEGER(iwp) ::  k        !<
127    INTEGER(iwp) ::  kk       !<
128    INTEGER(iwp) ::  n        !<
129    INTEGER(iwp) ::  pse      !<
130    INTEGER(iwp) ::  psi      !<
131    INTEGER(iwp) ::  rclass_l !<
132    INTEGER(iwp) ::  rclass_s !<
133
134    INTEGER(iwp), DIMENSION(prt_count(k,j,i)) ::  rclass_v !<
135
136    LOGICAL, SAVE ::  first_flag = .TRUE. !<
137
138    TYPE(particle_type) :: tmp_particle   !<
139
140    REAL(wp) ::  aa       !<
141    REAL(wp) ::  auxn     !< temporary variables
142    REAL(wp) ::  auxs     !< temporary variables
143    REAL(wp) ::  bb       !<
144    REAL(wp) ::  cc       !<
145    REAL(wp) ::  dd       !<
146    REAL(wp) ::  ddV      !<
147    REAL(wp) ::  delta_r  !<
148    REAL(wp) ::  delta_v  !<
149    REAL(wp) ::  epsilon  !<
150    REAL(wp) ::  gg       !<
151    REAL(wp) ::  mean_r   !<
152    REAL(wp) ::  ql_int   !<
153    REAL(wp) ::  ql_int_l !<
154    REAL(wp) ::  ql_int_u !<
155    REAL(wp) ::  r3       !<
156    REAL(wp) ::  sl_r3    !<
157    REAL(wp) ::  sl_r4    !<
158    REAL(wp) ::  sum1     !<
159    REAL(wp) ::  sum2     !<
160    REAL(wp) ::  sum3     !<
161    REAL(wp) ::  u_int    !<
162    REAL(wp) ::  u_int_l  !<
163    REAL(wp) ::  u_int_u  !<
164    REAL(wp) ::  v_int    !<
165    REAL(wp) ::  v_int_l  !<
166    REAL(wp) ::  v_int_u  !<
167    REAL(wp) ::  w_int    !<
168    REAL(wp) ::  w_int_l  !<
169    REAL(wp) ::  w_int_u  !<
170    REAL(wp) ::  x        !<
171    REAL(wp) ::  y        !<
172
173    REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad    !<
174    REAL(wp), DIMENSION(:), ALLOCATABLE ::  weight !<
175
176    REAL, DIMENSION(prt_count(k,j,i))    :: ck
177    REAL, DIMENSION(prt_count(k,j,i))    :: r3v
178    REAL, DIMENSION(prt_count(k,j,i))    :: sum1v
179    REAL, DIMENSION(prt_count(k,j,i))    :: sum2v
180
181    CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'start' )
182
183!
184!-- Collision requires at least two particles in the box
185    IF ( prt_count(k,j,i) > 1 )  THEN
186!
187!--    First, sort particles within the gridbox by their size,
188!--    using Shell's method (see Numerical Recipes)
189!--    NOTE: In case of using particle tails, the re-sorting of
190!--    ----  tails would have to be included here!
191       IF ( .NOT. ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
192       use_kernel_tables ) )  THEN
193          psi = 0
194          inc = 1
195          DO WHILE ( inc <= prt_count(k,j,i) )
196             inc = 3 * inc + 1
197          ENDDO
198
199          DO WHILE ( inc > 1 )
200             inc = inc / 3
201             DO  is = inc+1, prt_count(k,j,i)
202                tmp_particle = particles(psi+is)
203                js = is
204                DO WHILE ( particles(psi+js-inc)%radius >             &
205                tmp_particle%radius )
206                   particles(psi+js) = particles(psi+js-inc)
207                   js = js - inc
208                   IF ( js <= inc )  EXIT
209                ENDDO
210                particles(psi+js) = tmp_particle
211             ENDDO
212          ENDDO
213       ENDIF
214
215       psi = 1
216       pse = prt_count(k,j,i)
217
218!
219!--    Now apply the different kernels
220       IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
221            use_kernel_tables )  THEN
222!
223!--       Fast method with pre-calculated efficiencies for
224!--       discrete radius- and dissipation-classes.
225!--
226!--       Determine dissipation class index of this gridbox
227          IF ( wang_kernel )  THEN
228             eclass = INT( diss(k,j,i) * 1.0E4_wp / 1000.0_wp * &
229                           dissipation_classes ) + 1
230             epsilon = diss(k,j,i)
231          ELSE
232             epsilon = 0.0_wp
233          ENDIF
234          IF ( hall_kernel  .OR.  epsilon * 1.0E4_wp < 0.001_wp )  THEN
235             eclass = 0   ! Hall kernel is used
236          ELSE
237             eclass = MIN( dissipation_classes, eclass )
238          ENDIF
239
240!
241!--       Droplet collision are calculated using collision-coalescence
242!--       formulation proposed by Wang (see PALM documentation)
243!--       Since new radii after collision are defined by radii of all
244!--       droplets before collision, temporary fields for new radii and
245!--       weighting factors are needed
246          ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i)))
247
248          rad    = 0.0_wp
249          weight = 0.0_wp
250
251          sum1v(1:prt_count(k,j,i)) = 0.0_wp
252          sum2v(1:prt_count(k,j,i)) = 0.0_wp
253
254          DO  n = 1, prt_count(k,j,i)
255
256             rclass_l = particles(n)%class
257!
258!--          Mass added due to collisions with smaller droplets
259             DO  is = n+1, prt_count(k,j,i)
260                rclass_s = particles(is)%class
261                auxs = ckernel(rclass_l,rclass_s,eclass) * particles(is)%weight_factor
262                auxn = ckernel(rclass_s,rclass_l,eclass) * particles(n)%weight_factor
263                IF ( particles(is)%radius <  particles(n)%radius )  THEN
264                   sum1v(n) =  sum1v(n)  + particles(is)%radius**3 * auxs
265                   sum2v(is) = sum2v(is) + auxn
266                ELSE
267                   sum2v(n)  = sum2v(n)  + auxs
268                   sum1v(is) = sum1v(is) + particles(n)%radius**3 * auxn
269                ENDIF
270             ENDDO
271          ENDDO
272          rclass_v = particles(1:prt_count(k,j,i))%class
273          DO  n = 1, prt_count(k,j,i)
274            ck(n)       = ckernel(rclass_v(n),rclass_v(n),eclass)
275          ENDDO
276          r3v      = particles(1:prt_count(k,j,i))%radius**3
277          DO  n = 1, prt_count(k,j,i)
278             sum3 = 0.0_wp
279             ddV = ddx * ddy / dz
280!
281!--          Change of the current weighting factor
282             sum3 = 1 - dt_3d * ddV *                                 &
283                        ck(n) *                                       &
284                        ( particles(n)%weight_factor - 1 ) * 0.5_wp - &
285                    dt_3d * ddV * sum2v(n)
286             weight(n) = particles(n)%weight_factor * sum3
287!
288!--          Change of the current droplet radius
289             rad(n) = ( (r3v(n) + dt_3d * ddV * (sum1v(n) - sum2v(n) * r3v(n)) )/&
290                              sum3 )**0.33333333333333_wp
291
292             ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n) &
293                                         * rad(n)**3
294
295          ENDDO
296          IF ( ANY(weight < 0.0_wp) )  THEN
297                WRITE( message_string, * ) 'negative weighting'
298                CALL message( 'lpm_droplet_collision', 'PA0028',      &
299                               2, 2, -1, 6, 1 )
300          ENDIF
301
302          particles(psi:pse)%radius = rad(1:prt_count(k,j,i))
303          particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i))
304
305          DEALLOCATE(rad, weight)
306
307       ELSEIF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
308                .NOT. use_kernel_tables )  THEN
309!
310!--       Collision efficiencies are calculated for every new
311!--       grid box. First, allocate memory for kernel table.
312!--       Third dimension is 1, because table is re-calculated for
313!--       every new dissipation value.
314          ALLOCATE( ckernel(1:prt_count(k,j,i),1:prt_count(k,j,i),1:1) )
315!
316!--       Now calculate collision efficiencies for this box
317          CALL recalculate_kernel( i, j, k )
318
319!
320!--       Droplet collision are calculated using collision-coalescence
321!--       formulation proposed by Wang (see PALM documentation)
322!--       Since new radii after collision are defined by radii of all
323!--       droplets before collision, temporary fields for new radii and
324!--       weighting factors are needed
325          ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i)))
326
327          rad = 0.0_wp
328          weight = 0.0_wp
329
330          DO  n = psi, pse, 1
331
332             sum1 = 0.0_wp
333             sum2 = 0.0_wp
334             sum3 = 0.0_wp
335!
336!--          Mass added due to collisions with smaller droplets
337             DO  is = psi, n-1
338                sum1 = sum1 + ( particles(is)%radius**3 *            &
339                                ckernel(n,is,1) *  &
340                                particles(is)%weight_factor )
341             ENDDO
342!
343!--          Rate of collisions with larger droplets
344             DO  is = n+1, pse
345                sum2 = sum2 + ( ckernel(n,is,1) *  &
346                                particles(is)%weight_factor )
347             ENDDO
348
349             r3 = particles(n)%radius**3
350             ddV = ddx * ddy / dz
351             is = 1
352!
353!--                   Change of the current weighting factor
354             sum3 = 1 - dt_3d * ddV *                                 &
355                        ckernel(n,n,1) *           &
356                        ( particles(n)%weight_factor - 1 ) * 0.5_wp - &
357                    dt_3d * ddV * sum2
358             weight(n-is+1) = particles(n)%weight_factor * sum3
359!
360!--                   Change of the current droplet radius
361             rad(n-is+1) = ( (r3 + dt_3d * ddV * (sum1 - sum2 * r3) )/&
362                              sum3 )**0.33333333333333_wp
363
364             IF ( weight(n-is+1) < 0.0_wp )  THEN
365                WRITE( message_string, * ) 'negative weighting',      &
366                                           'factor: ', weight(n-is+1)
367                CALL message( 'lpm_droplet_collision', 'PA0037',      &
368                               2, 2, -1, 6, 1 )
369             ENDIF
370
371             ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n-is+1) &
372                                         * rad(n-is+1)**3
373
374          ENDDO
375
376          particles(psi:pse)%radius = rad(1:prt_count(k,j,i))
377          particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i))
378
379          DEALLOCATE( rad, weight, ckernel )
380
381       ELSEIF ( palm_kernel )  THEN
382!
383!--       PALM collision kernel
384!
385!--       Calculate the mean radius of all those particles which
386!--       are of smaller size than the current particle and
387!--       use this radius for calculating the collision efficiency
388          DO  n = psi+prt_count(k,j,i)-1, psi+1, -1
389
390             sl_r3 = 0.0_wp
391             sl_r4 = 0.0_wp
392
393             DO is = n-1, psi, -1
394                IF ( particles(is)%radius < particles(n)%radius )  &
395                THEN
396                   sl_r3 = sl_r3 + particles(is)%weight_factor     &
397                                   * particles(is)%radius**3
398                   sl_r4 = sl_r4 + particles(is)%weight_factor     &
399                                   * particles(is)%radius**4
400                ENDIF
401             ENDDO
402
403             IF ( ( sl_r3 ) > 0.0_wp )  THEN
404                mean_r = ( sl_r4 ) / ( sl_r3 )
405
406                CALL collision_efficiency_rogers( mean_r,             &
407                                           particles(n)%radius,       &
408                                           effective_coll_efficiency )
409
410             ELSE
411                effective_coll_efficiency = 0.0_wp
412             ENDIF
413
414             IF ( effective_coll_efficiency > 1.0_wp  .OR.            &
415                  effective_coll_efficiency < 0.0_wp )                &
416             THEN
417                WRITE( message_string, * )  'collision_efficien' , &
418                          'cy out of range:' ,effective_coll_efficiency
419                CALL message( 'lpm_droplet_collision', 'PA0145', 2, &
420                              2, -1, 6, 1 )
421             ENDIF
422
423!
424!--          Interpolation of liquid water content
425             ii = particles(n)%x * ddx
426             jj = particles(n)%y * ddy
427             kk = ( particles(n)%z + 0.5_wp * dz ) / dz
428
429             x  = particles(n)%x - ii * dx
430             y  = particles(n)%y - jj * dy
431             aa = x**2          + y**2
432             bb = ( dx - x )**2 + y**2
433             cc = x**2          + ( dy - y )**2
434             dd = ( dx - x )**2 + ( dy - y )**2
435             gg = aa + bb + cc + dd
436
437             ql_int_l = ( (gg-aa) * ql(kk,jj,ii)   + (gg-bb) *     &
438                                                  ql(kk,jj,ii+1)   &
439                        + (gg-cc) * ql(kk,jj+1,ii) + ( gg-dd ) *   &
440                                                  ql(kk,jj+1,ii+1) &
441                        ) / ( 3.0_wp * gg )
442
443             ql_int_u = ( (gg-aa) * ql(kk+1,jj,ii)   + (gg-bb) *   &
444                                                ql(kk+1,jj,ii+1)   &
445                        + (gg-cc) * ql(kk+1,jj+1,ii) + (gg-dd) *   &
446                                                ql(kk+1,jj+1,ii+1) &
447                        ) / ( 3.0_wp * gg )
448
449             ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz *&
450                                 ( ql_int_u - ql_int_l )
451
452!
453!--          Interpolate u velocity-component
454             ii = ( particles(n)%x + 0.5_wp * dx ) * ddx
455             jj =   particles(n)%y * ddy
456             kk = ( particles(n)%z + 0.5_wp * dz ) / dz ! only if equidistant
457
458             IF ( ( particles(n)%z - zu(kk) ) > ( 0.5_wp * dz ) ) kk = kk+1
459
460             x  = particles(n)%x + ( 0.5_wp - ii ) * dx
461             y  = particles(n)%y - jj * dy
462             aa = x**2          + y**2
463             bb = ( dx - x )**2 + y**2
464             cc = x**2          + ( dy - y )**2
465             dd = ( dx - x )**2 + ( dy - y )**2
466             gg = aa + bb + cc + dd
467
468             u_int_l = ( (gg-aa) * u(kk,jj,ii)   + (gg-bb) *       &
469                                                  u(kk,jj,ii+1)    &
470                       + (gg-cc) * u(kk,jj+1,ii) + (gg-dd) *       &
471                                                  u(kk,jj+1,ii+1)  &
472                       ) / ( 3.0_wp * gg ) - u_gtrans
473             IF ( kk+1 == nzt+1 )  THEN
474                u_int = u_int_l
475             ELSE
476                u_int_u = ( (gg-aa) * u(kk+1,jj,ii)   + (gg-bb) *  &
477                                                 u(kk+1,jj,ii+1)   &
478                          + (gg-cc) * u(kk+1,jj+1,ii) + (gg-dd) *  &
479                                                 u(kk+1,jj+1,ii+1) &
480                          ) / ( 3.0_wp * gg ) - u_gtrans
481                u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz &
482                                  * ( u_int_u - u_int_l )
483             ENDIF
484
485!
486!--          Same procedure for interpolation of the v velocity-component
487!--          (adopt index k from u velocity-component)
488             ii =   particles(n)%x * ddx
489             jj = ( particles(n)%y + 0.5_wp * dy ) * ddy
490
491             x  = particles(n)%x - ii * dx
492             y  = particles(n)%y + ( 0.5_wp - jj ) * dy
493             aa = x**2          + y**2
494             bb = ( dx - x )**2 + y**2
495             cc = x**2          + ( dy - y )**2
496             dd = ( dx - x )**2 + ( dy - y )**2
497             gg = aa + bb + cc + dd
498
499             v_int_l = ( ( gg-aa ) * v(kk,jj,ii)   + ( gg-bb ) *   &
500                                                   v(kk,jj,ii+1)   &
501                       + ( gg-cc ) * v(kk,jj+1,ii) + ( gg-dd ) *   &
502                                                   v(kk,jj+1,ii+1) &
503                       ) / ( 3.0_wp * gg ) - v_gtrans
504             IF ( kk+1 == nzt+1 )  THEN
505                v_int = v_int_l
506             ELSE
507                v_int_u = ( (gg-aa) * v(kk+1,jj,ii)   + (gg-bb) *  &
508                                                   v(kk+1,jj,ii+1) &
509                          + (gg-cc) * v(kk+1,jj+1,ii) + (gg-dd) *  &
510                                                 v(kk+1,jj+1,ii+1) &
511                          ) / ( 3.0_wp * gg ) - v_gtrans
512                v_int = v_int_l + ( particles(n)%z - zu(kk) ) / dz &
513                                  * ( v_int_u - v_int_l )
514             ENDIF
515
516!
517!--          Same procedure for interpolation of the w velocity-component
518!--          (adopt index i from v velocity-component)
519             jj = particles(n)%y * ddy
520             kk = particles(n)%z / dz
521
522             x  = particles(n)%x - ii * dx
523             y  = particles(n)%y - jj * dy
524             aa = x**2          + y**2
525             bb = ( dx - x )**2 + y**2
526             cc = x**2          + ( dy - y )**2
527             dd = ( dx - x )**2 + ( dy - y )**2
528             gg = aa + bb + cc + dd
529
530             w_int_l = ( ( gg-aa ) * w(kk,jj,ii)   + ( gg-bb ) *   &
531                                                     w(kk,jj,ii+1) &
532                       + ( gg-cc ) * w(kk,jj+1,ii) + ( gg-dd ) *   &
533                                                   w(kk,jj+1,ii+1) &
534                       ) / ( 3.0_wp * gg )
535             IF ( kk+1 == nzt+1 )  THEN
536                w_int = w_int_l
537             ELSE
538                w_int_u = ( (gg-aa) * w(kk+1,jj,ii)   + (gg-bb) *  &
539                                                   w(kk+1,jj,ii+1) &
540                          + (gg-cc) * w(kk+1,jj+1,ii) + (gg-dd) *  &
541                                                 w(kk+1,jj+1,ii+1) &
542                          ) / ( 3.0_wp * gg )
543                w_int = w_int_l + ( particles(n)%z - zw(kk) ) / dz &
544                                  * ( w_int_u - w_int_l )
545             ENDIF
546
547!
548!--          Change in radius due to collision
549             delta_r = effective_coll_efficiency / 3.0_wp          &
550                       * pi * sl_r3 * ddx * ddy / dz               &
551                       * SQRT( ( u_int - particles(n)%speed_x )**2 &
552                             + ( v_int - particles(n)%speed_y )**2 &
553                             + ( w_int - particles(n)%speed_z )**2 &
554                             ) * dt_3d
555!
556!--          Change in volume due to collision
557             delta_v = particles(n)%weight_factor                  &
558                       * ( ( particles(n)%radius + delta_r )**3    &
559                           - particles(n)%radius**3 )
560
561!
562!--          Check if collected particles provide enough LWC for
563!--          volume change of collector particle
564             IF ( delta_v >= sl_r3  .AND.  sl_r3 > 0.0_wp )  THEN
565
566                delta_r = ( ( sl_r3/particles(n)%weight_factor )               &
567                            + particles(n)%radius**3 )**( 1.0_wp / 3.0_wp )    &
568                            - particles(n)%radius
569
570                DO  is = n-1, psi, -1
571                   IF ( particles(is)%radius < particles(n)%radius )  THEN
572                      particles(is)%weight_factor = 0.0_wp
573                      particles(is)%particle_mask = .FALSE.
574                      deleted_particles = deleted_particles + 1
575                   ENDIF
576                ENDDO
577
578             ELSE IF ( delta_v < sl_r3  .AND.  sl_r3 > 0.0_wp )  THEN
579
580                DO  is = n-1, psi, -1
581                   IF ( particles(is)%radius < particles(n)%radius &
582                        .AND.  sl_r3 > 0.0_wp )  THEN
583                      particles(is)%weight_factor =                &
584                                   ( ( particles(is)%weight_factor &
585                                   * ( particles(is)%radius**3 ) ) &
586                                   - ( delta_v                     &
587                                   * particles(is)%weight_factor   &
588                                   * ( particles(is)%radius**3 )   &
589                                   / sl_r3 ) )                     &
590                                   / ( particles(is)%radius**3 )
591
592                      IF ( particles(is)%weight_factor < 0.0_wp )  THEN
593                         WRITE( message_string, * ) 'negative ',   &
594                                         'weighting factor: ',     &
595                                         particles(is)%weight_factor
596                         CALL message( 'lpm_droplet_collision', &
597                                       'PA0039',                &
598                                       2, 2, -1, 6, 1 )
599                      ENDIF
600                   ENDIF
601                ENDDO
602
603             ENDIF
604
605             particles(n)%radius = particles(n)%radius + delta_r
606             ql_vp(k,j,i) = ql_vp(k,j,i) +               &
607                            particles(n)%weight_factor * &
608                            ( particles(n)%radius**3 )
609          ENDDO
610
611          ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor  &
612                                        * particles(psi)%radius**3
613
614       ENDIF  ! collision kernel
615
616    ELSE IF ( prt_count(k,j,i) == 1 )  THEN
617
618       psi = 1
619
620!
621!--    Calculate change of weighting factor due to self collision
622       IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
623            use_kernel_tables )  THEN
624
625          IF ( wang_kernel )  THEN
626             eclass = INT( diss(k,j,i) * 1.0E4_wp / 1000.0_wp * &
627                           dissipation_classes ) + 1
628             epsilon = diss(k,j,i)
629          ELSE
630             epsilon = 0.0_wp
631          ENDIF
632          IF ( hall_kernel  .OR.  epsilon * 1.0E4_wp < 0.001_wp )  THEN
633             eclass = 0   ! Hall kernel is used
634          ELSE
635             eclass = MIN( dissipation_classes, eclass )
636          ENDIF
637
638          ddV = ddx * ddy / dz
639          rclass_l = particles(psi)%class
640          sum3 = 1 - dt_3d * ddV *                                 &
641                     ( ckernel(rclass_l,rclass_l,eclass) *         &
642                       ( particles(psi)%weight_factor-1 ) * 0.5_wp )
643
644          particles(psi)%radius = ( particles(psi)%radius**3 / &
645                                    sum3 )**0.33333333333333_wp
646          particles(psi)%weight_factor = particles(psi)%weight_factor &
647                                         * sum3
648
649       ELSE IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
650                  .NOT. use_kernel_tables )  THEN
651!
652!--       Collision efficiencies are calculated for every new
653!--       grid box. First, allocate memory for kernel table.
654!--       Third dimension is 1, because table is re-calculated for
655!--       every new dissipation value.
656          ALLOCATE( ckernel(psi:psi, psi:psi, 1:1) )
657!
658!--       Now calculate collision efficiencies for this box
659          CALL recalculate_kernel( i, j, k )
660
661          ddV = ddx * ddy / dz
662          sum3 = 1 - dt_3d * ddV * ( ckernel(psi,psi,1) *  &
663                       ( particles(psi)%weight_factor - 1 ) * 0.5_wp ) 
664
665          particles(psi)%radius = ( particles(psi)%radius**3 / &
666                                    sum3 )**0.33333333333333_wp
667          particles(psi)%weight_factor = particles(psi)%weight_factor &
668                                         * sum3
669
670          DEALLOCATE( ckernel )
671       ENDIF
672
673      ql_vp(k,j,i) =  particles(psi)%weight_factor *              &
674                       particles(psi)%radius**3
675    ENDIF
676
677!
678!-- Check if condensation of LWC was conserved during collision process
679    IF ( ql_v(k,j,i) /= 0.0_wp )  THEN
680       IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001_wp  .OR.             &
681            ql_vp(k,j,i) / ql_v(k,j,i) <= 0.9999_wp )  THEN
682          WRITE( message_string, * ) 'LWC is not conserved during',&
683                                     ' collision! ',               &
684                                     'LWC after condensation: ',   &
685                                     ql_v(k,j,i),                  &
686                                     ' LWC after collision: ',     &
687                                     ql_vp(k,j,i)
688          CALL message( 'lpm_droplet_collision', 'PA0040',         &
689                        2, 2, -1, 6, 1 )
690       ENDIF
691    ENDIF
692
693    CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'stop' )
694
695 END SUBROUTINE lpm_droplet_collision
Note: See TracBrowser for help on using the repository browser.