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

Last change on this file since 1071 was 1071, checked in by franke, 11 years ago

ventilation effect for cloud droplets included, calculation of cloud droplet collisions now uses formulation by Wang, bugfixes in Rosenbrock method and in calculation of collision efficiencies

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