Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Ignore:
Timestamp:
Apr 7, 2016 7:49:42 AM (8 years ago)
Author:
hoffmann
Message:

changes in LPM and bulk cloud microphysics

File:
1 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/lpm_droplet_collision.f90

    r1818 r1822  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Integration of a new collision algortithm based on Shima et al. (2009) and
     22! Soelch and Kaercher (2010) called all_or_nothing. The previous implemented
     23! collision algorithm is called average_impact. Moreover, both algorithms are
     24! now positive definit due to their construction, i.e., no negative weighting
     25! factors should occur.
    2226!
    2327! Former revisions:
    … …  
    6367!> Calculates change in droplet radius by collision. Droplet collision is
    6468!> calculated for each grid box seperately. Collision is parameterized by
    65 !> using collision kernels. Three different kernels are available:
    66 !> PALM kernel: Kernel is approximated using a method from Rogers and
    67 !>              Yau (1989, A Short Course in Cloud Physics, Pergamon Press).
    68 !>              All droplets smaller than the treated one are represented by
    69 !>              one droplet with mean features. Collision efficiencies are taken
    70 !>              from the respective table in Rogers and Yau.
     69!> using collision kernels. Two different kernels are available:
    7170!> Hall kernel: Kernel from Hall (1980, J. Atmos. Sci., 2486-2507), which
    7271!>              considers collision due to pure gravitational effects.
    … …  
    8584
    8685    USE arrays_3d,                                                             &
    87         ONLY:  diss, ql, ql_v, ql_vp, u, v, w, zu, zw
     86        ONLY:  diss, ql_v, ql_vp
    8887
    8988    USE cloud_parameters,                                                      &
    90         ONLY:  effective_coll_efficiency
     89        ONLY:  rho_l
    9190
    9291    USE constants,                                                             &
    … …  
    9493
    9594    USE control_parameters,                                                    &
    96         ONLY:  dt_3d, message_string, u_gtrans, v_gtrans, dz
     95        ONLY:  dt_3d, message_string, dz
    9796
    9897    USE cpulog,                                                                &
    … …  
    10099
    101100    USE grid_variables,                                                        &
    102         ONLY:  ddx, dx, ddy, dy
    103 
    104     USE indices,                                                               &
    105         ONLY:  nxl, nxr, nyn, nys, nzb, nzt
     101        ONLY:  dx, dy
    106102
    107103    USE kinds
    108104
    109105    USE lpm_collision_kernels_mod,                                             &
    110         ONLY:  ckernel, collision_efficiency_rogers, recalculate_kernel
     106        ONLY:  ckernel, recalculate_kernel
    111107
    112108    USE particle_attributes,                                                   &
    113         ONLY:  deleted_particles, dissipation_classes, hall_kernel,            &
    114                palm_kernel, particles, particle_type,                          &
    115                prt_count, use_kernel_tables, wang_kernel
     109        ONLY:  all_or_nothing, average_impact, dissipation_classes,            &
     110               hall_kernel, iran_part, number_of_particles, particles,         &
     111               particle_type, prt_count, use_kernel_tables, wang_kernel
     112
     113    USE random_function_mod,                                                   &
     114        ONLY:  random_function
    116115
    117116    USE pegrid
    … …  
    121120    INTEGER(iwp) ::  eclass   !<
    122121    INTEGER(iwp) ::  i        !<
    123     INTEGER(iwp) ::  ii       !<
    124     INTEGER(iwp) ::  inc      !<
    125     INTEGER(iwp) ::  is       !<
    126122    INTEGER(iwp) ::  j        !<
    127     INTEGER(iwp) ::  jj       !<
    128     INTEGER(iwp) ::  js       !<
    129123    INTEGER(iwp) ::  k        !<
    130     INTEGER(iwp) ::  kk       !<
    131124    INTEGER(iwp) ::  n        !<
    132     INTEGER(iwp) ::  pse      !<
    133     INTEGER(iwp) ::  psi      !<
     125    INTEGER(iwp) ::  m        !<
    134126    INTEGER(iwp) ::  rclass_l !<
    135127    INTEGER(iwp) ::  rclass_s !<
    136128
    137     INTEGER(iwp), DIMENSION(prt_count(k,j,i)) ::  rclass_v !<
    138 
    139     LOGICAL, SAVE ::  first_flag = .TRUE. !<
    140 
    141     TYPE(particle_type) :: tmp_particle   !<
    142 
    143     REAL(wp) ::  aa       !<
    144     REAL(wp) ::  auxn     !< temporary variables
    145     REAL(wp) ::  auxs     !< temporary variables
    146     REAL(wp) ::  bb       !<
    147     REAL(wp) ::  cc       !<
    148     REAL(wp) ::  dd       !<
    149     REAL(wp) ::  ddV      !<
    150     REAL(wp) ::  delta_r  !<
    151     REAL(wp) ::  delta_v  !<
    152     REAL(wp) ::  epsilon  !<
    153     REAL(wp) ::  gg       !<
    154     REAL(wp) ::  mean_r   !<
    155     REAL(wp) ::  ql_int   !<
    156     REAL(wp) ::  ql_int_l !<
    157     REAL(wp) ::  ql_int_u !<
    158     REAL(wp) ::  r3       !<
    159     REAL(wp) ::  sl_r3    !<
    160     REAL(wp) ::  sl_r4    !<
    161     REAL(wp) ::  sum1     !<
    162     REAL(wp) ::  sum2     !<
    163     REAL(wp) ::  sum3     !<
    164     REAL(wp) ::  u_int    !<
    165     REAL(wp) ::  u_int_l  !<
    166     REAL(wp) ::  u_int_u  !<
    167     REAL(wp) ::  v_int    !<
    168     REAL(wp) ::  v_int_l  !<
    169     REAL(wp) ::  v_int_u  !<
    170     REAL(wp) ::  w_int    !<
    171     REAL(wp) ::  w_int_l  !<
    172     REAL(wp) ::  w_int_u  !<
    173     REAL(wp) ::  x        !<
    174     REAL(wp) ::  y        !<
    175 
    176     REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad    !<
    177     REAL(wp), DIMENSION(:), ALLOCATABLE ::  weight !<
    178 
    179     REAL, DIMENSION(prt_count(k,j,i))    :: ck
    180     REAL, DIMENSION(prt_count(k,j,i))    :: r3v
    181     REAL, DIMENSION(prt_count(k,j,i))    :: sum1v
    182     REAL, DIMENSION(prt_count(k,j,i))    :: sum2v
     129    REAL(wp) ::  collection_probability  !< probability for collection
     130    REAL(wp) ::  ddV                     !< inverse grid box volume
     131    REAL(wp) ::  epsilon                 !< dissipation rate
     132    REAL(wp) ::  factor_volume_to_mass   !< 4.0 / 3.0 * pi * rho_l
     133    REAL(wp) ::  xm                      !< mean mass of droplet m
     134    REAL(wp) ::  xn                      !< mean mass of droplet n
     135
     136    REAL(wp), DIMENSION(:), ALLOCATABLE ::  weight  !< weighting factor
     137    REAL(wp), DIMENSION(:), ALLOCATABLE ::  mass    !< total mass of super droplet
    183138
    184139    CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'start' )
    185140
    186 !
    187 !-- Collision requires at least two particles in the box
    188     IF ( prt_count(k,j,i) > 1 )  THEN
    189 !
    190 !--    First, sort particles within the gridbox by their size,
    191 !--    using Shell's method (see Numerical Recipes)
    192 !--    NOTE: In case of using particle tails, the re-sorting of
    193 !--    ----  tails would have to be included here!
    194        IF ( .NOT. ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
    195        use_kernel_tables ) )  THEN
    196           psi = 0
    197           inc = 1
    198           DO WHILE ( inc <= prt_count(k,j,i) )
    199              inc = 3 * inc + 1
    200           ENDDO
    201 
    202           DO WHILE ( inc > 1 )
    203              inc = inc / 3
    204              DO  is = inc+1, prt_count(k,j,i)
    205                 tmp_particle = particles(psi+is)
    206                 js = is
    207                 DO WHILE ( particles(psi+js-inc)%radius >             &
    208                 tmp_particle%radius )
    209                    particles(psi+js) = particles(psi+js-inc)
    210                    js = js - inc
    211                    IF ( js <= inc )  EXIT
    212                 ENDDO
    213                 particles(psi+js) = tmp_particle
    214              ENDDO
    215           ENDDO
    216        ENDIF
    217 
    218        psi = 1
    219        pse = prt_count(k,j,i)
     141    number_of_particles   = prt_count(k,j,i)
     142    factor_volume_to_mass = 4.0_wp / 3.0_wp * pi * rho_l
     143    ddV                   = 1 / ( dx * dy * dz )
     144!
     145!-- Collision requires at least one super droplet inside the box
     146    IF ( number_of_particles > 0 )  THEN
    220147
    221148!
    222149!--    Now apply the different kernels
    223        IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
    224             use_kernel_tables )  THEN
    225 !
    226 !--       Fast method with pre-calculated efficiencies for
     150       IF ( use_kernel_tables )  THEN
     151!
     152!--       Fast method with pre-calculated collection kernels for
    227153!--       discrete radius- and dissipation-classes.
    228154!--
    … …  
    244170!--       Droplet collision are calculated using collision-coalescence
    245171!--       formulation proposed by Wang (see PALM documentation)
    246 !--       Since new radii after collision are defined by radii of all
    247 !--       droplets before collision, temporary fields for new radii and
    248 !--       weighting factors are needed
    249           ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i)))
    250 
    251           rad    = 0.0_wp
    252           weight = 0.0_wp
    253 
    254           sum1v(1:prt_count(k,j,i)) = 0.0_wp
    255           sum2v(1:prt_count(k,j,i)) = 0.0_wp
    256 
    257           DO  n = 1, prt_count(k,j,i)
    258 
    259              rclass_l = particles(n)%class
    260 !
    261 !--          Mass added due to collisions with smaller droplets
    262              DO  is = n+1, prt_count(k,j,i)
    263                 rclass_s = particles(is)%class
    264                 auxs = ckernel(rclass_l,rclass_s,eclass) * particles(is)%weight_factor
    265                 auxn = ckernel(rclass_s,rclass_l,eclass) * particles(n)%weight_factor
    266                 IF ( particles(is)%radius <  particles(n)%radius )  THEN
    267                    sum1v(n) =  sum1v(n)  + particles(is)%radius**3 * auxs
    268                    sum2v(is) = sum2v(is) + auxn
    269                 ELSE
    270                    sum2v(n)  = sum2v(n)  + auxs
    271                    sum1v(is) = sum1v(is) + particles(n)%radius**3 * auxn
    272                 ENDIF
     172!--       Temporary fields for total mass of super-droplet and weighting factors
     173!--       are allocated.
     174          ALLOCATE(mass(1:number_of_particles), weight(1:number_of_particles))
     175
     176          mass(1:number_of_particles)   = particles(1:number_of_particles)%weight_factor * &
     177                                          particles(1:number_of_particles)%radius**3     * &
     178                                          factor_volume_to_mass
     179          weight(1:number_of_particles) = particles(1:number_of_particles)%weight_factor
     180
     181          IF ( average_impact )  THEN  ! select collision algorithm
     182
     183             DO  n = 1, number_of_particles
     184
     185                rclass_l = particles(n)%class
     186                xn       = mass(n) / weight(n)
     187
     188                DO  m = n, number_of_particles
     189
     190                   rclass_s = particles(m)%class
     191                   xm       = mass(m) / weight(m)
     192
     193                   IF ( xm .LT. xn )  THEN
     194                     
     195!
     196!--                   Particle n collects smaller particle m
     197                      collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
     198                                               weight(n) * ddV * dt_3d
     199
     200                      mass(n)   = mass(n)   + mass(m)   * collection_probability
     201                      weight(m) = weight(m) - weight(m) * collection_probability
     202                      mass(m)   = mass(m)   - mass(m)   * collection_probability
     203                   ELSEIF ( xm .GT. xn )  THEN
     204!
     205!--                   Particle m collects smaller particle n
     206                      collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
     207                                               weight(m) * ddV * dt_3d
     208
     209                      mass(m)   = mass(m)   + mass(n)   * collection_probability
     210                      weight(n) = weight(n) - weight(n) * collection_probability
     211                      mass(n)   = mass(n)   - mass(n)   * collection_probability
     212                   ELSE
     213!
     214!--                   Same-size collections. If n = m, weight is reduced by the
     215!--                   number of possible same-size collections; the total mass
     216!--                   is not changed during same-size collection.
     217!--                   Same-size collections of different
     218!--                   particles ( n /= m ) are treated as same-size collections
     219!--                   of ONE partilce with weight = weight(n) + weight(m) and
     220!--                   mass = mass(n) + mass(m).
     221!--                   Accordingly, each particle loses the same number of
     222!--                   droplets to the other particle, but this has no effect on
     223!--                   total mass mass, since the exchanged droplets have the
     224!--                   same radius.
     225
     226!--                   Note: For m = n this equation is an approximation only
     227!--                   valid for weight >> 1 (which is usually the case). The
     228!--                   approximation is weight(n)-1 = weight(n).
     229                      weight(n) = weight(n) - 0.5_wp * weight(n) *                &
     230                                              ckernel(rclass_l,rclass_s,eclass) * &
     231                                              weight(m) * ddV * dt_3d
     232                      IF ( n .NE. m )  THEN
     233                         weight(m) = weight(m) - 0.5_wp * weight(m) *                &
     234                                                 ckernel(rclass_l,rclass_s,eclass) * &
     235                                                 weight(n) * ddV * dt_3d
     236                      ENDIF
     237                   ENDIF
     238
     239                ENDDO
     240
     241                ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
     242
    273243             ENDDO
    274           ENDDO
    275           rclass_v = particles(1:prt_count(k,j,i))%class
    276           DO  n = 1, prt_count(k,j,i)
    277             ck(n)       = ckernel(rclass_v(n),rclass_v(n),eclass)
    278           ENDDO
    279           r3v      = particles(1:prt_count(k,j,i))%radius**3
    280           DO  n = 1, prt_count(k,j,i)
    281              sum3 = 0.0_wp
    282              ddV = ddx * ddy / dz
    283 !
    284 !--          Change of the current weighting factor
    285              sum3 = 1 - dt_3d * ddV *                                 &
    286                         ck(n) *                                       &
    287                         ( particles(n)%weight_factor - 1 ) * 0.5_wp - &
    288                     dt_3d * ddV * sum2v(n)
    289              weight(n) = particles(n)%weight_factor * sum3
    290 !
    291 !--          Change of the current droplet radius
    292              rad(n) = ( (r3v(n) + dt_3d * ddV * (sum1v(n) - sum2v(n) * r3v(n)) )/&
    293                               sum3 )**0.33333333333333_wp
    294 
    295              ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n) &
    296                                          * rad(n)**3
    297 
    298           ENDDO
     244
     245          ELSEIF ( all_or_nothing )  THEN  ! select collision algorithm
     246
     247             DO  n = 1, number_of_particles
     248
     249                rclass_l = particles(n)%class
     250                xn       = mass(n) / weight(n) ! mean mass of droplet n
     251
     252                DO  m = n, number_of_particles
     253
     254                   rclass_s = particles(m)%class
     255                   xm = mass(m) / weight(m) ! mean mass of droplet m
     256
     257                   IF ( weight(n) .LT. weight(m) )  THEN
     258!
     259!--                   Particle n collects weight(n) droplets of particle m 
     260                      collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
     261                                               weight(m) * ddV * dt_3d
     262
     263                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
     264                         mass(n)   = mass(n)   + weight(n) * xm
     265                         weight(m) = weight(m) - weight(n)
     266                         mass(m)   = mass(m)   - weight(n) * xm
     267                      ENDIF
     268
     269                   ELSEIF ( weight(m) .LT. weight(n) )  THEN
     270!
     271!--                   Particle m collects weight(m) droplets of particle n 
     272                      collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
     273                                               weight(n) * ddV * dt_3d
     274
     275                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
     276                         mass(m)   = mass(m)   + weight(m) * xn
     277                         weight(n) = weight(n) - weight(m)
     278                         mass(n)   = mass(n)   - weight(m) * xn
     279                      ENDIF
     280                   ELSE
     281!
     282!--                   Collisions of particles of the same weighting factor.
     283!--                   Particle n collects 1/2 weight(n) droplets of particle m,
     284!--                   particle m collects 1/2 weight(m) droplets of particle n.
     285!--                   The total mass mass changes accordingly.
     286!--                   If n = m, the first half of the droplets coalesces with the
     287!--                   second half of the droplets; mass is unchanged because
     288!--                   xm = xn for n = m.
     289
     290!--                   Note: For m = n this equation is an approximation only
     291!--                   valid for weight >> 1 (which is usually the case). The
     292!--                   approximation is weight(n)-1 = weight(n).
     293                      collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
     294                                               weight(n) * ddV * dt_3d
     295
     296                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
     297                         mass(n)   = mass(n)   + 0.5_wp * weight(n) * ( xm - xn )
     298                         mass(m)   = mass(m)   + 0.5_wp * weight(m) * ( xn - xm )
     299                         weight(n) = weight(n) - 0.5_wp * weight(m)
     300                         weight(m) = weight(n)
     301                      ENDIF
     302                   ENDIF
     303
     304                ENDDO
     305
     306                ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
     307
     308             ENDDO
     309
     310          ENDIF
     311
     312
     313
     314
    299315          IF ( ANY(weight < 0.0_wp) )  THEN
    300316                WRITE( message_string, * ) 'negative weighting'
    … …  
    303319          ENDIF
    304320
    305           particles(psi:pse)%radius = rad(1:prt_count(k,j,i))
    306           particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i))
    307 
    308           DEALLOCATE(rad, weight)
    309 
    310        ELSEIF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
    311                 .NOT. use_kernel_tables )  THEN
    312 !
    313 !--       Collision efficiencies are calculated for every new
     321          particles(1:number_of_particles)%radius = ( mass(1:number_of_particles) /   &
     322                                                      ( weight(1:number_of_particles) &
     323                                                        * factor_volume_to_mass       &
     324                                                      )                               &
     325                                                    )**0.33333333333333_wp
     326
     327          particles(1:number_of_particles)%weight_factor = weight(1:number_of_particles)
     328
     329          DEALLOCATE(weight, mass)
     330
     331       ELSEIF ( .NOT. use_kernel_tables )  THEN
     332!
     333!--       Collection kernels are calculated for every new
    314334!--       grid box. First, allocate memory for kernel table.
    315335!--       Third dimension is 1, because table is re-calculated for
    316336!--       every new dissipation value.
    317           ALLOCATE( ckernel(1:prt_count(k,j,i),1:prt_count(k,j,i),1:1) )
    318 !
    319 !--       Now calculate collision efficiencies for this box
     337          ALLOCATE( ckernel(1:number_of_particles,1:number_of_particles,1:1) )
     338!
     339!--       Now calculate collection kernel for this box. Note that
     340!--       the kernel is based on the previous time step
    320341          CALL recalculate_kernel( i, j, k )
    321 
    322342!
    323343!--       Droplet collision are calculated using collision-coalescence
    324344!--       formulation proposed by Wang (see PALM documentation)
    325 !--       Since new radii after collision are defined by radii of all
    326 !--       droplets before collision, temporary fields for new radii and
    327 !--       weighting factors are needed
    328           ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i)))
    329 
    330           rad = 0.0_wp
    331           weight = 0.0_wp
    332 
    333           DO  n = psi, pse, 1
    334 
    335              sum1 = 0.0_wp
    336              sum2 = 0.0_wp
    337              sum3 = 0.0_wp
    338 !
    339 !--          Mass added due to collisions with smaller droplets
    340              DO  is = psi, n-1
    341                 sum1 = sum1 + ( particles(is)%radius**3 *            &
    342                                 ckernel(n,is,1) *  &
    343                                 particles(is)%weight_factor )
     345!--       Temporary fields for total mass of super-droplet and weighting factors
     346!--       are allocated.
     347          ALLOCATE(mass(1:number_of_particles), weight(1:number_of_particles))
     348
     349          mass(1:number_of_particles) = particles(1:number_of_particles)%weight_factor * &
     350                                        particles(1:number_of_particles)%radius**3     * &
     351                                        factor_volume_to_mass
     352
     353          weight(1:number_of_particles) = particles(1:number_of_particles)%weight_factor
     354
     355          IF ( average_impact )  THEN  ! select collision algorithm
     356
     357             DO  n = 1, number_of_particles
     358
     359                xn = mass(n) / weight(n) ! mean mass of droplet n
     360
     361                DO  m = n, number_of_particles
     362
     363                   xm = mass(m) / weight(m) !mean mass of droplet m
     364
     365                   IF ( xm .LT. xn )  THEN
     366!
     367!--                   Particle n collects smaller particle m
     368                      collection_probability = ckernel(n,m,1) * weight(n) *    &
     369                                               ddV * dt_3d
     370
     371                      mass(n)   = mass(n)   + mass(m)   * collection_probability
     372                      weight(m) = weight(m) - weight(m) * collection_probability
     373                      mass(m)   = mass(m)   - mass(m)   * collection_probability
     374                   ELSEIF ( xm .GT. xn )  THEN
     375!
     376!--                   Particle m collects smaller particle n
     377                      collection_probability = ckernel(n,m,1) * weight(m) *    &
     378                                               ddV * dt_3d
     379
     380                      mass(m)   = mass(m)   + mass(n)   * collection_probability
     381                      weight(n) = weight(n) - weight(n) * collection_probability
     382                      mass(n)   = mass(n)   - mass(n)   * collection_probability
     383                   ELSE
     384!
     385!--                   Same-size collections. If n = m, weight is reduced by the
     386!--                   number of possible same-size collections; the total mass
     387!--                   mass is not changed during same-size collection.
     388!--                   Same-size collections of different
     389!--                   particles ( n /= m ) are treated as same-size collections
     390!--                   of ONE partilce with weight = weight(n) + weight(m) and
     391!--                   mass = mass(n) + mass(m).
     392!--                   Accordingly, each particle loses the same number of
     393!--                   droplets to the other particle, but this has no effect on
     394!--                   total mass mass, since the exchanged droplets have the
     395!--                   same radius.
     396!--
     397!--                   Note: For m = n this equation is an approximation only
     398!--                   valid for weight >> 1 (which is usually the case). The
     399!--                   approximation is weight(n)-1 = weight(n).
     400                      weight(n) = weight(n) - 0.5_wp * weight(n) *             &
     401                                              ckernel(n,m,1) * weight(m) *     &
     402                                              ddV * dt_3d
     403                      IF ( n .NE. m )  THEN
     404                         weight(m) = weight(m) - 0.5_wp * weight(m) *          &
     405                                                 ckernel(n,m,1) * weight(n) *  &
     406                                                 ddV * dt_3d
     407                      ENDIF
     408                   ENDIF
     409
     410
     411                ENDDO
     412
     413                ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
     414
    344415             ENDDO
    345 !
    346 !--          Rate of collisions with larger droplets
    347              DO  is = n+1, pse
    348                 sum2 = sum2 + ( ckernel(n,is,1) *  &
    349                                 particles(is)%weight_factor )
     416
     417          ELSEIF ( all_or_nothing )  THEN  ! select collision algorithm
     418
     419             DO  n = 1, number_of_particles
     420
     421                xn = mass(n) / weight(n) ! mean mass of droplet n
     422
     423                DO  m = n, number_of_particles
     424
     425                   xm = mass(m) / weight(m) !mean mass of droplet m
     426
     427                   IF ( weight(n) .LT. weight(m) )  THEN
     428!
     429!--                   Particle n collects smaller particle m
     430                      collection_probability = ckernel(n,m,1) * weight(m) *    &
     431                                               ddV * dt_3d
     432
     433                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
     434                         mass(n) = mass(n) + weight(n) * xm 
     435                         weight(m)    = weight(m)    - weight(n)
     436                         mass(m) = mass(m) - weight(n) * xm
     437                      ENDIF
     438
     439                   ELSEIF ( weight(m) .LT. weight(n) )  THEN
     440!
     441!--                   Particle m collects smaller particle n
     442                      collection_probability = ckernel(n,m,1) * weight(n) *    &
     443                                               ddV * dt_3d
     444
     445                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
     446                         mass(m) = mass(m) + weight(m) * xn
     447                         weight(n)    = weight(n)    - weight(m)
     448                         mass(n) = mass(n) - weight(m) * xn
     449                      ENDIF
     450                   ELSE
     451!
     452!--                   Collisions of particles of the same weighting factor.
     453!--                   Particle n collects 1/2 weight(n) droplets of particle m,
     454!--                   particle m collects 1/2 weight(m) droplets of particle n.
     455!--                   The total mass mass changes accordingly.
     456!--                   If n = m, the first half of the droplets coalesces with the
     457!--                   second half of the droplets; mass is unchanged because
     458!--                   xm = xn for n = m.
     459!--
     460!--                   Note: For m = n this equation is an approximation only
     461!--                   valid for weight >> 1 (which is usually the case). The
     462!--                   approximation is weight(n)-1 = weight(n).
     463                      collection_probability = ckernel(n,m,1) * weight(n) *    &
     464                                               ddV * dt_3d
     465
     466                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
     467                         mass(n)   = mass(n)   + 0.5_wp * weight(n) * ( xm - xn )
     468                         mass(m)   = mass(m)   + 0.5_wp * weight(m) * ( xn - xm )
     469                         weight(n) = weight(n) - 0.5_wp * weight(m)
     470                         weight(m) = weight(n)
     471                      ENDIF
     472                   ENDIF
     473
     474
     475                ENDDO
     476
     477                ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
     478
    350479             ENDDO
    351480
    352              r3 = particles(n)%radius**3
    353              ddV = ddx * ddy / dz
    354              is = 1
    355 !
    356 !--                   Change of the current weighting factor
    357              sum3 = 1 - dt_3d * ddV *                                 &
    358                         ckernel(n,n,1) *           &
    359                         ( particles(n)%weight_factor - 1 ) * 0.5_wp - &
    360                     dt_3d * ddV * sum2
    361              weight(n-is+1) = particles(n)%weight_factor * sum3
    362 !
    363 !--                   Change of the current droplet radius
    364              rad(n-is+1) = ( (r3 + dt_3d * ddV * (sum1 - sum2 * r3) )/&
    365                               sum3 )**0.33333333333333_wp
    366 
    367              IF ( weight(n-is+1) < 0.0_wp )  THEN
    368                 WRITE( message_string, * ) 'negative weighting',      &
    369                                            'factor: ', weight(n-is+1)
    370                 CALL message( 'lpm_droplet_collision', 'PA0037',      &
     481          ENDIF
     482
     483          IF ( ANY(weight < 0.0_wp) )  THEN
     484                WRITE( message_string, * ) 'negative weighting'
     485                CALL message( 'lpm_droplet_collision', 'PA0028',      &
    371486                               2, 2, -1, 6, 1 )
    372              ENDIF
    373 
    374              ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n-is+1) &
    375                                          * rad(n-is+1)**3
    376 
    377           ENDDO
    378 
    379           particles(psi:pse)%radius = rad(1:prt_count(k,j,i))
    380           particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i))
    381 
    382           DEALLOCATE( rad, weight, ckernel )
    383 
    384        ELSEIF ( palm_kernel )  THEN
    385 !
    386 !--       PALM collision kernel
    387 !
    388 !--       Calculate the mean radius of all those particles which
    389 !--       are of smaller size than the current particle and
    390 !--       use this radius for calculating the collision efficiency
    391           DO  n = psi+prt_count(k,j,i)-1, psi+1, -1
    392 
    393              sl_r3 = 0.0_wp
    394              sl_r4 = 0.0_wp
    395 
    396              DO is = n-1, psi, -1
    397                 IF ( particles(is)%radius < particles(n)%radius )  &
    398                 THEN
    399                    sl_r3 = sl_r3 + particles(is)%weight_factor     &
    400                                    * particles(is)%radius**3
    401                    sl_r4 = sl_r4 + particles(is)%weight_factor     &
    402                                    * particles(is)%radius**4
    403                 ENDIF
    404              ENDDO
    405 
    406              IF ( ( sl_r3 ) > 0.0_wp )  THEN
    407                 mean_r = ( sl_r4 ) / ( sl_r3 )
    408 
    409                 CALL collision_efficiency_rogers( mean_r,             &
    410                                            particles(n)%radius,       &
    411                                            effective_coll_efficiency )
    412 
    413              ELSE
    414                 effective_coll_efficiency = 0.0_wp
    415              ENDIF
    416 
    417              IF ( effective_coll_efficiency > 1.0_wp  .OR.            &
    418                   effective_coll_efficiency < 0.0_wp )                &
    419              THEN
    420                 WRITE( message_string, * )  'collision_efficien' , &
    421                           'cy out of range:' ,effective_coll_efficiency
    422                 CALL message( 'lpm_droplet_collision', 'PA0145', 2, &
    423                               2, -1, 6, 1 )
    424              ENDIF
    425 
    426 !
    427 !--          Interpolation of liquid water content
    428              ii = particles(n)%x * ddx
    429              jj = particles(n)%y * ddy
    430              kk = ( particles(n)%z + 0.5_wp * dz ) / dz
    431 
    432              x  = particles(n)%x - ii * dx
    433              y  = particles(n)%y - jj * dy
    434              aa = x**2          + y**2
    435              bb = ( dx - x )**2 + y**2
    436              cc = x**2          + ( dy - y )**2
    437              dd = ( dx - x )**2 + ( dy - y )**2
    438              gg = aa + bb + cc + dd
    439 
    440              ql_int_l = ( (gg-aa) * ql(kk,jj,ii)   + (gg-bb) *     &
    441                                                   ql(kk,jj,ii+1)   &
    442                         + (gg-cc) * ql(kk,jj+1,ii) + ( gg-dd ) *   &
    443                                                   ql(kk,jj+1,ii+1) &
    444                         ) / ( 3.0_wp * gg )
    445 
    446              ql_int_u = ( (gg-aa) * ql(kk+1,jj,ii)   + (gg-bb) *   &
    447                                                 ql(kk+1,jj,ii+1)   &
    448                         + (gg-cc) * ql(kk+1,jj+1,ii) + (gg-dd) *   &
    449                                                 ql(kk+1,jj+1,ii+1) &
    450                         ) / ( 3.0_wp * gg )
    451 
    452              ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz *&
    453                                  ( ql_int_u - ql_int_l )
    454 
    455 !
    456 !--          Interpolate u velocity-component
    457              ii = ( particles(n)%x + 0.5_wp * dx ) * ddx
    458              jj =   particles(n)%y * ddy
    459              kk = ( particles(n)%z + 0.5_wp * dz ) / dz ! only if equidistant
    460 
    461              IF ( ( particles(n)%z - zu(kk) ) > ( 0.5_wp * dz ) ) kk = kk+1
    462 
    463              x  = particles(n)%x + ( 0.5_wp - ii ) * dx
    464              y  = particles(n)%y - jj * dy
    465              aa = x**2          + y**2
    466              bb = ( dx - x )**2 + y**2
    467              cc = x**2          + ( dy - y )**2
    468              dd = ( dx - x )**2 + ( dy - y )**2
    469              gg = aa + bb + cc + dd
    470 
    471              u_int_l = ( (gg-aa) * u(kk,jj,ii)   + (gg-bb) *       &
    472                                                   u(kk,jj,ii+1)    &
    473                        + (gg-cc) * u(kk,jj+1,ii) + (gg-dd) *       &
    474                                                   u(kk,jj+1,ii+1)  &
    475                        ) / ( 3.0_wp * gg ) - u_gtrans
    476              IF ( kk+1 == nzt+1 )  THEN
    477                 u_int = u_int_l
    478              ELSE
    479                 u_int_u = ( (gg-aa) * u(kk+1,jj,ii)   + (gg-bb) *  &
    480                                                  u(kk+1,jj,ii+1)   &
    481                           + (gg-cc) * u(kk+1,jj+1,ii) + (gg-dd) *  &
    482                                                  u(kk+1,jj+1,ii+1) &
    483                           ) / ( 3.0_wp * gg ) - u_gtrans
    484                 u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz &
    485                                   * ( u_int_u - u_int_l )
    486              ENDIF
    487 
    488 !
    489 !--          Same procedure for interpolation of the v velocity-component
    490 !--          (adopt index k from u velocity-component)
    491              ii =   particles(n)%x * ddx
    492              jj = ( particles(n)%y + 0.5_wp * dy ) * ddy
    493 
    494              x  = particles(n)%x - ii * dx
    495              y  = particles(n)%y + ( 0.5_wp - jj ) * dy
    496              aa = x**2          + y**2
    497              bb = ( dx - x )**2 + y**2
    498              cc = x**2          + ( dy - y )**2
    499              dd = ( dx - x )**2 + ( dy - y )**2
    500              gg = aa + bb + cc + dd
    501 
    502              v_int_l = ( ( gg-aa ) * v(kk,jj,ii)   + ( gg-bb ) *   &
    503                                                    v(kk,jj,ii+1)   &
    504                        + ( gg-cc ) * v(kk,jj+1,ii) + ( gg-dd ) *   &
    505                                                    v(kk,jj+1,ii+1) &
    506                        ) / ( 3.0_wp * gg ) - v_gtrans
    507              IF ( kk+1 == nzt+1 )  THEN
    508                 v_int = v_int_l
    509              ELSE
    510                 v_int_u = ( (gg-aa) * v(kk+1,jj,ii)   + (gg-bb) *  &
    511                                                    v(kk+1,jj,ii+1) &
    512                           + (gg-cc) * v(kk+1,jj+1,ii) + (gg-dd) *  &
    513                                                  v(kk+1,jj+1,ii+1) &
    514                           ) / ( 3.0_wp * gg ) - v_gtrans
    515                 v_int = v_int_l + ( particles(n)%z - zu(kk) ) / dz &
    516                                   * ( v_int_u - v_int_l )
    517              ENDIF
    518 
    519 !
    520 !--          Same procedure for interpolation of the w velocity-component
    521 !--          (adopt index i from v velocity-component)
    522              jj = particles(n)%y * ddy
    523              kk = particles(n)%z / dz
    524 
    525              x  = particles(n)%x - ii * dx
    526              y  = particles(n)%y - jj * dy
    527              aa = x**2          + y**2
    528              bb = ( dx - x )**2 + y**2
    529              cc = x**2          + ( dy - y )**2
    530              dd = ( dx - x )**2 + ( dy - y )**2
    531              gg = aa + bb + cc + dd
    532 
    533              w_int_l = ( ( gg-aa ) * w(kk,jj,ii)   + ( gg-bb ) *   &
    534                                                      w(kk,jj,ii+1) &
    535                        + ( gg-cc ) * w(kk,jj+1,ii) + ( gg-dd ) *   &
    536                                                    w(kk,jj+1,ii+1) &
    537                        ) / ( 3.0_wp * gg )
    538              IF ( kk+1 == nzt+1 )  THEN
    539                 w_int = w_int_l
    540              ELSE
    541                 w_int_u = ( (gg-aa) * w(kk+1,jj,ii)   + (gg-bb) *  &
    542                                                    w(kk+1,jj,ii+1) &
    543                           + (gg-cc) * w(kk+1,jj+1,ii) + (gg-dd) *  &
    544                                                  w(kk+1,jj+1,ii+1) &
    545                           ) / ( 3.0_wp * gg )
    546                 w_int = w_int_l + ( particles(n)%z - zw(kk) ) / dz &
    547                                   * ( w_int_u - w_int_l )
    548              ENDIF
    549 
    550 !
    551 !--          Change in radius due to collision
    552              delta_r = effective_coll_efficiency / 3.0_wp          &
    553                        * pi * sl_r3 * ddx * ddy / dz               &
    554                        * SQRT( ( u_int - particles(n)%speed_x )**2 &
    555                              + ( v_int - particles(n)%speed_y )**2 &
    556                              + ( w_int - particles(n)%speed_z )**2 &
    557                              ) * dt_3d
    558 !
    559 !--          Change in volume due to collision
    560              delta_v = particles(n)%weight_factor                  &
    561                        * ( ( particles(n)%radius + delta_r )**3    &
    562                            - particles(n)%radius**3 )
    563 
    564 !
    565 !--          Check if collected particles provide enough LWC for
    566 !--          volume change of collector particle
    567              IF ( delta_v >= sl_r3  .AND.  sl_r3 > 0.0_wp )  THEN
    568 
    569                 delta_r = ( ( sl_r3/particles(n)%weight_factor )               &
    570                             + particles(n)%radius**3 )**( 1.0_wp / 3.0_wp )    &
    571                             - particles(n)%radius
    572 
    573                 DO  is = n-1, psi, -1
    574                    IF ( particles(is)%radius < particles(n)%radius )  THEN
    575                       particles(is)%weight_factor = 0.0_wp
    576                       particles(is)%particle_mask = .FALSE.
    577                       deleted_particles = deleted_particles + 1
    578                    ENDIF
    579                 ENDDO
    580 
    581              ELSE IF ( delta_v < sl_r3  .AND.  sl_r3 > 0.0_wp )  THEN
    582 
    583                 DO  is = n-1, psi, -1
    584                    IF ( particles(is)%radius < particles(n)%radius &
    585                         .AND.  sl_r3 > 0.0_wp )  THEN
    586                       particles(is)%weight_factor =                &
    587                                    ( ( particles(is)%weight_factor &
    588                                    * ( particles(is)%radius**3 ) ) &
    589                                    - ( delta_v                     &
    590                                    * particles(is)%weight_factor   &
    591                                    * ( particles(is)%radius**3 )   &
    592                                    / sl_r3 ) )                     &
    593                                    / ( particles(is)%radius**3 )
    594 
    595                       IF ( particles(is)%weight_factor < 0.0_wp )  THEN
    596                          WRITE( message_string, * ) 'negative ',   &
    597                                          'weighting factor: ',     &
    598                                          particles(is)%weight_factor
    599                          CALL message( 'lpm_droplet_collision', &
    600                                        'PA0039',                &
    601                                        2, 2, -1, 6, 1 )
    602                       ENDIF
    603                    ENDIF
    604                 ENDDO
    605 
    606              ENDIF
    607 
    608              particles(n)%radius = particles(n)%radius + delta_r
    609              ql_vp(k,j,i) = ql_vp(k,j,i) +               &
    610                             particles(n)%weight_factor * &
    611                             ( particles(n)%radius**3 )
    612           ENDDO
    613 
    614           ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor  &
    615                                         * particles(psi)%radius**3
    616 
    617        ENDIF  ! collision kernel
    618 
    619     ELSE IF ( prt_count(k,j,i) == 1 )  THEN
    620 
    621        psi = 1
    622 
    623 !
    624 !--    Calculate change of weighting factor due to self collision
    625        IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
    626             use_kernel_tables )  THEN
    627 
    628           IF ( wang_kernel )  THEN
    629              eclass = INT( diss(k,j,i) * 1.0E4_wp / 1000.0_wp * &
    630                            dissipation_classes ) + 1
    631              epsilon = diss(k,j,i)
    632           ELSE
    633              epsilon = 0.0_wp
    634487          ENDIF
    635           IF ( hall_kernel  .OR.  epsilon * 1.0E4_wp < 0.001_wp )  THEN
    636              eclass = 0   ! Hall kernel is used
    637           ELSE
    638              eclass = MIN( dissipation_classes, eclass )
    639           ENDIF
    640 
    641           ddV = ddx * ddy / dz
    642           rclass_l = particles(psi)%class
    643           sum3 = 1 - dt_3d * ddV *                                 &
    644                      ( ckernel(rclass_l,rclass_l,eclass) *         &
    645                        ( particles(psi)%weight_factor-1 ) * 0.5_wp )
    646 
    647           particles(psi)%radius = ( particles(psi)%radius**3 / &
    648                                     sum3 )**0.33333333333333_wp
    649           particles(psi)%weight_factor = particles(psi)%weight_factor &
    650                                          * sum3
    651 
    652        ELSE IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
    653                   .NOT. use_kernel_tables )  THEN
    654 !
    655 !--       Collision efficiencies are calculated for every new
    656 !--       grid box. First, allocate memory for kernel table.
    657 !--       Third dimension is 1, because table is re-calculated for
    658 !--       every new dissipation value.
    659           ALLOCATE( ckernel(psi:psi, psi:psi, 1:1) )
    660 !
    661 !--       Now calculate collision efficiencies for this box
    662           CALL recalculate_kernel( i, j, k )
    663 
    664           ddV = ddx * ddy / dz
    665           sum3 = 1 - dt_3d * ddV * ( ckernel(psi,psi,1) *  &
    666                        ( particles(psi)%weight_factor - 1 ) * 0.5_wp )
    667 
    668           particles(psi)%radius = ( particles(psi)%radius**3 / &
    669                                     sum3 )**0.33333333333333_wp
    670           particles(psi)%weight_factor = particles(psi)%weight_factor &
    671                                          * sum3
    672 
    673           DEALLOCATE( ckernel )
    674        ENDIF
    675 
    676       ql_vp(k,j,i) =  particles(psi)%weight_factor *              &
    677                        particles(psi)%radius**3
     488
     489          particles(1:number_of_particles)%radius = ( mass(1:number_of_particles) /   &
     490                                                      ( weight(1:number_of_particles) &
     491                                                        * factor_volume_to_mass       &
     492                                                      )                               &
     493                                                    )**0.33333333333333_wp
     494
     495          particles(1:number_of_particles)%weight_factor = weight(1:number_of_particles)
     496
     497          DEALLOCATE( weight, mass, ckernel )
     498
     499       ENDIF
     500 
    678501    ENDIF
    679 
    680 !
    681 !-- Check if condensation of LWC was conserved during collision process
     502   
     503
     504!
     505!-- Check if LWC is conserved during collision process
    682506    IF ( ql_v(k,j,i) /= 0.0_wp )  THEN
    683        IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001_wp  .OR.             &
     507       IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001_wp  .OR.                      &
    684508            ql_vp(k,j,i) / ql_v(k,j,i) <= 0.9999_wp )  THEN
    685           WRITE( message_string, * ) 'LWC is not conserved during',&
    686                                      ' collision! ',               &
    687                                      'LWC after condensation: ',   &
    688                                      ql_v(k,j,i),                  &
    689                                      ' LWC after collision: ',     &
    690                                      ql_vp(k,j,i)
    691           CALL message( 'lpm_droplet_collision', 'PA0040',         &
    692                         2, 2, -1, 6, 1 )
     509          WRITE( message_string, * ) ' LWC is not conserved during',           &
     510                                     ' collision! ',                           &
     511                                     ' LWC after condensation: ', ql_v(k,j,i), &
     512                                     ' LWC after collision: ', ql_vp(k,j,i)
     513          CALL message( 'lpm_droplet_collision', 'PA0040', 2, 2, -1, 6, 1 )
    693514       ENDIF
    694515    ENDIF
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |