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

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

former files/routines cpu_log and cpu_statistics combined to one module,
which also includes the former data module cpulog from the modules-file,
module interfaces removed

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