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

Last change on this file since 1818 was 1818, checked in by maronga, 5 years ago

last commit documented / copyright update

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