source: palm/tags/release-4.0/SOURCE/lpm_droplet_collision.f90 @ 1983

Last change on this file since 1983 was 1360, checked in by hoffmann, 10 years ago

last commit documented

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