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

Last change on this file since 850 was 850, checked in by raasch, 9 years ago

last commit documented

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