Ignore:
Timestamp:
Apr 11, 2014 5:15:14 PM (7 years ago)
Author:
hoffmann
Message:

new Lagrangian particle structure integrated

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/lpm_droplet_collision.f90

    r1323 r1359  
    1  SUBROUTINE lpm_droplet_collision
     1 SUBROUTINE lpm_droplet_collision (i,j,k)
    22
    33!--------------------------------------------------------------------------------!
     
    2020! Current revisions:
    2121! ------------------
    22 !
     22! New particle structure integrated.
     23! Kind definition added to all floating point numbers.
    2324!
    2425! Former revisions:
     
    7576!------------------------------------------------------------------------------!
    7677
     78
    7779    USE arrays_3d,                                                             &
    7880        ONLY:  diss, ql, ql_v, ql_vp, u, v, w, zu, zw
     
    103105    USE particle_attributes,                                                   &
    104106        ONLY:  deleted_particles, dissipation_classes, hall_kernel,            &
    105                palm_kernel, particles, particle_mask, particle_type,           &
    106                prt_count, prt_start_index, use_kernel_tables, wang_kernel
     107               palm_kernel, particles, particle_type,                          &
     108               prt_count, use_kernel_tables, wang_kernel
     109
     110    USE pegrid
    107111
    108112    IMPLICIT NONE
     
    124128    INTEGER(iwp) ::  rclass_s !:
    125129
     130    INTEGER(iwp), DIMENSION(prt_count(k,j,i)) ::  rclass_v !:
     131
     132    LOGICAL, SAVE ::  first_flag = .TRUE. !:
     133
     134    TYPE(particle_type) :: tmp_particle   !:
     135
    126136    REAL(wp) ::  aa       !:
     137    REAL(wp) ::  auxn     !: temporary variables
     138    REAL(wp) ::  auxs     !: temporary variables
    127139    REAL(wp) ::  bb       !:
    128140    REAL(wp) ::  cc       !:
     
    158170    REAL(wp), DIMENSION(:), ALLOCATABLE ::  weight !:
    159171
    160 
    161     TYPE(particle_type) ::  tmp_particle           !:
    162 
    163 
     172    REAL, DIMENSION(prt_count(k,j,i))    :: ck
     173    REAL, DIMENSION(prt_count(k,j,i))    :: r3v
     174    REAL, DIMENSION(prt_count(k,j,i))    :: sum1v
     175    REAL, DIMENSION(prt_count(k,j,i))    :: sum2v
    164176
    165177    CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'start' )
    166178
    167     DO  i = nxl, nxr
    168        DO  j = nys, nyn
    169           DO  k = nzb+1, nzt
    170 !
    171 !--          Collision requires at least two particles in the box
    172              IF ( prt_count(k,j,i) > 1 )  THEN
    173 !
    174 !--             First, sort particles within the gridbox by their size,
    175 !--             using Shell's method (see Numerical Recipes)
    176 !--             NOTE: In case of using particle tails, the re-sorting of
    177 !--             ----  tails would have to be included here!
    178                 psi = prt_start_index(k,j,i) - 1
    179                 inc = 1
    180                 DO WHILE ( inc <= prt_count(k,j,i) )
    181                    inc = 3 * inc + 1
     179!
     180!-- Collision requires at least two particles in the box
     181    IF ( prt_count(k,j,i) > 1 )  THEN
     182!
     183!--    First, sort particles within the gridbox by their size,
     184!--    using Shell's method (see Numerical Recipes)
     185!--    NOTE: In case of using particle tails, the re-sorting of
     186!--    ----  tails would have to be included here!
     187       IF ( .NOT. ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
     188       use_kernel_tables ) )  THEN
     189          psi = 0
     190          inc = 1
     191          DO WHILE ( inc <= prt_count(k,j,i) )
     192             inc = 3 * inc + 1
     193          ENDDO
     194
     195          DO WHILE ( inc > 1 )
     196             inc = inc / 3
     197             DO  is = inc+1, prt_count(k,j,i)
     198                tmp_particle = particles(psi+is)
     199                js = is
     200                DO WHILE ( particles(psi+js-inc)%radius >             &
     201                tmp_particle%radius )
     202                   particles(psi+js) = particles(psi+js-inc)
     203                   js = js - inc
     204                   IF ( js <= inc )  EXIT
    182205                ENDDO
    183 
    184                 DO WHILE ( inc > 1 )
    185                    inc = inc / 3
    186                    DO  is = inc+1, prt_count(k,j,i)
    187                       tmp_particle = particles(psi+is)
    188                       js = is
    189                       DO WHILE ( particles(psi+js-inc)%radius >             &
    190                                  tmp_particle%radius )
    191                          particles(psi+js) = particles(psi+js-inc)
    192                          js = js - inc
    193                          IF ( js <= inc )  EXIT
    194                       ENDDO
    195                       particles(psi+js) = tmp_particle
    196                    ENDDO
     206                particles(psi+js) = tmp_particle
     207             ENDDO
     208          ENDDO
     209       ENDIF
     210
     211       psi = 1
     212       pse = prt_count(k,j,i)
     213
     214!
     215!--    Now apply the different kernels
     216       IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
     217            use_kernel_tables )  THEN
     218!
     219!--       Fast method with pre-calculated efficiencies for
     220!--       discrete radius- and dissipation-classes.
     221!--
     222!--       Determine dissipation class index of this gridbox
     223          IF ( wang_kernel )  THEN
     224             eclass = INT( diss(k,j,i) * 1.0E4_wp / 1000.0_wp * &
     225                           dissipation_classes ) + 1
     226             epsilon = diss(k,j,i)
     227          ELSE
     228             epsilon = 0.0_wp
     229          ENDIF
     230          IF ( hall_kernel  .OR.  epsilon * 1.0E4_wp < 0.001_wp )  THEN
     231             eclass = 0   ! Hall kernel is used
     232          ELSE
     233             eclass = MIN( dissipation_classes, eclass )
     234          ENDIF
     235
     236!
     237!--       Droplet collision are calculated using collision-coalescence
     238!--       formulation proposed by Wang (see PALM documentation)
     239!--       Since new radii after collision are defined by radii of all
     240!--       droplets before collision, temporary fields for new radii and
     241!--       weighting factors are needed
     242          ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i)))
     243
     244          rad    = 0.0_wp
     245          weight = 0.0_wp
     246
     247          sum1v(1:prt_count(k,j,i)) = 0.0_wp
     248          sum2v(1:prt_count(k,j,i)) = 0.0_wp
     249
     250          DO  n = 1, prt_count(k,j,i)
     251
     252             rclass_l = particles(n)%class
     253!
     254!--          Mass added due to collisions with smaller droplets
     255             DO  is = n+1, prt_count(k,j,i)
     256                rclass_s = particles(is)%class
     257                auxs = ckernel(rclass_l,rclass_s,eclass) * particles(is)%weight_factor
     258                auxn = ckernel(rclass_s,rclass_l,eclass) * particles(n)%weight_factor
     259                IF ( particles(is)%radius <  particles(n)%radius )  THEN
     260                   sum1v(n) =  sum1v(n)  + particles(is)%radius**3 * auxs
     261                   sum2v(is) = sum2v(is) + auxn
     262                ELSE
     263                   sum2v(n)  = sum2v(n)  + auxs
     264                   sum1v(is) = sum1v(is) + particles(n)%radius**3 * auxn
     265                ENDIF
     266             ENDDO
     267          ENDDO
     268          rclass_v = particles(1:prt_count(k,j,i))%class
     269          DO  n = 1, prt_count(k,j,i)
     270            ck(n)       = ckernel(rclass_v(n),rclass_v(n),eclass)
     271          ENDDO
     272          r3v      = particles(1:prt_count(k,j,i))%radius**3
     273          DO  n = 1, prt_count(k,j,i)
     274             sum3 = 0.0_wp
     275             ddV = ddx * ddy / dz
     276!
     277!--          Change of the current weighting factor
     278             sum3 = 1 - dt_3d * ddV *                                 &
     279                        ck(n) *                                       &
     280                        ( particles(n)%weight_factor - 1 ) * 0.5_wp - &
     281                    dt_3d * ddV * sum2v(n)
     282             weight(n) = particles(n)%weight_factor * sum3
     283!
     284!--          Change of the current droplet radius
     285             rad(n) = ( (r3v(n) + dt_3d * ddV * (sum1v(n) - sum2v(n) * r3v(n)) )/&
     286                              sum3 )**0.33333333333333_wp
     287
     288             ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n) &
     289                                         * rad(n)**3
     290
     291          ENDDO
     292          IF ( ANY(weight < 0.0_wp) )  THEN
     293                WRITE( message_string, * ) 'negative weighting'
     294                CALL message( 'lpm_droplet_collision', 'PA0028',      &
     295                               2, 2, -1, 6, 1 )
     296          ENDIF
     297
     298          particles(psi:pse)%radius = rad(1:prt_count(k,j,i))
     299          particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i))
     300
     301          DEALLOCATE(rad, weight)
     302
     303       ELSEIF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
     304                .NOT. use_kernel_tables )  THEN
     305!
     306!--       Collision efficiencies are calculated for every new
     307!--       grid box. First, allocate memory for kernel table.
     308!--       Third dimension is 1, because table is re-calculated for
     309!--       every new dissipation value.
     310          ALLOCATE( ckernel(1:prt_count(k,j,i),1:prt_count(k,j,i),1:1) )
     311!
     312!--       Now calculate collision efficiencies for this box
     313          CALL recalculate_kernel( i, j, k )
     314
     315!
     316!--       Droplet collision are calculated using collision-coalescence
     317!--       formulation proposed by Wang (see PALM documentation)
     318!--       Since new radii after collision are defined by radii of all
     319!--       droplets before collision, temporary fields for new radii and
     320!--       weighting factors are needed
     321          ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i)))
     322
     323          rad = 0.0_wp
     324          weight = 0.0_wp
     325
     326          DO  n = psi, pse, 1
     327
     328             sum1 = 0.0_wp
     329             sum2 = 0.0_wp
     330             sum3 = 0.0_wp
     331!
     332!--          Mass added due to collisions with smaller droplets
     333             DO  is = psi, n-1
     334                sum1 = sum1 + ( particles(is)%radius**3 *            &
     335                                ckernel(n,is,1) *  &
     336                                particles(is)%weight_factor )
     337             ENDDO
     338!
     339!--          Rate of collisions with larger droplets
     340             DO  is = n+1, pse
     341                sum2 = sum2 + ( ckernel(n,is,1) *  &
     342                                particles(is)%weight_factor )
     343             ENDDO
     344
     345             r3 = particles(n)%radius**3
     346             ddV = ddx * ddy / dz
     347             is = 1
     348!
     349!--                   Change of the current weighting factor
     350             sum3 = 1 - dt_3d * ddV *                                 &
     351                        ckernel(n,n,1) *           &
     352                        ( particles(n)%weight_factor - 1 ) * 0.5_wp - &
     353                    dt_3d * ddV * sum2
     354             weight(n-is+1) = particles(n)%weight_factor * sum3
     355!
     356!--                   Change of the current droplet radius
     357             rad(n-is+1) = ( (r3 + dt_3d * ddV * (sum1 - sum2 * r3) )/&
     358                              sum3 )**0.33333333333333_wp
     359
     360             IF ( weight(n-is+1) < 0.0_wp )  THEN
     361                WRITE( message_string, * ) 'negative weighting',      &
     362                                           'factor: ', weight(n-is+1)
     363                CALL message( 'lpm_droplet_collision', 'PA0037',      &
     364                               2, 2, -1, 6, 1 )
     365             ENDIF
     366
     367             ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n-is+1) &
     368                                         * rad(n-is+1)**3
     369
     370          ENDDO
     371
     372          particles(psi:pse)%radius = rad(1:prt_count(k,j,i))
     373          particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i))
     374
     375          DEALLOCATE( rad, weight, ckernel )
     376
     377       ELSEIF ( palm_kernel )  THEN
     378!
     379!--       PALM collision kernel
     380!
     381!--       Calculate the mean radius of all those particles which
     382!--       are of smaller size than the current particle and
     383!--       use this radius for calculating the collision efficiency
     384          DO  n = psi+prt_count(k,j,i)-1, psi+1, -1
     385
     386             sl_r3 = 0.0_wp
     387             sl_r4 = 0.0_wp
     388
     389             DO is = n-1, psi, -1
     390                IF ( particles(is)%radius < particles(n)%radius )  &
     391                THEN
     392                   sl_r3 = sl_r3 + particles(is)%weight_factor     &
     393                                   * particles(is)%radius**3
     394                   sl_r4 = sl_r4 + particles(is)%weight_factor     &
     395                                   * particles(is)%radius**4
     396                ENDIF
     397             ENDDO
     398
     399             IF ( ( sl_r3 ) > 0.0_wp )  THEN
     400                mean_r = ( sl_r4 ) / ( sl_r3 )
     401
     402                CALL collision_efficiency_rogers( mean_r,             &
     403                                           particles(n)%radius,       &
     404                                           effective_coll_efficiency )
     405
     406             ELSE
     407                effective_coll_efficiency = 0.0_wp
     408             ENDIF
     409
     410             IF ( effective_coll_efficiency > 1.0_wp  .OR.            &
     411                  effective_coll_efficiency < 0.0_wp )                &
     412             THEN
     413                WRITE( message_string, * )  'collision_efficien' , &
     414                          'cy out of range:' ,effective_coll_efficiency
     415                CALL message( 'lpm_droplet_collision', 'PA0145', 2, &
     416                              2, -1, 6, 1 )
     417             ENDIF
     418
     419!
     420!--          Interpolation of liquid water content
     421             ii = particles(n)%x * ddx
     422             jj = particles(n)%y * ddy
     423             kk = ( particles(n)%z + 0.5_wp * dz ) / dz
     424
     425             x  = particles(n)%x - ii * dx
     426             y  = particles(n)%y - jj * dy
     427             aa = x**2          + y**2
     428             bb = ( dx - x )**2 + y**2
     429             cc = x**2          + ( dy - y )**2
     430             dd = ( dx - x )**2 + ( dy - y )**2
     431             gg = aa + bb + cc + dd
     432
     433             ql_int_l = ( (gg-aa) * ql(kk,jj,ii)   + (gg-bb) *     &
     434                                                  ql(kk,jj,ii+1)   &
     435                        + (gg-cc) * ql(kk,jj+1,ii) + ( gg-dd ) *   &
     436                                                  ql(kk,jj+1,ii+1) &
     437                        ) / ( 3.0_wp * gg )
     438
     439             ql_int_u = ( (gg-aa) * ql(kk+1,jj,ii)   + (gg-bb) *   &
     440                                                ql(kk+1,jj,ii+1)   &
     441                        + (gg-cc) * ql(kk+1,jj+1,ii) + (gg-dd) *   &
     442                                                ql(kk+1,jj+1,ii+1) &
     443                        ) / ( 3.0_wp * gg )
     444
     445             ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz *&
     446                                 ( ql_int_u - ql_int_l )
     447
     448!
     449!--          Interpolate u velocity-component
     450             ii = ( particles(n)%x + 0.5_wp * dx ) * ddx
     451             jj =   particles(n)%y * ddy
     452             kk = ( particles(n)%z + 0.5_wp * dz ) / dz ! only if equidistant
     453
     454             IF ( ( particles(n)%z - zu(kk) ) > ( 0.5_wp * dz ) ) kk = kk+1
     455
     456             x  = particles(n)%x + ( 0.5_wp - ii ) * dx
     457             y  = particles(n)%y - jj * dy
     458             aa = x**2          + y**2
     459             bb = ( dx - x )**2 + y**2
     460             cc = x**2          + ( dy - y )**2
     461             dd = ( dx - x )**2 + ( dy - y )**2
     462             gg = aa + bb + cc + dd
     463
     464             u_int_l = ( (gg-aa) * u(kk,jj,ii)   + (gg-bb) *       &
     465                                                  u(kk,jj,ii+1)    &
     466                       + (gg-cc) * u(kk,jj+1,ii) + (gg-dd) *       &
     467                                                  u(kk,jj+1,ii+1)  &
     468                       ) / ( 3.0_wp * gg ) - u_gtrans
     469             IF ( kk+1 == nzt+1 )  THEN
     470                u_int = u_int_l
     471             ELSE
     472                u_int_u = ( (gg-aa) * u(kk+1,jj,ii)   + (gg-bb) *  &
     473                                                 u(kk+1,jj,ii+1)   &
     474                          + (gg-cc) * u(kk+1,jj+1,ii) + (gg-dd) *  &
     475                                                 u(kk+1,jj+1,ii+1) &
     476                          ) / ( 3.0_wp * gg ) - u_gtrans
     477                u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz &
     478                                  * ( u_int_u - u_int_l )
     479             ENDIF
     480
     481!
     482!--          Same procedure for interpolation of the v velocity-component
     483!--          (adopt index k from u velocity-component)
     484             ii =   particles(n)%x * ddx
     485             jj = ( particles(n)%y + 0.5_wp * dy ) * ddy
     486
     487             x  = particles(n)%x - ii * dx
     488             y  = particles(n)%y + ( 0.5_wp - jj ) * dy
     489             aa = x**2          + y**2
     490             bb = ( dx - x )**2 + y**2
     491             cc = x**2          + ( dy - y )**2
     492             dd = ( dx - x )**2 + ( dy - y )**2
     493             gg = aa + bb + cc + dd
     494
     495             v_int_l = ( ( gg-aa ) * v(kk,jj,ii)   + ( gg-bb ) *   &
     496                                                   v(kk,jj,ii+1)   &
     497                       + ( gg-cc ) * v(kk,jj+1,ii) + ( gg-dd ) *   &
     498                                                   v(kk,jj+1,ii+1) &
     499                       ) / ( 3.0_wp * gg ) - v_gtrans
     500             IF ( kk+1 == nzt+1 )  THEN
     501                v_int = v_int_l
     502             ELSE
     503                v_int_u = ( (gg-aa) * v(kk+1,jj,ii)   + (gg-bb) *  &
     504                                                   v(kk+1,jj,ii+1) &
     505                          + (gg-cc) * v(kk+1,jj+1,ii) + (gg-dd) *  &
     506                                                 v(kk+1,jj+1,ii+1) &
     507                          ) / ( 3.0_wp * gg ) - v_gtrans
     508                v_int = v_int_l + ( particles(n)%z - zu(kk) ) / dz &
     509                                  * ( v_int_u - v_int_l )
     510             ENDIF
     511
     512!
     513!--          Same procedure for interpolation of the w velocity-component
     514!--          (adopt index i from v velocity-component)
     515             jj = particles(n)%y * ddy
     516             kk = particles(n)%z / dz
     517
     518             x  = particles(n)%x - ii * dx
     519             y  = particles(n)%y - jj * dy
     520             aa = x**2          + y**2
     521             bb = ( dx - x )**2 + y**2
     522             cc = x**2          + ( dy - y )**2
     523             dd = ( dx - x )**2 + ( dy - y )**2
     524             gg = aa + bb + cc + dd
     525
     526             w_int_l = ( ( gg-aa ) * w(kk,jj,ii)   + ( gg-bb ) *   &
     527                                                     w(kk,jj,ii+1) &
     528                       + ( gg-cc ) * w(kk,jj+1,ii) + ( gg-dd ) *   &
     529                                                   w(kk,jj+1,ii+1) &
     530                       ) / ( 3.0_wp * gg )
     531             IF ( kk+1 == nzt+1 )  THEN
     532                w_int = w_int_l
     533             ELSE
     534                w_int_u = ( (gg-aa) * w(kk+1,jj,ii)   + (gg-bb) *  &
     535                                                   w(kk+1,jj,ii+1) &
     536                          + (gg-cc) * w(kk+1,jj+1,ii) + (gg-dd) *  &
     537                                                 w(kk+1,jj+1,ii+1) &
     538                          ) / ( 3.0_wp * gg )
     539                w_int = w_int_l + ( particles(n)%z - zw(kk) ) / dz &
     540                                  * ( w_int_u - w_int_l )
     541             ENDIF
     542
     543!
     544!--          Change in radius due to collision
     545             delta_r = effective_coll_efficiency / 3.0_wp          &
     546                       * pi * sl_r3 * ddx * ddy / dz               &
     547                       * SQRT( ( u_int - particles(n)%speed_x )**2 &
     548                             + ( v_int - particles(n)%speed_y )**2 &
     549                             + ( w_int - particles(n)%speed_z )**2 &
     550                             ) * dt_3d
     551!
     552!--          Change in volume due to collision
     553             delta_v = particles(n)%weight_factor                  &
     554                       * ( ( particles(n)%radius + delta_r )**3    &
     555                           - particles(n)%radius**3 )
     556
     557!
     558!--          Check if collected particles provide enough LWC for
     559!--          volume change of collector particle
     560             IF ( delta_v >= sl_r3  .AND.  sl_r3 > 0.0_wp )  THEN
     561
     562                delta_r = ( ( sl_r3/particles(n)%weight_factor )               &
     563                            + particles(n)%radius**3 )**( 1.0_wp / 3.0_wp )    &
     564                            - particles(n)%radius
     565
     566                DO  is = n-1, psi, -1
     567                   IF ( particles(is)%radius < particles(n)%radius )  THEN
     568                      particles(is)%weight_factor = 0.0_wp
     569                      particles(is)%particle_mask = .FALSE.
     570                      deleted_particles = deleted_particles + 1
     571                   ENDIF
    197572                ENDDO
    198573
    199                 psi = prt_start_index(k,j,i)
    200                 pse = psi + prt_count(k,j,i)-1
    201 
    202 !
    203 !--             Now apply the different kernels
    204                 IF ( ( hall_kernel .OR. wang_kernel )  .AND.  &
    205                      use_kernel_tables )  THEN
    206 !
    207 !--                Fast method with pre-calculated efficiencies for
    208 !--                discrete radius- and dissipation-classes.
    209 !
    210 !--                Determine dissipation class index of this gridbox
    211                    IF ( wang_kernel )  THEN
    212                       eclass = INT( diss(k,j,i) * 1.0E4_wp / 1000.0_wp * &
    213                                     dissipation_classes ) + 1
    214                       epsilon = diss(k,j,i)
    215                    ELSE
    216                       epsilon = 0.0
     574             ELSE IF ( delta_v < sl_r3  .AND.  sl_r3 > 0.0_wp )  THEN
     575
     576                DO  is = n-1, psi, -1
     577                   IF ( particles(is)%radius < particles(n)%radius &
     578                        .AND.  sl_r3 > 0.0_wp )  THEN
     579                      particles(is)%weight_factor =                &
     580                                   ( ( particles(is)%weight_factor &
     581                                   * ( particles(is)%radius**3 ) ) &
     582                                   - ( delta_v                     &
     583                                   * particles(is)%weight_factor   &
     584                                   * ( particles(is)%radius**3 )   &
     585                                   / sl_r3 ) )                     &
     586                                   / ( particles(is)%radius**3 )
     587
     588                      IF ( particles(is)%weight_factor < 0.0_wp )  THEN
     589                         WRITE( message_string, * ) 'negative ',   &
     590                                         'weighting factor: ',     &
     591                                         particles(is)%weight_factor
     592                         CALL message( 'lpm_droplet_collision', &
     593                                       'PA0039',                &
     594                                       2, 2, -1, 6, 1 )
     595                      ENDIF
    217596                   ENDIF
    218                    IF ( hall_kernel .OR. epsilon * 1.0E4_wp < 0.001 )  THEN
    219                       eclass = 0   ! Hall kernel is used
    220                    ELSE
    221                       eclass = MIN( dissipation_classes, eclass )
    222                    ENDIF
    223 
    224  !
    225 !--                Droplet collision are calculated using collision-coalescence
    226 !--                formulation proposed by Wang (see PALM documentation)
    227 !--                Since new radii after collision are defined by radii of all
    228 !--                droplets before collision, temporary fields for new radii and
    229 !--                weighting factors are needed
    230                    ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i)))
    231 
    232                    rad = 0.0
    233                    weight = 0.0
    234 
    235                    DO  n = psi, pse, 1
    236 
    237                       sum1 = 0.0
    238                       sum2 = 0.0
    239                       sum3 = 0.0
    240 
    241                       rclass_l = particles(n)%class
    242 !
    243 !--                   Mass added due to collisions with smaller droplets
    244                       DO  is = psi, n-1
    245 
    246                          rclass_s = particles(is)%class
    247                          sum1 = sum1 + ( particles(is)%radius**3 *            &
    248                                          ckernel(rclass_l,rclass_s,eclass) *  &
    249                                          particles(is)%weight_factor )
    250 
    251                       ENDDO
    252 !
    253 !--                   Rate of collisions with larger droplets
    254                       DO  is = n+1, pse
    255 
    256                          rclass_s = particles(is)%class
    257                          sum2 = sum2 + ( ckernel(rclass_l,rclass_s,eclass) *  &
    258                                          particles(is)%weight_factor )
    259 
    260                       ENDDO
    261 
    262                       r3 = particles(n)%radius**3
    263                       ddV = ddx * ddy / dz
    264                       is = prt_start_index(k,j,i)
    265 !
    266 !--                   Change of the current weighting factor
    267                       sum3 = 1 - dt_3d * ddV *                                 &
    268                                  ckernel(rclass_l,rclass_l,eclass) *           &
    269                                  ( particles(n)%weight_factor - 1 ) * 0.5 -    &
    270                              dt_3d * ddV * sum2
    271                       weight(n-is+1) = particles(n)%weight_factor * sum3
    272 !
    273 !--                   Change of the current droplet radius
    274                       rad(n-is+1) = ( (r3 + dt_3d * ddV * (sum1 - sum2 * r3) )/&
    275                                        sum3 )**0.33333333333333_wp
    276 
    277                       IF ( weight(n-is+1) < 0.0 )  THEN
    278                          WRITE( message_string, * ) 'negative weighting',      &
    279                                                     'factor: ', weight(n-is+1)
    280                          CALL message( 'lpm_droplet_collision', 'PA0028',      &
    281                                         2, 2, -1, 6, 1 )
    282                       ENDIF
    283 
    284                       ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n-is+1) &
    285                                                   * rad(n-is+1)**3
    286 
    287                    ENDDO
    288 
    289                    particles(psi:pse)%radius = rad(1:prt_count(k,j,i))
    290                    particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i))
    291 
    292                    DEALLOCATE(rad, weight)
    293 
    294                 ELSEIF ( ( hall_kernel .OR. wang_kernel )  .AND.  &
    295                          .NOT. use_kernel_tables )  THEN
    296 !
    297 !--                Collision efficiencies are calculated for every new
    298 !--                grid box. First, allocate memory for kernel table.
    299 !--                Third dimension is 1, because table is re-calculated for
    300 !--                every new dissipation value.
    301                    ALLOCATE( ckernel(prt_start_index(k,j,i):                 &
    302                              prt_start_index(k,j,i)+prt_count(k,j,i)-1,      &
    303                              prt_start_index(k,j,i):                         &
    304                              prt_start_index(k,j,i)+prt_count(k,j,i)-1,1:1) )
    305 !
    306 !--                Now calculate collision efficiencies for this box
    307                    CALL recalculate_kernel( i, j, k )
    308 
    309 !
    310 !--                Droplet collision are calculated using collision-coalescence
    311 !--                formulation proposed by Wang (see PALM documentation)
    312 !--                Since new radii after collision are defined by radii of all
    313 !--                droplets before collision, temporary fields for new radii and
    314 !--                weighting factors are needed
    315                    ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i)))
    316 
    317                    rad = 0.0
    318                    weight = 0.0
    319 
    320                    DO  n = psi, pse, 1
    321 
    322                       sum1 = 0.0
    323                       sum2 = 0.0
    324                       sum3 = 0.0
    325 !
    326 !--                   Mass added due to collisions with smaller droplets
    327                       DO  is = psi, n-1
    328                          sum1 = sum1 + ( particles(is)%radius**3 *            &
    329                                          ckernel(n,is,1) *  &
    330                                          particles(is)%weight_factor )
    331                       ENDDO
    332 !
    333 !--                   Rate of collisions with larger droplets
    334                       DO  is = n+1, pse
    335                          sum2 = sum2 + ( ckernel(n,is,1) *  &
    336                                          particles(is)%weight_factor )
    337                       ENDDO
    338 
    339                       r3 = particles(n)%radius**3
    340                       ddV = ddx * ddy / dz
    341                       is = prt_start_index(k,j,i)
    342 !
    343 !--                   Change of the current weighting factor
    344                       sum3 = 1 - dt_3d * ddV *                                 &
    345                                  ckernel(n,n,1) *           &
    346                                  ( particles(n)%weight_factor - 1 ) * 0.5 -    &
    347                              dt_3d * ddV * sum2
    348                       weight(n-is+1) = particles(n)%weight_factor * sum3
    349 !
    350 !--                   Change of the current droplet radius
    351                       rad(n-is+1) = ( (r3 + dt_3d * ddV * (sum1 - sum2 * r3) )/&
    352                                        sum3 )**0.33333333333333_wp
    353 
    354                       IF ( weight(n-is+1) < 0.0 )  THEN
    355                          WRITE( message_string, * ) 'negative weighting',      &
    356                                                     'factor: ', weight(n-is+1)
    357                          CALL message( 'lpm_droplet_collision', 'PA0037',      &
    358                                         2, 2, -1, 6, 1 )
    359                       ENDIF
    360 
    361                       ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n-is+1) &
    362                                                   * rad(n-is+1)**3
    363 
    364                    ENDDO
    365 
    366                    particles(psi:pse)%radius = rad(1:prt_count(k,j,i))
    367                    particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i))
    368 
    369                    DEALLOCATE( rad, weight, ckernel )
    370 
    371                 ELSEIF ( palm_kernel )  THEN
    372 !
    373 !--                PALM collision kernel
    374 !
    375 !--                Calculate the mean radius of all those particles which
    376 !--                are of smaller size than the current particle and
    377 !--                use this radius for calculating the collision efficiency
    378                    DO  n = psi+prt_count(k,j,i)-1, psi+1, -1
    379 
    380                       sl_r3 = 0.0
    381                       sl_r4 = 0.0
    382 
    383                       DO is = n-1, psi, -1
    384                          IF ( particles(is)%radius < particles(n)%radius )  &
    385                          THEN
    386                             sl_r3 = sl_r3 + particles(is)%weight_factor     &
    387                                             * particles(is)%radius**3
    388                             sl_r4 = sl_r4 + particles(is)%weight_factor     &
    389                                             * particles(is)%radius**4
    390                          ENDIF
    391                       ENDDO
    392 
    393                       IF ( ( sl_r3 ) > 0.0 )  THEN
    394                          mean_r = ( sl_r4 ) / ( sl_r3 )
    395 
    396                          CALL collision_efficiency_rogers( mean_r,             &
    397                                                     particles(n)%radius,       &
    398                                                     effective_coll_efficiency )
    399 
    400                       ELSE
    401                          effective_coll_efficiency = 0.0
    402                       ENDIF
    403 
    404                       IF ( effective_coll_efficiency > 1.0  .OR.            &
    405                            effective_coll_efficiency < 0.0 )                &
    406                       THEN
    407                          WRITE( message_string, * )  'collision_efficien' , &
    408                                    'cy out of range:' ,effective_coll_efficiency
    409                          CALL message( 'lpm_droplet_collision', 'PA0145', 2, &
    410                                        2, -1, 6, 1 )
    411                       ENDIF
    412 
    413 !
    414 !--                   Interpolation of ...
    415                       ii = particles(n)%x * ddx
    416                       jj = particles(n)%y * ddy
    417                       kk = ( particles(n)%z + 0.5 * dz ) / dz
    418 
    419                       x  = particles(n)%x - ii * dx
    420                       y  = particles(n)%y - jj * dy
    421                       aa = x**2          + y**2
    422                       bb = ( dx - x )**2 + y**2
    423                       cc = x**2          + ( dy - y )**2
    424                       dd = ( dx - x )**2 + ( dy - y )**2
    425                       gg = aa + bb + cc + dd
    426 
    427                       ql_int_l = ( (gg-aa) * ql(kk,jj,ii)   + (gg-bb) *     &
    428                                                            ql(kk,jj,ii+1)   &
    429                                  + (gg-cc) * ql(kk,jj+1,ii) + ( gg-dd ) *   &
    430                                                            ql(kk,jj+1,ii+1) &
    431                                  ) / ( 3.0 * gg )
    432 
    433                       ql_int_u = ( (gg-aa) * ql(kk+1,jj,ii)   + (gg-bb) *   &
    434                                                          ql(kk+1,jj,ii+1)   &
    435                                  + (gg-cc) * ql(kk+1,jj+1,ii) + (gg-dd) *   &
    436                                                          ql(kk+1,jj+1,ii+1) &
    437                                  ) / ( 3.0 * gg )
    438 
    439                       ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz *&
    440                                           ( ql_int_u - ql_int_l )
    441 
    442 !
    443 !--                   Interpolate u velocity-component
    444                       ii = ( particles(n)%x + 0.5 * dx ) * ddx
    445                       jj =   particles(n)%y * ddy
    446                       kk = ( particles(n)%z + 0.5 * dz ) / dz ! only if eqist
    447 
    448                       IF ( ( particles(n)%z - zu(kk) ) > (0.5*dz) ) kk = kk+1
    449 
    450                       x  = particles(n)%x + ( 0.5 - ii ) * dx
    451                       y  = particles(n)%y - jj * dy
    452                       aa = x**2          + y**2
    453                       bb = ( dx - x )**2 + y**2
    454                       cc = x**2          + ( dy - y )**2
    455                       dd = ( dx - x )**2 + ( dy - y )**2
    456                       gg = aa + bb + cc + dd
    457 
    458                       u_int_l = ( (gg-aa) * u(kk,jj,ii)   + (gg-bb) *       &
    459                                                            u(kk,jj,ii+1)    &
    460                                 + (gg-cc) * u(kk,jj+1,ii) + (gg-dd) *       &
    461                                                            u(kk,jj+1,ii+1)  &
    462                                 ) / ( 3.0 * gg ) - u_gtrans
    463                       IF ( kk+1 == nzt+1 )  THEN
    464                          u_int = u_int_l
    465                       ELSE
    466                          u_int_u = ( (gg-aa) * u(kk+1,jj,ii)   + (gg-bb) *  &
    467                                                           u(kk+1,jj,ii+1)   &
    468                                    + (gg-cc) * u(kk+1,jj+1,ii) + (gg-dd) *  &
    469                                                           u(kk+1,jj+1,ii+1) &
    470                                    ) / ( 3.0 * gg ) - u_gtrans
    471                          u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz &
    472                                            * ( u_int_u - u_int_l )
    473                       ENDIF
    474 
    475 !
    476 !--                   Same procedure for interpolation of the v velocity-com-
    477 !--                   ponent (adopt index k from u velocity-component)
    478                       ii =   particles(n)%x * ddx
    479                       jj = ( particles(n)%y + 0.5 * dy ) * ddy
    480 
    481                       x  = particles(n)%x - ii * dx
    482                       y  = particles(n)%y + ( 0.5 - jj ) * dy
    483                       aa = x**2          + y**2
    484                       bb = ( dx - x )**2 + y**2
    485                       cc = x**2          + ( dy - y )**2
    486                       dd = ( dx - x )**2 + ( dy - y )**2
    487                       gg = aa + bb + cc + dd
    488 
    489                       v_int_l = ( ( gg-aa ) * v(kk,jj,ii)   + ( gg-bb ) *   &
    490                                                             v(kk,jj,ii+1)   &
    491                                 + ( gg-cc ) * v(kk,jj+1,ii) + ( gg-dd ) *   &
    492                                                             v(kk,jj+1,ii+1) &
    493                                 ) / ( 3.0 * gg ) - v_gtrans
    494                       IF ( kk+1 == nzt+1 )  THEN
    495                          v_int = v_int_l
    496                       ELSE
    497                          v_int_u = ( (gg-aa) * v(kk+1,jj,ii)   + (gg-bb) *  &
    498                                                             v(kk+1,jj,ii+1) &
    499                                    + (gg-cc) * v(kk+1,jj+1,ii) + (gg-dd) *  &
    500                                                           v(kk+1,jj+1,ii+1) &
    501                                    ) / ( 3.0 * gg ) - v_gtrans
    502                          v_int = v_int_l + ( particles(n)%z - zu(kk) ) / dz &
    503                                            * ( v_int_u - v_int_l )
    504                       ENDIF
    505 
    506 !
    507 !--                   Same procedure for interpolation of the w velocity-com-
    508 !--                   ponent (adopt index i from v velocity-component)
    509                       jj = particles(n)%y * ddy
    510                       kk = particles(n)%z / dz
    511 
    512                       x  = particles(n)%x - ii * dx
    513                       y  = particles(n)%y - jj * dy
    514                       aa = x**2          + y**2
    515                       bb = ( dx - x )**2 + y**2
    516                       cc = x**2          + ( dy - y )**2
    517                       dd = ( dx - x )**2 + ( dy - y )**2
    518                       gg = aa + bb + cc + dd
    519 
    520                       w_int_l = ( ( gg-aa ) * w(kk,jj,ii)   + ( gg-bb ) *   &
    521                                                               w(kk,jj,ii+1) &
    522                                 + ( gg-cc ) * w(kk,jj+1,ii) + ( gg-dd ) *   &
    523                                                             w(kk,jj+1,ii+1) &
    524                                 ) / ( 3.0 * gg )
    525                       IF ( kk+1 == nzt+1 )  THEN
    526                          w_int = w_int_l
    527                       ELSE
    528                          w_int_u = ( (gg-aa) * w(kk+1,jj,ii)   + (gg-bb) *  &
    529                                                             w(kk+1,jj,ii+1) &
    530                                    + (gg-cc) * w(kk+1,jj+1,ii) + (gg-dd) *  &
    531                                                           w(kk+1,jj+1,ii+1) &
    532                                    ) / ( 3.0 * gg )
    533                          w_int = w_int_l + ( particles(n)%z - zw(kk) ) / dz &
    534                                            * ( w_int_u - w_int_l )
    535                       ENDIF
    536 
    537 !
    538 !--                   Change in radius due to collision
    539                       delta_r = effective_coll_efficiency / 3.0_wp          &
    540                                 * pi * sl_r3 * ddx * ddy / dz               &
    541                                 * SQRT( ( u_int - particles(n)%speed_x )**2 &
    542                                       + ( v_int - particles(n)%speed_y )**2 &
    543                                       + ( w_int - particles(n)%speed_z )**2 &
    544                                       ) * dt_3d
    545 !
    546 !--                   Change in volume due to collision
    547                       delta_v = particles(n)%weight_factor                  &
    548                                 * ( ( particles(n)%radius + delta_r )**3    &
    549                                     - particles(n)%radius**3 )
    550 
    551 !
    552 !--                   Check if collected particles provide enough LWC for
    553 !--                   volume change of collector particle
    554                       IF ( delta_v >= sl_r3  .AND.  sl_r3 > 0.0 )  THEN
    555 
    556                          delta_r = ( ( sl_r3/particles(n)%weight_factor )           &
    557                                      + particles(n)%radius**3 )**( 1.0_wp/3.0_wp )  &
    558                                      - particles(n)%radius
    559 
    560                          DO  is = n-1, psi, -1
    561                             IF ( particles(is)%radius < &
    562                                  particles(n)%radius )  THEN
    563                                particles(is)%weight_factor = 0.0
    564                                particle_mask(is) = .FALSE.
    565                                deleted_particles = deleted_particles + 1
    566                             ENDIF
    567                          ENDDO
    568 
    569                       ELSE IF ( delta_v < sl_r3  .AND.  sl_r3 > 0.0 )  THEN
    570 
    571                          DO  is = n-1, psi, -1
    572                             IF ( particles(is)%radius < particles(n)%radius &
    573                                  .AND.  sl_r3 > 0.0 )  THEN
    574                                particles(is)%weight_factor =                &
    575                                             ( ( particles(is)%weight_factor &
    576                                             * ( particles(is)%radius**3 ) ) &
    577                                             - ( delta_v                     &
    578                                             * particles(is)%weight_factor   &
    579                                             * ( particles(is)%radius**3 )   &
    580                                             / sl_r3 ) )                     &
    581                                             / ( particles(is)%radius**3 )
    582 
    583                                IF ( particles(is)%weight_factor < 0.0 )  THEN
    584                                   WRITE( message_string, * ) 'negative ',   &
    585                                                   'weighting factor: ',     &
    586                                                   particles(is)%weight_factor
    587                                   CALL message( 'lpm_droplet_collision', &
    588                                                 'PA0039',                &
    589                                                 2, 2, -1, 6, 1 )
    590                                ENDIF
    591                             ENDIF
    592                          ENDDO
    593 
    594                       ENDIF
    595 
    596                       particles(n)%radius = particles(n)%radius + delta_r
    597                       ql_vp(k,j,i) = ql_vp(k,j,i) +               &
    598                                      particles(n)%weight_factor * &
    599                                      ( particles(n)%radius**3 )
    600                    ENDDO
    601 
    602                 ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor  &
    603                                               * particles(psi)%radius**3
    604 
    605                 ENDIF  ! collision kernel
    606 
    607              ELSE IF ( prt_count(k,j,i) == 1 )  THEN
    608 
    609                 psi = prt_start_index(k,j,i)
    610 
    611 !
    612 !--             Calculate change of weighting factor due to self collision
    613                 IF ( ( hall_kernel .OR. wang_kernel )  .AND.  &
    614                      use_kernel_tables )  THEN
    615 
    616                    IF ( wang_kernel )  THEN
    617                       eclass = INT( diss(k,j,i) * 1.0E4_wp / 1000.0_wp * &
    618                                     dissipation_classes ) + 1
    619                       epsilon = diss(k,j,i)
    620                    ELSE
    621                       epsilon = 0.0
    622                    ENDIF
    623                    IF ( hall_kernel .OR. epsilon * 1.0E4_wp < 0.001 )  THEN
    624                       eclass = 0   ! Hall kernel is used
    625                    ELSE
    626                       eclass = MIN( dissipation_classes, eclass )
    627                    ENDIF
    628 
    629                    ddV = ddx * ddy / dz
    630                    rclass_l = particles(psi)%class
    631                    sum3 = 1 - dt_3d * ddV *                                 &
    632                               ( ckernel(rclass_l,rclass_l,eclass) *         &
    633                                 ( particles(psi)%weight_factor-1 ) * 0.5 )
    634 
    635                    particles(psi)%radius = ( particles(psi)%radius**3 / &
    636                                              sum3 )**0.33333333333333_wp
    637                    particles(psi)%weight_factor = particles(psi)%weight_factor &
    638                                                   * sum3
    639 
    640                 ELSE IF ( ( hall_kernel .OR. wang_kernel )  .AND.  &
    641                            .NOT. use_kernel_tables )  THEN
    642 !
    643 !--                Collision efficiencies are calculated for every new
    644 !--                grid box. First, allocate memory for kernel table.
    645 !--                Third dimension is 1, because table is re-calculated for
    646 !--                every new dissipation value.
    647                    ALLOCATE( ckernel(psi:psi, psi:psi, 1:1) )
    648 !
    649 !--                Now calculate collision efficiencies for this box
    650                    CALL recalculate_kernel( i, j, k )
    651 
    652                    ddV = ddx * ddy / dz
    653                    sum3 = 1 - dt_3d * ddV * ( ckernel(psi,psi,1) *  &
    654                                 ( particles(psi)%weight_factor - 1 ) * 0.5 )
    655 
    656                    particles(psi)%radius = ( particles(psi)%radius**3 / &
    657                                              sum3 )**0.33333333333333_wp
    658                    particles(psi)%weight_factor = particles(psi)%weight_factor &
    659                                                   * sum3
    660 
    661                    DEALLOCATE( ckernel )
    662                 ENDIF
    663 
    664                ql_vp(k,j,i) =  particles(psi)%weight_factor *              &
    665                                 particles(psi)%radius**3
     597                ENDDO
     598
    666599             ENDIF
    667600
    668 !
    669 !--          Check if condensation of LWC was conserved during collision
    670 !--          process
    671              IF ( ql_v(k,j,i) /= 0.0 )  THEN
    672                 IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001  .OR.             &
    673                      ql_vp(k,j,i) / ql_v(k,j,i) <= 0.9999 )  THEN
    674                    WRITE( message_string, * ) 'LWC is not conserved during',&
    675                                               ' collision! ',               &
    676                                               'LWC after condensation: ',   &
    677                                               ql_v(k,j,i),                  &
    678                                               ' LWC after collision: ',     &
    679                                               ql_vp(k,j,i)
    680                    CALL message( 'lpm_droplet_collision', 'PA0040',         &
    681                                  2, 2, -1, 6, 1 )
    682                 ENDIF
    683              ENDIF
    684 
     601             particles(n)%radius = particles(n)%radius + delta_r
     602             ql_vp(k,j,i) = ql_vp(k,j,i) +               &
     603                            particles(n)%weight_factor * &
     604                            ( particles(n)%radius**3 )
    685605          ENDDO
    686        ENDDO
    687     ENDDO
     606
     607          ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor  &
     608                                        * particles(psi)%radius**3
     609
     610       ENDIF  ! collision kernel
     611
     612    ELSE IF ( prt_count(k,j,i) == 1 )  THEN
     613
     614       psi = 1
     615
     616!
     617!--    Calculate change of weighting factor due to self collision
     618       IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
     619            use_kernel_tables )  THEN
     620
     621          IF ( wang_kernel )  THEN
     622             eclass = INT( diss(k,j,i) * 1.0E4_wp / 1000.0_wp * &
     623                           dissipation_classes ) + 1
     624             epsilon = diss(k,j,i)
     625          ELSE
     626             epsilon = 0.0_wp
     627          ENDIF
     628          IF ( hall_kernel  .OR.  epsilon * 1.0E4_wp < 0.001_wp )  THEN
     629             eclass = 0   ! Hall kernel is used
     630          ELSE
     631             eclass = MIN( dissipation_classes, eclass )
     632          ENDIF
     633
     634          ddV = ddx * ddy / dz
     635          rclass_l = particles(psi)%class
     636          sum3 = 1 - dt_3d * ddV *                                 &
     637                     ( ckernel(rclass_l,rclass_l,eclass) *         &
     638                       ( particles(psi)%weight_factor-1 ) * 0.5_wp )
     639
     640          particles(psi)%radius = ( particles(psi)%radius**3 / &
     641                                    sum3 )**0.33333333333333_wp
     642          particles(psi)%weight_factor = particles(psi)%weight_factor &
     643                                         * sum3
     644
     645       ELSE IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
     646                  .NOT. use_kernel_tables )  THEN
     647!
     648!--       Collision efficiencies are calculated for every new
     649!--       grid box. First, allocate memory for kernel table.
     650!--       Third dimension is 1, because table is re-calculated for
     651!--       every new dissipation value.
     652          ALLOCATE( ckernel(psi:psi, psi:psi, 1:1) )
     653!
     654!--       Now calculate collision efficiencies for this box
     655          CALL recalculate_kernel( i, j, k )
     656
     657          ddV = ddx * ddy / dz
     658          sum3 = 1 - dt_3d * ddV * ( ckernel(psi,psi,1) *  &
     659                       ( particles(psi)%weight_factor - 1 ) * 0.5_wp )
     660
     661          particles(psi)%radius = ( particles(psi)%radius**3 / &
     662                                    sum3 )**0.33333333333333_wp
     663          particles(psi)%weight_factor = particles(psi)%weight_factor &
     664                                         * sum3
     665
     666          DEALLOCATE( ckernel )
     667       ENDIF
     668
     669      ql_vp(k,j,i) =  particles(psi)%weight_factor *              &
     670                       particles(psi)%radius**3
     671    ENDIF
     672
     673!
     674!-- Check if condensation of LWC was conserved during collision process
     675    IF ( ql_v(k,j,i) /= 0.0_wp )  THEN
     676       IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001_wp  .OR.             &
     677            ql_vp(k,j,i) / ql_v(k,j,i) <= 0.9999_wp )  THEN
     678          WRITE( message_string, * ) 'LWC is not conserved during',&
     679                                     ' collision! ',               &
     680                                     'LWC after condensation: ',   &
     681                                     ql_v(k,j,i),                  &
     682                                     ' LWC after collision: ',     &
     683                                     ql_vp(k,j,i)
     684          CALL message( 'lpm_droplet_collision', 'PA0040',         &
     685                        2, 2, -1, 6, 1 )
     686       ENDIF
     687    ENDIF
    688688
    689689    CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'stop' )
    690690
    691 
    692691 END SUBROUTINE lpm_droplet_collision
Note: See TracChangeset for help on using the changeset viewer.