Ignore:
Timestamp:
Nov 29, 2012 4:54:55 PM (9 years ago)
Author:
franke
Message:

ventilation effect for cloud droplets included, calculation of cloud droplet collisions now uses formulation by Wang, bugfixes in Rosenbrock method and in calculation of collision efficiencies

File:
1 edited

Legend:

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

    r1037 r1071  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Calculation of Hall and Wang kernel now uses collision-coalescence formulation
     23! proposed by Wang instead of the continuous collection equation (for more
     24! information about new method see PALM documentation)
     25! Bugfix: message identifiers added
    2326!
    2427! Former revisions:
     
    3538! Description:
    3639! ------------
    37 ! Calculates chang in droplet radius by collision. Droplet collision is
     40! Calculates change in droplet radius by collision. Droplet collision is
    3841! calculated for each grid box seperately. Collision is parameterized by
    3942! using collision kernels. Three different kernels are available:
     
    7477             mean_r, ql_int, ql_int_l, ql_int_u, u_int, u_int_l, u_int_u,     &
    7578             v_int, v_int_l, v_int_u, w_int, w_int_l, w_int_u, sl_r3, sl_r4,  &
    76              x, y
     79             x, y, sum1, sum2, sum3, r3, ddV
    7780
    7881    TYPE(particle_type) ::  tmp_particle
     82    REAL, DIMENSION(:), ALLOCATABLE :: rad, weight
    7983
    8084
     
    138142                   ENDIF
    139143
    140                    DO  n = pse, psi+1, -1
    141 
    142                       integral = 0.0
    143                       lw_max   = 0.0
     144 !
     145!--                Droplet collision are calculated using collision-coalescence
     146!--                formulation proposed by Wang (see PALM documentation)
     147!--                Since new radii after collision are defined by radii of all
     148!--                droplets before collision, temporary fields for new radii and
     149!--                weighting factors are needed
     150                   ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i)))
     151
     152                   rad = 0.0
     153                   weight = 0.0
     154
     155                   DO  n = psi, pse, 1
     156
     157                      sum1 = 0.0
     158                      sum2 = 0.0
     159                      sum3 = 0.0
     160
    144161                      rclass_l = particles(n)%class
    145162!
    146 !--                   Calculate growth of collector particle radius using all
    147 !--                   droplets smaller than current droplet
     163!--                   Mass added due to collisions with smaller droplets
    148164                      DO  is = psi, n-1
    149165
    150166                         rclass_s = particles(is)%class
    151                          integral = integral +                             &
    152                                     ( particles(is)%radius**3 *            &
    153                                       ckernel(rclass_l,rclass_s,eclass) *  &
    154                                       particles(is)%weight_factor )
    155 !
    156 !--                      Calculate volume of liquid water of the collected
    157 !--                      droplets which is the maximum liquid water available
    158 !--                      for droplet growth   
    159                          lw_max = lw_max + ( particles(is)%radius**3 * &
    160                                              particles(is)%weight_factor )
     167                         sum1 = sum1 + ( particles(is)%radius**3 *            &
     168                                         ckernel(rclass_l,rclass_s,eclass) *  &
     169                                         particles(is)%weight_factor )
    161170
    162171                      ENDDO
    163 
    164 !
    165 !--                   Change in radius of collector droplet due to collision
    166                       delta_r = 1.0 / ( 3.0 * particles(n)%radius**2 ) * &
    167                                 integral * dt_3d * ddx * ddy / dz
    168 
    169 !
    170 !--                   Change in volume of collector droplet due to collision
    171                       delta_v = particles(n)%weight_factor                  &
    172                                 * ( ( particles(n)%radius + delta_r )**3    &
    173                                     - particles(n)%radius**3 )
    174 
    175                       IF ( lw_max < delta_v  .AND.  delta_v > 0.0 )  THEN
    176 !-- replace by message call
    177                          PRINT*, 'Particle has grown to much because the',  &
    178                                  ' change of volume of particle is larger', &
    179                                  ' than liquid water available!'
    180 
    181                       ELSEIF ( lw_max == delta_v  .AND.  delta_v > 0.0 ) THEN
    182 !-- can this case really happen??
    183                          DO  is = psi, n-1
    184 
    185                             particles(is)%weight_factor = 0.0
    186                             particle_mask(is)  = .FALSE.
    187                             deleted_particles = deleted_particles + 1
    188 
    189                          ENDDO
    190 
    191                       ELSEIF ( lw_max > delta_v  .AND.  delta_v > 0.0 )  THEN
    192 !
    193 !--                      Calculate new weighting factor of collected droplets
    194                          DO  is = psi, n-1
    195 
    196                             rclass_s = particles(is)%class
    197                             particles(is)%weight_factor = &
    198                                particles(is)%weight_factor - &
    199                                ( ( ckernel(rclass_l,rclass_s,eclass) * particles(is)%weight_factor &
    200                                  / integral ) * delta_v )
    201 
    202                             IF ( particles(is)%weight_factor < 0.0 )  THEN
    203                                   WRITE( message_string, * ) 'negative ',   &
    204                                                      'weighting factor: ',  &
    205                                                   particles(is)%weight_factor
    206                                   CALL message( 'lpm_droplet_collision', '', &
    207                                                 2, 2, -1, 6, 1 )
    208 
    209                             ELSEIF ( particles(is)%weight_factor == 0.0 )  &
    210                             THEN
    211 
    212                                 particles(is)%weight_factor = 0.0
    213                                 particle_mask(is)  = .FALSE.
    214                                 deleted_particles = deleted_particles + 1
    215 
    216                             ENDIF
    217 
    218                          ENDDO
    219 
    220                       ENDIF
    221 
    222                       particles(n)%radius = particles(n)%radius + delta_r
    223                       ql_vp(k,j,i) = ql_vp(k,j,i) + particles(n)%weight_factor &
    224                                                     * particles(n)%radius**3
     172!
     173!--                   Rate of collisions with larger droplets
     174                      DO  is = n+1, pse
     175
     176                         rclass_s = particles(is)%class
     177                         sum2 = sum2 + ( ckernel(rclass_l,rclass_s,eclass) *  &
     178                                         particles(is)%weight_factor )
     179
     180                      ENDDO
     181
     182                      r3 = particles(n)%radius**3
     183                      ddV = ddx * ddy / dz
     184                      is = prt_start_index(k,j,i)
     185!
     186!--                   Change of the current weighting factor
     187                      sum3 = 1 - dt_3d * ddV *                                 &
     188                                 ckernel(rclass_l,rclass_l,eclass) *           &
     189                                 ( particles(n)%weight_factor - 1 ) * 0.5 -    &
     190                             dt_3d * ddV * sum2
     191                      weight(n-is+1) = particles(n)%weight_factor * sum3
     192!
     193!--                   Change of the current droplet radius
     194                      rad(n-is+1) = ( (r3 + dt_3d * ddV * (sum1 - sum2 * r3) )/&
     195                                       sum3 )**0.33333333333333
     196
     197                      IF ( weight(n-is+1) < 0.0 )  THEN
     198                         WRITE( message_string, * ) 'negative weighting',      &
     199                                                    'factor: ', weight(n-is+1)
     200                         CALL message( 'lpm_droplet_collision', 'PA0028',      &
     201                                        2, 2, -1, 6, 1 )
     202                      ENDIF
     203
     204                      ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n-is+1) &
     205                                                  * rad(n-is+1)**3
    225206
    226207                   ENDDO
     208
     209                   particles(psi:pse)%radius = rad(1:prt_count(k,j,i))
     210                   particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i))
     211
     212                   DEALLOCATE(rad, weight)
    227213
    228214                ELSEIF ( ( hall_kernel .OR. wang_kernel )  .AND.  &
     
    241227                   CALL recalculate_kernel( i, j, k )
    242228
    243                    DO  n = pse, psi+1, -1
    244 
    245                       integral = 0.0
    246                       lw_max   = 0.0
    247 !
    248 !--                   Calculate growth of collector particle radius using all
    249 !--                   droplets smaller than current droplet
     229!
     230!--                Droplet collision are calculated using collision-coalescence
     231!--                formulation proposed by Wang (see PALM documentation)
     232!--                Since new radii after collision are defined by radii of all
     233!--                droplets before collision, temporary fields for new radii and
     234!--                weighting factors are needed
     235                   ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i)))
     236
     237                   rad = 0.0
     238                   weight = 0.0
     239
     240                   DO  n = psi, pse, 1
     241
     242                      sum1 = 0.0
     243                      sum2 = 0.0
     244                      sum3 = 0.0
     245!
     246!--                   Mass added due to collisions with smaller droplets
    250247                      DO  is = psi, n-1
    251 
    252                          integral = integral + particles(is)%radius**3 * &
    253                                                ckernel(n,is,1) *         &
    254                                                particles(is)%weight_factor
    255 !
    256 !--                      Calculate volume of liquid water of the collected
    257 !--                      droplets which is the maximum liquid water available
    258 !--                      for droplet growth   
    259                          lw_max = lw_max + ( particles(is)%radius**3 * &
    260                                              particles(is)%weight_factor )
    261 
     248                         sum1 = sum1 + ( particles(is)%radius**3 *            &
     249                                         ckernel(n,is,1) *  &
     250                                         particles(is)%weight_factor )
    262251                      ENDDO
    263 
    264 !
    265 !--                   Change in radius of collector droplet due to collision
    266                       delta_r = 1.0 / ( 3.0 * particles(n)%radius**2 ) * &
    267                                 integral * dt_3d * ddx * ddy / dz
    268 
    269 !
    270 !--                   Change in volume of collector droplet due to collision
    271                       delta_v = particles(n)%weight_factor                  &
    272                                 * ( ( particles(n)%radius + delta_r )**3    &
    273                                     - particles(n)%radius**3 )
    274 
    275                       IF ( lw_max < delta_v  .AND.  delta_v > 0.0 )  THEN
    276 !-- replace by message call
    277                          PRINT*, 'Particle has grown to much because the',  &
    278                                  ' change of volume of particle is larger', &
    279                                  ' than liquid water available!'
    280 
    281                       ELSEIF ( lw_max == delta_v  .AND.  delta_v > 0.0 ) THEN
    282 !-- can this case really happen??
    283                          DO  is = psi, n-1
    284 
    285                             particles(is)%weight_factor = 0.0
    286                             particle_mask(is)  = .FALSE.
    287                             deleted_particles = deleted_particles + 1
    288 
    289                          ENDDO
    290 
    291                       ELSEIF ( lw_max > delta_v  .AND.  delta_v > 0.0 )  THEN
    292 !
    293 !--                      Calculate new weighting factor of collected droplets
    294                          DO  is = psi, n-1
    295 
    296                             particles(is)%weight_factor =                   &
    297                                    particles(is)%weight_factor -            &
    298                                    ( ckernel(n,is,1) / integral * delta_v * &
    299                                      particles(is)%weight_factor )
    300 
    301                             IF ( particles(is)%weight_factor < 0.0 )  THEN
    302                                   WRITE( message_string, * ) 'negative ',   &
    303                                                      'weighting factor: ',  &
    304                                                   particles(is)%weight_factor
    305                                   CALL message( 'lpm_droplet_collision', '', &
    306                                                 2, 2, -1, 6, 1 )
    307 
    308                             ELSEIF ( particles(is)%weight_factor == 0.0 )  &
    309                             THEN
    310 
    311                                 particles(is)%weight_factor = 0.0
    312                                 particle_mask(is)  = .FALSE.
    313                                 deleted_particles = deleted_particles + 1
    314 
    315                             ENDIF
    316 
    317                          ENDDO
    318 
    319                       ENDIF
    320 
    321                       particles(n)%radius = particles(n)%radius + delta_r
    322                       ql_vp(k,j,i) = ql_vp(k,j,i) + particles(n)%weight_factor &
    323                                                     * particles(n)%radius**3
     252!
     253!--                   Rate of collisions with larger droplets
     254                      DO  is = n+1, pse
     255                         sum2 = sum2 + ( ckernel(n,is,1) *  &
     256                                         particles(is)%weight_factor )
     257                      ENDDO
     258
     259                      r3 = particles(n)%radius**3
     260                      ddV = ddx * ddy / dz
     261                      is = prt_start_index(k,j,i)
     262!
     263!--                   Change of the current weighting factor
     264                      sum3 = 1 - dt_3d * ddV *                                 &
     265                                 ckernel(n,n,1) *           &
     266                                 ( particles(n)%weight_factor - 1 ) * 0.5 -    &
     267                             dt_3d * ddV * sum2
     268                      weight(n-is+1) = particles(n)%weight_factor * sum3
     269!
     270!--                   Change of the current droplet radius
     271                      rad(n-is+1) = ( (r3 + dt_3d * ddV * (sum1 - sum2 * r3) )/&
     272                                       sum3 )**0.33333333333333
     273
     274                      IF ( weight(n-is+1) < 0.0 )  THEN
     275                         WRITE( message_string, * ) 'negative weighting',      &
     276                                                    'factor: ', weight(n-is+1)
     277                         CALL message( 'lpm_droplet_collision', 'PA0037',      &
     278                                        2, 2, -1, 6, 1 )
     279                      ENDIF
     280
     281                      ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n-is+1) &
     282                                                  * rad(n-is+1)**3
    324283
    325284                   ENDDO
    326285
    327                    DEALLOCATE( ckernel )
     286                   particles(psi:pse)%radius = rad(1:prt_count(k,j,i))
     287                   particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i))
     288
     289                   DEALLOCATE( rad, weight, ckernel )
    328290
    329291                ELSEIF ( palm_kernel )  THEN
     
    543505                                                  'weighting factor: ',     &
    544506                                                  particles(is)%weight_factor
    545                                   CALL message( 'lpm_droplet_collision', '', &
     507                                  CALL message( 'lpm_droplet_collision', &
     508                                                'PA0039',                &
    546509                                                2, 2, -1, 6, 1 )
    547510                               ENDIF
     
    557520                   ENDDO
    558521
    559                 ENDIF  ! collision kernel
    560 
    561522                ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor  &
    562523                                              * particles(psi)%radius**3
    563524
     525                ENDIF  ! collision kernel
    564526
    565527             ELSE IF ( prt_count(k,j,i) == 1 )  THEN
    566528
    567529                psi = prt_start_index(k,j,i)
    568                 ql_vp(k,j,i) =  particles(psi)%weight_factor *              &
     530
     531!
     532!--             Calculate change of weighting factor due to self collision
     533                IF ( ( hall_kernel .OR. wang_kernel )  .AND.  &
     534                     use_kernel_tables )  THEN
     535
     536                   IF ( wang_kernel )  THEN
     537                      eclass = INT( diss(k,j,i) * 1.0E4 / 1000.0 * &
     538                                    dissipation_classes ) + 1
     539                      epsilon = diss(k,j,i)
     540                   ELSE
     541                      epsilon = 0.0
     542                   ENDIF
     543                   IF ( hall_kernel .OR. epsilon * 1.0E4 < 0.001 )  THEN
     544                      eclass = 0   ! Hall kernel is used
     545                   ELSE
     546                      eclass = MIN( dissipation_classes, eclass )
     547                   ENDIF
     548
     549                   ddV = ddx * ddy / dz
     550                   rclass_l = particles(psi)%class
     551                   sum3 = 1 - dt_3d * ddV *                                 &
     552                              ( ckernel(rclass_l,rclass_l,eclass) *         &
     553                                ( particles(psi)%weight_factor-1 ) * 0.5 )
     554
     555                   particles(psi)%radius = ( particles(psi)%radius**3 / &
     556                                             sum3 )**0.33333333333333
     557                   particles(psi)%weight_factor = particles(psi)%weight_factor &
     558                                                  * sum3
     559
     560                ELSE IF ( ( hall_kernel .OR. wang_kernel )  .AND.  &
     561                           .NOT. use_kernel_tables )  THEN
     562!
     563!--                Collision efficiencies are calculated for every new
     564!--                grid box. First, allocate memory for kernel table.
     565!--                Third dimension is 1, because table is re-calculated for
     566!--                every new dissipation value.
     567                   ALLOCATE( ckernel(psi:psi, psi:psi, 1:1) )
     568!
     569!--                Now calculate collision efficiencies for this box
     570                   CALL recalculate_kernel( i, j, k )
     571
     572                   ddV = ddx * ddy / dz
     573                   sum3 = 1 - dt_3d * ddV * ( ckernel(psi,psi,1) *  &
     574                                ( particles(psi)%weight_factor - 1 ) * 0.5 )
     575
     576                   particles(psi)%radius = ( particles(psi)%radius**3 / &
     577                                             sum3 )**0.33333333333333
     578                   particles(psi)%weight_factor = particles(psi)%weight_factor &
     579                                                  * sum3
     580
     581                   DEALLOCATE( ckernel )
     582                ENDIF
     583
     584               ql_vp(k,j,i) =  particles(psi)%weight_factor *              &
    569585                                particles(psi)%radius**3
    570586             ENDIF
     
    582598                                              ' LWC after collision: ',     &
    583599                                              ql_vp(k,j,i)
    584                    CALL message( 'lpm_droplet_collision', '', 2, 2, -1, 6, 1 )
     600                   CALL message( 'lpm_droplet_collision', 'PA0040',         &
     601                                 2, 2, -1, 6, 1 )
    585602                ENDIF
    586603             ENDIF
Note: See TracChangeset for help on using the changeset viewer.