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

Last change on this file since 1359 was 1359, checked in by hoffmann, 7 years ago

new Lagrangian particle structure integrated

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