Changeset 1071 for palm/trunk/SOURCE/lpm_droplet_collision.f90
- Timestamp:
- Nov 29, 2012 4:54:55 PM (11 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/lpm_droplet_collision.f90 (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lpm_droplet_collision.f90
r1037 r1071 20 20 ! Current revisions: 21 21 ! ------------------ 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 23 26 ! 24 27 ! Former revisions: … … 35 38 ! Description: 36 39 ! ------------ 37 ! Calculates chang in droplet radius by collision. Droplet collision is40 ! Calculates change in droplet radius by collision. Droplet collision is 38 41 ! calculated for each grid box seperately. Collision is parameterized by 39 42 ! using collision kernels. Three different kernels are available: … … 74 77 mean_r, ql_int, ql_int_l, ql_int_u, u_int, u_int_l, u_int_u, & 75 78 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 77 80 78 81 TYPE(particle_type) :: tmp_particle 82 REAL, DIMENSION(:), ALLOCATABLE :: rad, weight 79 83 80 84 … … 138 142 ENDIF 139 143 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 144 161 rclass_l = particles(n)%class 145 162 ! 146 !-- Calculate growth of collector particle radius using all 147 !-- droplets smaller than current droplet 163 !-- Mass added due to collisions with smaller droplets 148 164 DO is = psi, n-1 149 165 150 166 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 ) 161 170 162 171 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 225 206 226 207 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) 227 213 228 214 ELSEIF ( ( hall_kernel .OR. wang_kernel ) .AND. & … … 241 227 CALL recalculate_kernel( i, j, k ) 242 228 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 250 247 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 ) 262 251 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 324 283 325 284 ENDDO 326 285 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 ) 328 290 329 291 ELSEIF ( palm_kernel ) THEN … … 543 505 'weighting factor: ', & 544 506 particles(is)%weight_factor 545 CALL message( 'lpm_droplet_collision', '', & 507 CALL message( 'lpm_droplet_collision', & 508 'PA0039', & 546 509 2, 2, -1, 6, 1 ) 547 510 ENDIF … … 557 520 ENDDO 558 521 559 ENDIF ! collision kernel560 561 522 ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor & 562 523 * particles(psi)%radius**3 563 524 525 ENDIF ! collision kernel 564 526 565 527 ELSE IF ( prt_count(k,j,i) == 1 ) THEN 566 528 567 529 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 * & 569 585 particles(psi)%radius**3 570 586 ENDIF … … 582 598 ' LWC after collision: ', & 583 599 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 ) 585 602 ENDIF 586 603 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
