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

Last change on this file since 1323 was 1323, checked in by raasch, 8 years ago

last commit documented

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