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

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

last commit documented

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