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

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

last commit documented

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