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

Last change on this file since 849 was 849, checked in by raasch, 10 years ago

Changed:


Original routine advec_particles split into several new subroutines and renamed
lpm.
init_particles renamed lpm_init
user_advec_particles renamed user_lpm_advec,
particle_boundary_conds renamed lpm_boundary_conds,
set_particle_attributes renamed lpm_set_attributes,
user_init_particles renamed user_lpm_init,
user_particle_attributes renamed user_lpm_set_attributes
(Makefile, lpm_droplet_collision, lpm_droplet_condensation, init_3d_model, modules, palm, read_var_list, time_integration, write_var_list, deleted: advec_particles, init_particles, particle_boundary_conds, set_particle_attributes, user_advec_particles, user_init_particles, user_particle_attributes, new: lpm, lpm_advec, lpm_boundary_conds, lpm_calc_liquid_water_content, lpm_data_output_particles, lpm_droplet_collision, lpm_drollet_condensation, lpm_exchange_horiz, lpm_extend_particle_array, lpm_extend_tails, lpm_extend_tail_array, lpm_init, lpm_init_sgs_tke, lpm_pack_arrays, lpm_read_restart_file, lpm_release_set, lpm_set_attributes, lpm_sort_arrays, lpm_write_exchange_statistics, lpm_write_restart_file, user_lpm_advec, user_lpm_init, user_lpm_set_attributes

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