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

Last change on this file since 1046 was 1037, checked in by raasch, 11 years ago

last commit documented

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