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

Last change on this file since 1092 was 1092, checked in by raasch, 9 years ago

unused variables remove from several routines

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