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

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

update of GPL copyright

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