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

Last change on this file since 1320 was 1320, checked in by raasch, 7 years ago

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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