Changeset 825 for palm/trunk
- Timestamp:
- Feb 19, 2012 3:03:44 AM (13 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 10 edited
- 1 moved
-
Makefile (modified) (8 diffs)
-
advec_particles.f90 (modified) (11 diffs)
-
check_parameters.f90 (modified) (2 diffs)
-
data_output_ptseries.f90 (modified) (14 diffs)
-
diffusion_e.f90 (modified) (5 diffs)
-
header.f90 (modified) (2 diffs)
-
init_3d_model.f90 (modified) (2 diffs)
-
lpm_collision_kernels.f90 (moved) (moved from palm/trunk/SOURCE/wang_kernel.f90) (3 diffs)
-
modules.f90 (modified) (5 diffs)
-
package_parin.f90 (modified) (3 diffs)
-
time_integration.f90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/Makefile
r791 r825 4 4 # Current revisions: 5 5 # ------------------ 6 # 6 # wang_kernel renamed lpm_collision_kernels 7 7 # 8 8 # Former revisions: … … 83 83 init_masks.f90 init_ocean.f90 init_particles.f90 init_pegrid.f90 \ 84 84 init_pt_anomaly.f90 init_rankine.f90 init_slope.f90 \ 85 interaction_droplets_ptq.f90 local_flush.f90 \86 local_getenv.f90 local_stop.f90 local_system.f90 local_tremain.f90 \87 local_tremain_ini.f90 message.f90 modules.f90 netcdf.f90 \85 interaction_droplets_ptq.f90 local_flush.f90 local_getenv.f90 \ 86 local_stop.f90 local_system.f90 local_tremain.f90 local_tremain_ini.f90 \ 87 lpm_collision_kernels.f90 message.f90 modules.f90 netcdf.f90 \ 88 88 package_parin.f90 palm.f90 parin.f90 particle_boundary_conds.f90 \ 89 89 plant_canopy_model.f90 poisfft.f90 \ … … 108 108 user_particle_attributes.f90 user_read_restart_data.f90 \ 109 109 user_spectra.f90 user_statistics.f90 wall_fluxes.f90 \ 110 w ang_kernel.f90 write_3d_binary.f90 write_compressed.f90 \110 write_3d_binary.f90 write_compressed.f90 \ 111 111 write_var_list.f90 112 112 … … 124 124 data_output_3d.o diffusion_e.o diffusion_s.o diffusion_u.o \ 125 125 diffusion_v.o diffusion_w.o diffusivities.o disturb_field.o \ 126 disturb_heatflux.o eqn_state_seawater.o exchange_horiz.o exchange_horiz_2d.o fft_xy.o \ 127 flow_statistics.o global_min_max.o header.o impact_of_latent_heat.o \ 128 inflow_turbulence.o init_1d_model.o init_3d_model.o init_advec.o init_cloud_physics.o \ 129 init_coupling.o init_dvrp.o init_grid.o init_masks.o init_ocean.o init_particles.o init_pegrid.o \ 130 init_pt_anomaly.o init_rankine.o init_slope.o \ 126 disturb_heatflux.o eqn_state_seawater.o exchange_horiz.o \ 127 exchange_horiz_2d.o fft_xy.o flow_statistics.o global_min_max.o \ 128 header.o impact_of_latent_heat.o inflow_turbulence.o init_1d_model.o \ 129 init_3d_model.o init_advec.o init_cloud_physics.o init_coupling.o \ 130 init_dvrp.o init_grid.o init_masks.o init_ocean.o init_particles.o \ 131 init_pegrid.o init_pt_anomaly.o init_rankine.o init_slope.o \ 131 132 interaction_droplets_ptq.o local_flush.o local_getenv.o local_stop.o \ 132 local_system.o local_tremain.o local_tremain_ini.o message.o modules.o \ 133 local_system.o local_tremain.o local_tremain_ini.o \ 134 lpm_collision_kernels.f90 message.o modules.o \ 133 135 netcdf.o package_parin.o palm.o parin.o particle_boundary_conds.o \ 134 136 plant_canopy_model.o poisfft.o \ … … 151 153 user_module.o user_parin.o user_particle_attributes.o \ 152 154 user_read_restart_data.o user_spectra.o user_statistics.o \ 153 wall_fluxes.o w ang_kernel.o write_3d_binary.o write_compressed.o \155 wall_fluxes.o write_3d_binary.o write_compressed.o \ 154 156 write_var_list.o 155 157 … … 181 183 182 184 183 advec_particles.o: modules.o random_function.o wang_kernel.o185 advec_particles.o: modules.o lpm_collision_kernels.o random_function.o 184 186 advec_s_bc.o: modules.o 185 187 advec_s_pw.o: modules.o … … 257 259 local_tremain.o: modules.o 258 260 local_tremain_ini.o: modules.o 261 lpm_collision_kernels.o: modules.o user_module.o 259 262 message.o: modules.o 260 263 modules.o: modules.f90 … … 328 331 user_statistics.o: modules.o user_module.o 329 332 wall_fluxes.o: modules.o 330 wang_kernel.o: modules.o user_module.o331 333 write_3d_binary.o: modules.o random_function.o 332 334 write_compressed.o: modules.o -
palm/trunk/SOURCE/advec_particles.f90
r824 r825 4 4 ! Current revisions: 5 5 ! ------------------ 6 ! droplet growth by condensation may include curvature and solution effects. 7 ! initialisation of temporary particle array for resorting removed. 8 ! particle attributes speed_x|y|z_sgs renamed rvar1|2|3. 6 ! droplet growth by condensation may include curvature and solution effects, 7 ! initialisation of temporary particle array for resorting removed, 8 ! particle attributes speed_x|y|z_sgs renamed rvar1|2|3, 9 ! module wang_kernel_mod renamed lpm_collision_kernels_mod, 10 ! wang_collision_kernel renamed wang_kernel 9 11 ! 10 12 ! Former revisions: … … 130 132 USE indices 131 133 USE interfaces 134 USE lpm_collision_kernels_mod 132 135 USE netcdf_control 133 136 USE particle_attributes … … 135 138 USE random_function_mod 136 139 USE statistics 137 USE wang_kernel_mod138 140 139 141 IMPLICIT NONE … … 200 202 CALL cpu_log( log_point(25), 'advec_particles', 'start' ) 201 203 202 ! IF ( number_of_particles /= number_of_tails ) THEN 203 ! WRITE (9,*) '--- advec_particles: #1' 204 ! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails 205 ! CALL local_flush( 9 ) 206 ! ENDIF 204 207 205 ! 208 206 !-- Write particle data on file for later analysis. … … 270 268 CALL cpu_log( log_point_s(41), 'advec_part_exp', 'stop' ) 271 269 272 ! WRITE ( 9, * ) '*** advec_particles: ##0.3'273 ! CALL local_flush( 9 )274 ! nd = 0275 ! DO n = 1, number_of_particles276 ! IF ( .NOT. particle_mask(n) ) nd = nd + 1277 ! ENDDO278 ! IF ( nd /= deleted_particles ) THEN279 ! WRITE (9,*) '*** nd=',nd,' deleted_particles=',deleted_particles280 ! CALL local_flush( 9 )281 ! CALL MPI_ABORT( comm2d, 9999, ierr )282 ! ENDIF283 284 270 ! 285 271 !-- Particle (droplet) growth by condensation/evaporation and collision … … 357 343 ! 358 344 !-- Change in radius by condensation/evaporation 359 IF ( curvature_solution_effects .AND. particles(n)%radius < 1.0E-6 ) & 360 THEN 361 ! 362 !-- Curvature and solutions effects are included in growth equation. 363 !-- Change in Radius is calculated with 4th-order Rosenbrock method 364 !-- for stiff o.d.e's with monitoring local truncation error to adjust 365 !-- stepsize (see Numerical recipes in FORTRAN, 2nd edition, p. 731). 366 !-- For larger radii the simple analytic method (see ELSE) gives 367 !-- almost the same results. 368 ! 369 !-- Surface tension after (Straka, 2009) 370 sigma = 0.0761 - 0.00155 * ( t_int - 273.15 ) 371 372 r_ros = particles(n)%radius 373 dt_ros_sum = 0.0 ! internal integrated time (s) 374 internal_timestep_count = 0 375 376 ddenom = 1.0 / ( rho_l * r_v * t_int / ( e_s * diff_coeff_l ) + & 377 ( l_v / ( r_v * t_int ) - 1.0 ) * & 378 rho_l * l_v / ( thermal_conductivity_l * t_int )& 379 ) 380 381 afactor = 2.0 * sigma / ( rho_l * r_v * t_int ) 382 383 IF ( particles(n)%rvar3 == -9999999.9 ) THEN 384 ! 385 !-- First particle timestep. Derivative has to be calculated. 386 drdt_m = ddenom / r_ros * ( e_a / e_s - 1.0 - & 387 afactor / r_ros + & 388 bfactor / r_ros**3 ) 389 ELSE 390 ! 391 !-- Take value from last PALM timestep 392 drdt_m = particles(n)%rvar3 393 ENDIF 394 ! 395 !-- Take internal timestep values from the end of last PALM timestep 396 dt_ros_last = particles(n)%rvar1 397 dt_ros_next = particles(n)%rvar2 398 ! 399 !-- Internal timestep must not be larger than PALM timestep 400 dt_ros = MIN( dt_ros_next, dt_3d ) 401 ! 402 !-- Integrate growth equation in time unless PALM timestep is reached 403 DO WHILE ( dt_ros_sum < dt_3d ) 404 405 internal_timestep_count = internal_timestep_count + 1 406 407 ! 408 !-- Derivative at starting value 409 drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - afactor / r_ros + & 410 bfactor / r_ros**3 ) 411 drdt_ini = drdt 412 dt_ros_sum_ini = dt_ros_sum 413 r_ros_ini = r_ros 414 415 ! 416 !-- Calculate time derivative of dr/dt 417 d2rdt2 = ( drdt - drdt_m ) / dt_ros_last 418 419 ! 420 !-- Calculate radial derivative of dr/dt 421 d2rdtdr = ddenom * ( ( 1.0 - e_a/e_s ) / r_ros**2 + & 422 2.0 * afactor / r_ros**3 - & 423 4.0 * bfactor / r_ros**5 ) 424 ! 425 !-- Adjust stepsize unless required accuracy is reached 426 DO jtry = 1, maxtry+1 427 428 IF ( jtry == maxtry+1 ) THEN 429 message_string = 'maxtry > 40 in Rosenbrock method' 430 CALL message( 'advec_particles', 'PA0347', 2, 2, 0, 6, 0 ) 431 ENDIF 432 433 aa = 1.0 / ( gam * dt_ros ) - d2rdtdr 434 g1 = ( drdt_ini + dt_ros * c1x * d2rdt2 ) / aa 435 r_ros = r_ros_ini + a21 * g1 436 drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - & 437 afactor / r_ros + & 438 bfactor / r_ros**3 ) 439 440 g2 = ( drdt + dt_ros * c2x * d2rdt2 + c21 * g1 / dt_ros )& 441 / aa 442 r_ros = r_ros_ini + a31 * g1 + a32 * g2 443 drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - & 444 afactor / r_ros + & 445 bfactor / r_ros**3 ) 446 447 g3 = ( drdt + dt_ros * c3x * d2rdt2 + & 448 ( c31 * g1 + c32 * g2 ) / dt_ros ) / aa 449 g4 = ( drdt + dt_ros * c4x * d2rdt2 + & 450 ( c41 * g1 + c42 * g2 + c43 * g3 ) / dt_ros ) / aa 451 r_ros = r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + b4 * g4 452 453 dt_ros_sum = dt_ros_sum_ini + dt_ros 454 455 IF ( dt_ros_sum == dt_ros_sum_ini ) THEN 456 message_string = 'zero stepsize in Rosenbrock method' 457 CALL message( 'advec_particles', 'PA0348', 2, 2, 0, 6, 0 ) 458 ENDIF 459 ! 460 !-- Calculate error 461 err_ros = e1*g1 + e2*g2 + e3*g3 + e4*g4 462 errmax = 0.0 463 errmax = MAX( errmax, ABS( err_ros / r_ros_ini ) ) / eps_ros 464 ! 465 !-- Leave loop if accuracy is sufficient, otherwise try again 466 !-- with a reduced stepsize 467 IF ( errmax <= 1.0 ) THEN 468 EXIT 469 ELSE 470 dt_ros = SIGN( & 471 MAX( ABS( 0.9 * dt_ros * errmax**pshrnk ), & 472 shrnk * ABS( dt_ros ) ), & 473 dt_ros & 474 ) 475 ENDIF 476 477 ENDDO ! loop for stepsize adjustment 478 479 ! 480 !-- Calculate next internal timestep 481 IF ( errmax > errcon ) THEN 482 dt_ros_next = 0.9 * dt_ros * errmax**pgrow 483 ELSE 484 dt_ros_next = grow * dt_ros 485 ENDIF 486 487 ! 488 !-- Estimated timestep is reduced if the PALM time step is exceeded 489 dt_ros_last = dt_ros 490 IF ( ( dt_ros_next + dt_ros_sum ) >= dt_3d ) THEN 491 dt_ros = dt_3d - dt_ros_sum 492 ELSE 493 dt_ros = dt_ros_next 494 ENDIF 495 496 drdt_m = drdt 497 498 ENDDO 499 ! 500 !-- Store derivative and internal timestep values for next PALM step 501 particles(n)%rvar1 = dt_ros_last 502 particles(n)%rvar2 = dt_ros_next 503 particles(n)%rvar3 = drdt_m 504 505 new_r = r_ros 506 ! 507 !-- Radius should not fall below 1E-8 because Rosenbrock method may 508 !-- lead to errors otherwise 509 new_r = MAX( new_r, 1.0E-8 ) 510 511 ELSE 345 IF ( particles(n)%radius >= 1.0E-6 .OR. & 346 .NOT. curvature_solution_effects ) THEN 512 347 ! 513 348 !-- Approximation for large radii, where curvature and solution … … 523 358 new_r = SQRT( arg ) 524 359 ENDIF 360 ENDIF 361 362 IF ( curvature_solution_effects .AND. & 363 ( ( particles(n)%radius < 1.0E-6 ) .OR. ( new_r < 1.0E-6 ) ) ) & 364 THEN 365 ! 366 !-- Curvature and solutions effects are included in growth equation. 367 !-- Change in Radius is calculated with 4th-order Rosenbrock method 368 !-- for stiff o.d.e's with monitoring local truncation error to adjust 369 !-- stepsize (see Numerical recipes in FORTRAN, 2nd edition, p. 731). 370 !-- For larger radii the simple analytic method (see ELSE) gives 371 !-- almost the same results. 372 ! 373 !-- Surface tension after (Straka, 2009) 374 sigma = 0.0761 - 0.00155 * ( t_int - 273.15 ) 375 376 r_ros = particles(n)%radius 377 dt_ros_sum = 0.0 ! internal integrated time (s) 378 internal_timestep_count = 0 379 380 ddenom = 1.0 / ( rho_l * r_v * t_int / ( e_s * diff_coeff_l ) + & 381 ( l_v / ( r_v * t_int ) - 1.0 ) * & 382 rho_l * l_v / ( thermal_conductivity_l * t_int )& 383 ) 384 385 afactor = 2.0 * sigma / ( rho_l * r_v * t_int ) 386 387 IF ( particles(n)%rvar3 == -9999999.9 ) THEN 388 ! 389 !-- First particle timestep. Derivative has to be calculated. 390 drdt_m = ddenom / r_ros * ( e_a / e_s - 1.0 - & 391 afactor / r_ros + & 392 bfactor / r_ros**3 ) 393 ELSE 394 ! 395 !-- Take value from last PALM timestep 396 drdt_m = particles(n)%rvar3 397 ENDIF 398 ! 399 !-- Take internal timestep values from the end of last PALM timestep 400 dt_ros_last = particles(n)%rvar1 401 dt_ros_next = particles(n)%rvar2 402 ! 403 !-- Internal timestep must not be larger than PALM timestep 404 dt_ros = MIN( dt_ros_next, dt_3d ) 405 ! 406 !-- Integrate growth equation in time unless PALM timestep is reached 407 DO WHILE ( dt_ros_sum < dt_3d ) 408 409 internal_timestep_count = internal_timestep_count + 1 410 411 ! 412 !-- Derivative at starting value 413 drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - afactor / r_ros + & 414 bfactor / r_ros**3 ) 415 drdt_ini = drdt 416 dt_ros_sum_ini = dt_ros_sum 417 r_ros_ini = r_ros 418 419 ! 420 !-- Calculate time derivative of dr/dt 421 d2rdt2 = ( drdt - drdt_m ) / dt_ros_last 422 423 ! 424 !-- Calculate radial derivative of dr/dt 425 d2rdtdr = ddenom * ( ( 1.0 - e_a/e_s ) / r_ros**2 + & 426 2.0 * afactor / r_ros**3 - & 427 4.0 * bfactor / r_ros**5 ) 428 ! 429 !-- Adjust stepsize unless required accuracy is reached 430 DO jtry = 1, maxtry+1 431 432 IF ( jtry == maxtry+1 ) THEN 433 message_string = 'maxtry > 40 in Rosenbrock method' 434 CALL message( 'advec_particles', 'PA0347', 2, 2, 0, 6, 0 ) 435 ENDIF 436 437 aa = 1.0 / ( gam * dt_ros ) - d2rdtdr 438 g1 = ( drdt_ini + dt_ros * c1x * d2rdt2 ) / aa 439 r_ros = r_ros_ini + a21 * g1 440 drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - & 441 afactor / r_ros + & 442 bfactor / r_ros**3 ) 443 444 g2 = ( drdt + dt_ros * c2x * d2rdt2 + c21 * g1 / dt_ros )& 445 / aa 446 r_ros = r_ros_ini + a31 * g1 + a32 * g2 447 drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - & 448 afactor / r_ros + & 449 bfactor / r_ros**3 ) 450 451 g3 = ( drdt + dt_ros * c3x * d2rdt2 + & 452 ( c31 * g1 + c32 * g2 ) / dt_ros ) / aa 453 g4 = ( drdt + dt_ros * c4x * d2rdt2 + & 454 ( c41 * g1 + c42 * g2 + c43 * g3 ) / dt_ros ) / aa 455 r_ros = r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + b4 * g4 456 457 dt_ros_sum = dt_ros_sum_ini + dt_ros 458 459 IF ( dt_ros_sum == dt_ros_sum_ini ) THEN 460 message_string = 'zero stepsize in Rosenbrock method' 461 CALL message( 'advec_particles', 'PA0348', 2, 2, 0, 6, 0 ) 462 ENDIF 463 ! 464 !-- Calculate error 465 err_ros = e1*g1 + e2*g2 + e3*g3 + e4*g4 466 errmax = 0.0 467 errmax = MAX( errmax, ABS( err_ros / r_ros_ini ) ) / eps_ros 468 ! 469 !-- Leave loop if accuracy is sufficient, otherwise try again 470 !-- with a reduced stepsize 471 IF ( errmax <= 1.0 ) THEN 472 EXIT 473 ELSE 474 dt_ros = SIGN( & 475 MAX( ABS( 0.9 * dt_ros * errmax**pshrnk ), & 476 shrnk * ABS( dt_ros ) ), & 477 dt_ros & 478 ) 479 ENDIF 480 481 ENDDO ! loop for stepsize adjustment 482 483 ! 484 !-- Calculate next internal timestep 485 IF ( errmax > errcon ) THEN 486 dt_ros_next = 0.9 * dt_ros * errmax**pgrow 487 ELSE 488 dt_ros_next = grow * dt_ros 489 ENDIF 490 491 ! 492 !-- Estimated timestep is reduced if the PALM time step is exceeded 493 dt_ros_last = dt_ros 494 IF ( ( dt_ros_next + dt_ros_sum ) >= dt_3d ) THEN 495 dt_ros = dt_3d - dt_ros_sum 496 ELSE 497 dt_ros = dt_ros_next 498 ENDIF 499 500 drdt_m = drdt 501 502 ENDDO 503 ! 504 !-- Store derivative and internal timestep values for next PALM step 505 particles(n)%rvar1 = dt_ros_last 506 particles(n)%rvar2 = dt_ros_next 507 particles(n)%rvar3 = drdt_m 508 509 new_r = r_ros 510 ! 511 !-- Radius should not fall below 1E-8 because Rosenbrock method may 512 !-- lead to errors otherwise 513 new_r = MAX( new_r, 1.0E-8 ) 525 514 526 515 ENDIF … … 570 559 ! 571 560 !-- Particle growth by collision 561 IF ( collision_kernel /= 'none' ) THEN 562 572 563 CALL cpu_log( log_point_s(43), 'advec_part_coll', 'start' ) 573 574 IF ( wang_collision_kernel ) THEN575 564 576 565 DO i = nxl, nxr … … 609 598 pse = psi + prt_count(k,j,i)-1 610 599 611 ALLOCATE( kern(psi:pse,psi:pse) ) 612 613 ! 614 !-- Calculate collision kernel for all particles in grid box 615 CALL colker( i, j, k, kern ) 616 ! 617 !-- collison kernel is calculated in cm**3/s but needed in m**3/s 618 kern = kern * 1.0E-6 619 620 DO n = pse, psi+1, -1 621 622 integral = 0.0 623 lw_max = 0.0 624 625 ! 626 !-- Calculate growth of collector particle radius using all 627 !-- droplets smaller than current droplet 628 DO is = psi, n-1 629 630 integral = integral + & 631 ( particles(is)%radius**3 * kern(n,is) * & 632 particles(is)%weight_factor ) 633 ! 634 !-- Calculate volume of liquid water of the collected 635 !-- droplets which is the maximum liquid water available 636 !-- for droplet growth 637 lw_max = lw_max + ( particles(is)%radius**3 * & 638 particles(is)%weight_factor ) 600 IF ( hall_kernel .OR. wang_kernel ) THEN 601 602 ALLOCATE( kern(psi:pse,psi:pse) ) 603 604 ! 605 !-- Calculate collision kernel for all particles in grid box 606 CALL colker( i, j, k, kern ) 607 ! 608 !-- Collison kernel is calculated in cm**3/s but needed 609 !-- in m**3/s 610 kern = kern * 1.0E-6 ! to be moved to colker 611 612 DO n = pse, psi+1, -1 613 614 integral = 0.0 615 lw_max = 0.0 616 617 ! 618 !-- Calculate growth of collector particle radius using all 619 !-- droplets smaller than current droplet 620 DO is = psi, n-1 621 622 integral = integral + & 623 ( particles(is)%radius**3 * kern(n,is) *& 624 particles(is)%weight_factor ) 625 ! 626 !-- Calculate volume of liquid water of the collected 627 !-- droplets which is the maximum liquid water available 628 !-- for droplet growth 629 lw_max = lw_max + ( particles(is)%radius**3 * & 630 particles(is)%weight_factor ) 631 632 ENDDO 633 634 ! 635 !-- Change in radius of collector droplet due to collision 636 delta_r = 1.0 / ( 3.0 * particles(n)%radius**2 ) * & 637 integral * dt_3d * ddx * ddy / dz 638 639 ! 640 !-- Change in volume of collector droplet due to collision 641 delta_v = particles(n)%weight_factor & 642 * ( ( particles(n)%radius + delta_r )**3 & 643 - particles(n)%radius**3 ) 644 645 IF ( lw_max < delta_v .AND. delta_v > 0.0 ) THEN 646 !-- replace by message call 647 PRINT*, 'Particle has grown to much because the', & 648 ' change of volume of particle is larger', & 649 ' than liquid water available!' 650 651 ELSEIF ( lw_max == delta_v .AND. delta_v > 0.0 ) THEN 652 653 DO is = psi, n-1 654 655 particles(is)%weight_factor = 0.0 656 particle_mask(is) = .FALSE. 657 deleted_particles = deleted_particles + 1 658 659 ENDDO 660 661 ELSEIF ( lw_max > delta_v .AND. delta_v > 0.0 ) THEN 662 ! 663 !-- Calculate new weighting factor of collected droplets 664 DO is = psi, n-1 665 666 particles(is)%weight_factor = & 667 particles(is)%weight_factor - & 668 ( ( kern(n,is) * particles(is)%weight_factor & 669 / integral ) * delta_v ) 670 671 IF ( particles(is)%weight_factor < 0.0 ) THEN 672 WRITE( message_string, * ) 'negative ', & 673 'weighting factor: ', & 674 particles(is)%weight_factor 675 CALL message( 'advec_particles', '', 2, 2,& 676 -1, 6, 1 ) 677 678 ELSEIF ( particles(is)%weight_factor == 0.0 ) & 679 THEN 680 681 particles(is)%weight_factor = 0.0 682 particle_mask(is) = .FALSE. 683 deleted_particles = deleted_particles + 1 684 685 ENDIF 686 687 ENDDO 688 689 ENDIF 690 691 particles(n)%radius = particles(n)%radius + delta_r 692 ql_vp(k,j,i) = ql_vp(k,j,i) + & 693 particles(n)%weight_factor & 694 * particles(n)%radius**3 639 695 640 696 ENDDO 641 697 642 ! 643 !-- Change in radius of collector droplet due to collision 644 delta_r = 1.0 / ( 3.0 * particles(n)%radius**2 ) * & 645 integral * dt_3d * ddx * ddy / dz 646 647 ! 648 !-- Change in volume of collector droplet due to collision 649 delta_v = particles(n)%weight_factor & 650 * ( ( particles(n)%radius + delta_r )**3 & 651 - particles(n)%radius**3 ) 652 653 IF ( lw_max < delta_v .AND. delta_v > 0.0 ) THEN 654 655 PRINT*, 'Particle has grown to much because the & 656 change of volume of particle is larger than & 657 liquid water available!' 658 659 ELSEIF ( lw_max == delta_v .AND. delta_v > 0.0 ) THEN 660 661 DO is = psi, n-1 662 663 particles(is)%weight_factor = 0.0 664 particle_mask(is) = .FALSE. 665 deleted_particles = deleted_particles + 1 666 698 DEALLOCATE( kern ) 699 700 ELSEIF ( palm_kernel ) THEN 701 ! 702 !-- PALM collision kernel 703 ! 704 !-- Calculate the mean radius of all those particles which 705 !-- are of smaller size than the current particle and 706 !-- use this radius for calculating the collision efficiency 707 DO n = psi+prt_count(k,j,i)-1, psi+1, -1 708 709 sl_r3 = 0.0 710 sl_r4 = 0.0 711 712 DO is = n-1, psi, -1 713 IF ( particles(is)%radius < particles(n)%radius ) & 714 THEN 715 sl_r3 = sl_r3 + particles(is)%weight_factor & 716 *( particles(is)%radius**3 ) 717 sl_r4 = sl_r4 + particles(is)%weight_factor & 718 *( particles(is)%radius**4 ) 719 ENDIF 667 720 ENDDO 668 721 669 ELSEIF ( lw_max > delta_v .AND. delta_v > 0.0 ) THEN 670 ! 671 !-- Calculate new weighting factor of collected droplets 672 DO is = psi, n-1 673 674 particles(is)%weight_factor = & 675 particles(is)%weight_factor - & 676 ( ( kern(n,is) * particles(is)%weight_factor / & 677 integral ) * delta_v ) 678 679 IF ( particles(is)%weight_factor < 0.0 ) THEN 680 WRITE( message_string, * ) 'negative ', & 722 IF ( ( sl_r3 ) > 0.0 ) THEN 723 mean_r = ( sl_r4 ) / ( sl_r3 ) 724 725 CALL collision_efficiency( mean_r, & 726 particles(n)%radius, & 727 effective_coll_efficiency ) 728 729 ELSE 730 effective_coll_efficiency = 0.0 731 ENDIF 732 733 IF ( effective_coll_efficiency > 1.0 .OR. & 734 effective_coll_efficiency < 0.0 ) & 735 THEN 736 WRITE( message_string, * ) 'collision_efficien' , & 737 'cy out of range:' ,effective_coll_efficiency 738 CALL message( 'advec_particles', 'PA0145', 2, 2, & 739 -1, 6, 1 ) 740 ENDIF 741 742 ! 743 !-- Interpolation of ... 744 ii = particles(n)%x * ddx 745 jj = particles(n)%y * ddy 746 kk = ( particles(n)%z + 0.5 * dz ) / dz 747 748 x = particles(n)%x - ii * dx 749 y = particles(n)%y - jj * dy 750 aa = x**2 + y**2 751 bb = ( dx - x )**2 + y**2 752 cc = x**2 + ( dy - y )**2 753 dd = ( dx - x )**2 + ( dy - y )**2 754 gg = aa + bb + cc + dd 755 756 ql_int_l = ( (gg-aa) * ql(kk,jj,ii) + (gg-bb) * & 757 ql(kk,jj,ii+1) & 758 + (gg-cc) * ql(kk,jj+1,ii) + ( gg-dd ) * & 759 ql(kk,jj+1,ii+1) & 760 ) / ( 3.0 * gg ) 761 762 ql_int_u = ( (gg-aa) * ql(kk+1,jj,ii) + (gg-bb) * & 763 ql(kk+1,jj,ii+1) & 764 + (gg-cc) * ql(kk+1,jj+1,ii) + (gg-dd) * & 765 ql(kk+1,jj+1,ii+1) & 766 ) / ( 3.0 * gg ) 767 768 ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz *& 769 ( ql_int_u - ql_int_l ) 770 771 ! 772 !-- Interpolate u velocity-component 773 ii = ( particles(n)%x + 0.5 * dx ) * ddx 774 jj = particles(n)%y * ddy 775 kk = ( particles(n)%z + 0.5 * dz ) / dz ! only if eqist 776 777 IF ( ( particles(n)%z - zu(kk) ) > (0.5*dz) ) kk = kk+1 778 779 x = particles(n)%x + ( 0.5 - ii ) * dx 780 y = particles(n)%y - jj * dy 781 aa = x**2 + y**2 782 bb = ( dx - x )**2 + y**2 783 cc = x**2 + ( dy - y )**2 784 dd = ( dx - x )**2 + ( dy - y )**2 785 gg = aa + bb + cc + dd 786 787 u_int_l = ( (gg-aa) * u(kk,jj,ii) + (gg-bb) * & 788 u(kk,jj,ii+1) & 789 + (gg-cc) * u(kk,jj+1,ii) + (gg-dd) * & 790 u(kk,jj+1,ii+1) & 791 ) / ( 3.0 * gg ) - u_gtrans 792 IF ( kk+1 == nzt+1 ) THEN 793 u_int = u_int_l 794 ELSE 795 u_int_u = ( (gg-aa) * u(kk+1,jj,ii) + (gg-bb) * & 796 u(kk+1,jj,ii+1) & 797 + (gg-cc) * u(kk+1,jj+1,ii) + (gg-dd) * & 798 u(kk+1,jj+1,ii+1) & 799 ) / ( 3.0 * gg ) - u_gtrans 800 u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz & 801 * ( u_int_u - u_int_l ) 802 ENDIF 803 804 ! 805 !-- Same procedure for interpolation of the v velocity-com- 806 !-- ponent (adopt index k from u velocity-component) 807 ii = particles(n)%x * ddx 808 jj = ( particles(n)%y + 0.5 * dy ) * ddy 809 810 x = particles(n)%x - ii * dx 811 y = particles(n)%y + ( 0.5 - jj ) * dy 812 aa = x**2 + y**2 813 bb = ( dx - x )**2 + y**2 814 cc = x**2 + ( dy - y )**2 815 dd = ( dx - x )**2 + ( dy - y )**2 816 gg = aa + bb + cc + dd 817 818 v_int_l = ( ( gg-aa ) * v(kk,jj,ii) + ( gg-bb ) * & 819 v(kk,jj,ii+1) & 820 + ( gg-cc ) * v(kk,jj+1,ii) + ( gg-dd ) * & 821 v(kk,jj+1,ii+1) & 822 ) / ( 3.0 * gg ) - v_gtrans 823 IF ( kk+1 == nzt+1 ) THEN 824 v_int = v_int_l 825 ELSE 826 v_int_u = ( (gg-aa) * v(kk+1,jj,ii) + (gg-bb) * & 827 v(kk+1,jj,ii+1) & 828 + (gg-cc) * v(kk+1,jj+1,ii) + (gg-dd) * & 829 v(kk+1,jj+1,ii+1) & 830 ) / ( 3.0 * gg ) - v_gtrans 831 v_int = v_int_l + ( particles(n)%z - zu(kk) ) / dz & 832 * ( v_int_u - v_int_l ) 833 ENDIF 834 835 ! 836 !-- Same procedure for interpolation of the w velocity-com- 837 !-- ponent (adopt index i from v velocity-component) 838 jj = particles(n)%y * ddy 839 kk = particles(n)%z / dz 840 841 x = particles(n)%x - ii * dx 842 y = particles(n)%y - jj * dy 843 aa = x**2 + y**2 844 bb = ( dx - x )**2 + y**2 845 cc = x**2 + ( dy - y )**2 846 dd = ( dx - x )**2 + ( dy - y )**2 847 gg = aa + bb + cc + dd 848 849 w_int_l = ( ( gg-aa ) * w(kk,jj,ii) + ( gg-bb ) * & 850 w(kk,jj,ii+1) & 851 + ( gg-cc ) * w(kk,jj+1,ii) + ( gg-dd ) * & 852 w(kk,jj+1,ii+1) & 853 ) / ( 3.0 * gg ) 854 IF ( kk+1 == nzt+1 ) THEN 855 w_int = w_int_l 856 ELSE 857 w_int_u = ( (gg-aa) * w(kk+1,jj,ii) + (gg-bb) * & 858 w(kk+1,jj,ii+1) & 859 + (gg-cc) * w(kk+1,jj+1,ii) + (gg-dd) * & 860 w(kk+1,jj+1,ii+1) & 861 ) / ( 3.0 * gg ) 862 w_int = w_int_l + ( particles(n)%z - zw(kk) ) / dz & 863 * ( w_int_u - w_int_l ) 864 ENDIF 865 866 ! 867 !-- Change in radius due to collision 868 delta_r = effective_coll_efficiency / 3.0 & 869 * pi * sl_r3 * ddx * ddy / dz & 870 * SQRT( ( u_int - particles(n)%speed_x )**2 & 871 + ( v_int - particles(n)%speed_y )**2 & 872 + ( w_int - particles(n)%speed_z )**2 & 873 ) * dt_3d 874 ! 875 !-- Change in volume due to collision 876 delta_v = particles(n)%weight_factor & 877 * ( ( particles(n)%radius + delta_r )**3 & 878 - particles(n)%radius**3 ) 879 880 ! 881 !-- Check if collected particles provide enough LWC for 882 !-- volume change of collector particle 883 IF ( delta_v >= sl_r3 .AND. sl_r3 > 0.0 ) THEN 884 885 delta_r = ( ( sl_r3/particles(n)%weight_factor ) & 886 + particles(n)%radius**3 )**( 1./3. ) & 887 - particles(n)%radius 888 889 DO is = n-1, psi, -1 890 IF ( particles(is)%radius < & 891 particles(n)%radius ) THEN 892 particles(is)%weight_factor = 0.0 893 particle_mask(is) = .FALSE. 894 deleted_particles = deleted_particles + 1 895 ENDIF 896 ENDDO 897 898 ELSE IF ( delta_v < sl_r3 .AND. sl_r3 > 0.0 ) THEN 899 900 DO is = n-1, psi, -1 901 IF ( particles(is)%radius < particles(n)%radius & 902 .AND. sl_r3 > 0.0 ) THEN 903 particles(is)%weight_factor = & 904 ( ( particles(is)%weight_factor & 905 * ( particles(is)%radius**3 ) ) & 906 - ( delta_v & 907 * particles(is)%weight_factor & 908 * ( particles(is)%radius**3 ) & 909 / sl_r3 ) ) & 910 / ( particles(is)%radius**3 ) 911 912 IF ( particles(is)%weight_factor < 0.0 ) THEN 913 WRITE( message_string, * ) 'negative ', & 681 914 'weighting factor: ', & 682 915 particles(is)%weight_factor 683 CALL message( 'advec_particles', '', 2, 2, & 684 -1, 6, 1 ) 685 686 ELSEIF ( particles(is)%weight_factor == 0.0 ) THEN 687 688 particles(is)%weight_factor = 0.0 689 particle_mask(is) = .FALSE. 690 deleted_particles = deleted_particles + 1 691 692 ENDIF 693 694 ENDDO 695 696 ENDIF 697 698 particles(n)%radius = particles(n)%radius + delta_r 699 ql_vp(k,j,i) = ql_vp(k,j,i) + particles(n)%weight_factor & 700 * particles(n)%radius**3 701 702 ENDDO 916 CALL message( 'advec_particles', '', 2, & 917 2, -1, 6, 1 ) 918 ENDIF 919 ENDIF 920 ENDDO 921 ENDIF 922 923 particles(n)%radius = particles(n)%radius + delta_r 924 ql_vp(k,j,i) = ql_vp(k,j,i) + & 925 particles(n)%weight_factor * & 926 ( particles(n)%radius**3 ) 927 ENDDO 928 929 ENDIF ! collision kernel 703 930 704 931 ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor & 705 932 * particles(psi)%radius**3 706 933 707 DEALLOCATE( kern )708 934 709 935 ELSE IF ( prt_count(k,j,i) == 1 ) THEN … … 735 961 ENDDO 736 962 737 ELSE 738 739 DO i = nxl, nxr 740 DO j = nys, nyn 741 DO k = nzb+1, nzt 742 ! 743 !-- Collision requires at least two particles in the box 744 IF ( prt_count(k,j,i) > 1 ) THEN 745 ! 746 !-- First, sort particles within the gridbox by their size, 747 !-- using Shell's method (see Numerical Recipes) 748 !-- NOTE: In case of using particle tails, the re-sorting of 749 !-- ---- tails would have to be included here! 750 psi = prt_start_index(k,j,i) - 1 751 inc = 1 752 DO WHILE ( inc <= prt_count(k,j,i) ) 753 inc = 3 * inc + 1 754 ENDDO 755 756 DO WHILE ( inc > 1 ) 757 inc = inc / 3 758 DO is = inc+1, prt_count(k,j,i) 759 tmp_particle = particles(psi+is) 760 js = is 761 DO WHILE ( particles(psi+js-inc)%radius > & 762 tmp_particle%radius ) 763 particles(psi+js) = particles(psi+js-inc) 764 js = js - inc 765 IF ( js <= inc ) EXIT 766 ENDDO 767 particles(psi+js) = tmp_particle 768 ENDDO 769 ENDDO 770 771 ! 772 !-- Calculate the mean radius of all those particles which 773 !-- are of smaller size than the current particle 774 !-- and use this radius for calculating the collision efficiency 775 psi = prt_start_index(k,j,i) 776 777 DO n = psi+prt_count(k,j,i)-1, psi+1, -1 778 sl_r3 = 0.0 779 sl_r4 = 0.0 780 781 DO is = n-1, psi, -1 782 IF ( particles(is)%radius < particles(n)%radius ) THEN 783 sl_r3 = sl_r3 + particles(is)%weight_factor & 784 *( particles(is)%radius**3 ) 785 sl_r4 = sl_r4 + particles(is)%weight_factor & 786 *( particles(is)%radius**4 ) 787 ENDIF 788 ENDDO 789 790 IF ( ( sl_r3 ) > 0.0 ) THEN 791 mean_r = ( sl_r4 ) / ( sl_r3 ) 792 793 CALL collision_efficiency( mean_r, & 794 particles(n)%radius, & 795 effective_coll_efficiency ) 796 797 ELSE 798 effective_coll_efficiency = 0.0 799 ENDIF 800 801 IF ( effective_coll_efficiency > 1.0 .OR. & 802 effective_coll_efficiency < 0.0 ) & 803 THEN 804 WRITE( message_string, * ) 'collision_efficiency ' , & 805 'out of range:' ,effective_coll_efficiency 806 CALL message( 'advec_particles', 'PA0145', 2, 2, -1, & 807 6, 1 ) 808 ENDIF 809 810 ! 811 !-- Interpolation of ... 812 ii = particles(n)%x * ddx 813 jj = particles(n)%y * ddy 814 kk = ( particles(n)%z + 0.5 * dz ) / dz 815 816 x = particles(n)%x - ii * dx 817 y = particles(n)%y - jj * dy 818 aa = x**2 + y**2 819 bb = ( dx - x )**2 + y**2 820 cc = x**2 + ( dy - y )**2 821 dd = ( dx - x )**2 + ( dy - y )**2 822 gg = aa + bb + cc + dd 823 824 ql_int_l = ( ( gg-aa ) * ql(kk,jj,ii) + ( gg-bb ) * & 825 ql(kk,jj,ii+1) & 826 + ( gg-cc ) * ql(kk,jj+1,ii) + ( gg-dd ) * & 827 ql(kk,jj+1,ii+1) & 828 ) / ( 3.0 * gg ) 829 830 ql_int_u = ( ( gg-aa ) * ql(kk+1,jj,ii) + ( gg-bb ) * & 831 ql(kk+1,jj,ii+1) & 832 + ( gg-cc ) * ql(kk+1,jj+1,ii) + ( gg-dd ) * & 833 ql(kk+1,jj+1,ii+1) & 834 ) / ( 3.0 * gg ) 835 836 ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz * & 837 ( ql_int_u - ql_int_l ) 838 839 ! 840 !-- Interpolate u velocity-component 841 ii = ( particles(n)%x + 0.5 * dx ) * ddx 842 jj = particles(n)%y * ddy 843 kk = ( particles(n)%z + 0.5 * dz ) / dz ! only if eq.dist 844 845 IF ( ( particles(n)%z - zu(kk) ) > ( 0.5*dz ) ) kk = kk+1 846 847 x = particles(n)%x + ( 0.5 - ii ) * dx 848 y = particles(n)%y - jj * dy 849 aa = x**2 + y**2 850 bb = ( dx - x )**2 + y**2 851 cc = x**2 + ( dy - y )**2 852 dd = ( dx - x )**2 + ( dy - y )**2 853 gg = aa + bb + cc + dd 854 855 u_int_l = ( ( gg-aa ) * u(kk,jj,ii) + ( gg-bb ) * & 856 u(kk,jj,ii+1) & 857 + ( gg-cc ) * u(kk,jj+1,ii) + ( gg-dd ) * & 858 u(kk,jj+1,ii+1) & 859 ) / ( 3.0 * gg ) - u_gtrans 860 IF ( kk+1 == nzt+1 ) THEN 861 u_int = u_int_l 862 ELSE 863 u_int_u = ( ( gg-aa ) * u(kk+1,jj,ii) + ( gg-bb ) * & 864 u(kk+1,jj,ii+1) & 865 + ( gg-cc ) * u(kk+1,jj+1,ii) + ( gg-dd ) * & 866 u(kk+1,jj+1,ii+1) & 867 ) / ( 3.0 * gg ) - u_gtrans 868 u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz * & 869 ( u_int_u - u_int_l ) 870 ENDIF 871 872 ! 873 !-- Same procedure for interpolation of the v velocity-compo- 874 !-- nent (adopt index k from u velocity-component) 875 ii = particles(n)%x * ddx 876 jj = ( particles(n)%y + 0.5 * dy ) * ddy 877 878 x = particles(n)%x - ii * dx 879 y = particles(n)%y + ( 0.5 - jj ) * dy 880 aa = x**2 + y**2 881 bb = ( dx - x )**2 + y**2 882 cc = x**2 + ( dy - y )**2 883 dd = ( dx - x )**2 + ( dy - y )**2 884 gg = aa + bb + cc + dd 885 886 v_int_l = ( ( gg-aa ) * v(kk,jj,ii) + ( gg-bb ) * & 887 v(kk,jj,ii+1) & 888 + ( gg-cc ) * v(kk,jj+1,ii) + ( gg-dd ) * & 889 v(kk,jj+1,ii+1) & 890 ) / ( 3.0 * gg ) - v_gtrans 891 IF ( kk+1 == nzt+1 ) THEN 892 v_int = v_int_l 893 ELSE 894 v_int_u = ( ( gg-aa ) * v(kk+1,jj,ii) + ( gg-bb ) * & 895 v(kk+1,jj,ii+1) & 896 + ( gg-cc ) * v(kk+1,jj+1,ii) + ( gg-dd ) * & 897 v(kk+1,jj+1,ii+1) & 898 ) / ( 3.0 * gg ) - v_gtrans 899 v_int = v_int_l + ( particles(n)%z - zu(kk) ) / dz * & 900 ( v_int_u - v_int_l ) 901 ENDIF 902 903 ! 904 !-- Same procedure for interpolation of the w velocity-compo- 905 !-- nent (adopt index i from v velocity-component) 906 jj = particles(n)%y * ddy 907 kk = particles(n)%z / dz 908 909 x = particles(n)%x - ii * dx 910 y = particles(n)%y - jj * dy 911 aa = x**2 + y**2 912 bb = ( dx - x )**2 + y**2 913 cc = x**2 + ( dy - y )**2 914 dd = ( dx - x )**2 + ( dy - y )**2 915 gg = aa + bb + cc + dd 916 917 w_int_l = ( ( gg-aa ) * w(kk,jj,ii) + ( gg-bb ) * & 918 w(kk,jj,ii+1) & 919 + ( gg-cc ) * w(kk,jj+1,ii) + ( gg-dd ) * & 920 w(kk,jj+1,ii+1) & 921 ) / ( 3.0 * gg ) 922 IF ( kk+1 == nzt+1 ) THEN 923 w_int = w_int_l 924 ELSE 925 w_int_u = ( ( gg-aa ) * w(kk+1,jj,ii) + ( gg-bb ) * & 926 w(kk+1,jj,ii+1) & 927 + ( gg-cc ) * w(kk+1,jj+1,ii) + ( gg-dd ) * & 928 w(kk+1,jj+1,ii+1) & 929 ) / ( 3.0 * gg ) 930 w_int = w_int_l + ( particles(n)%z - zw(kk) ) / dz * & 931 ( w_int_u - w_int_l ) 932 ENDIF 933 934 ! 935 !-- Change in radius due to collision 936 delta_r = effective_coll_efficiency / 3.0 & 937 * pi * sl_r3 * ddx * ddy / dz & 938 * SQRT( ( u_int - particles(n)%speed_x )**2 & 939 + ( v_int - particles(n)%speed_y )**2 & 940 + ( w_int - particles(n)%speed_z )**2 & 941 ) * dt_3d 942 ! 943 !-- Change in volume due to collision 944 delta_v = particles(n)%weight_factor & 945 * ( ( particles(n)%radius + delta_r )**3 & 946 - particles(n)%radius**3 ) 947 948 ! 949 !-- Check if collected particles provide enough LWC for volume 950 !-- change of collector particle 951 IF ( delta_v >= sl_r3 .and. sl_r3 > 0.0 ) THEN 952 953 delta_r = ( ( sl_r3/particles(n)%weight_factor ) & 954 + particles(n)%radius**3 )**( 1./3. ) & 955 - particles(n)%radius 956 957 DO is = n-1, psi, -1 958 IF ( particles(is)%radius < particles(n)%radius ) & 959 THEN 960 particles(is)%weight_factor = 0.0 961 particle_mask(is) = .FALSE. 962 deleted_particles = deleted_particles + 1 963 ENDIF 964 ENDDO 965 966 ELSE IF ( delta_v < sl_r3 .and. sl_r3 > 0.0 ) THEN 967 968 DO is = n-1, psi, -1 969 IF ( particles(is)%radius < particles(n)%radius & 970 .and. sl_r3 > 0.0 ) THEN 971 particles(is)%weight_factor = & 972 ( ( particles(is)%weight_factor & 973 * ( particles(is)%radius**3 ) ) & 974 - ( delta_v & 975 * particles(is)%weight_factor & 976 * ( particles(is)%radius**3 ) & 977 / sl_r3 ) ) & 978 / ( particles(is)%radius**3 ) 979 980 IF (particles(is)%weight_factor < 0.0) THEN 981 WRITE( message_string, * ) 'negative ', & 982 'weighting factor: ', & 983 particles(is)%weight_factor 984 CALL message( 'advec_particles', '', 2, 2, & 985 -1, 6, 1 ) 986 ENDIF 987 ENDIF 988 ENDDO 989 ENDIF 990 991 particles(n)%radius = particles(n)%radius + delta_r 992 ql_vp(k,j,i) = ql_vp(k,j,i) + particles(n)%weight_factor & 993 *( particles(n)%radius**3 ) 994 ENDDO 995 996 ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor & 997 *( particles(psi)%radius**3 ) 998 999 ELSE IF ( prt_count(k,j,i) == 1 ) THEN 1000 psi = prt_start_index(k,j,i) 1001 ql_vp(k,j,i) = particles(psi)%weight_factor & 1002 *( particles(psi)%radius**3 ) 1003 ENDIF 1004 ! 1005 !-- Check if condensation of LWC was conserved 1006 !-- during collision process 1007 IF (ql_v(k,j,i) /= 0.0 ) THEN 1008 IF (ql_vp(k,j,i)/ql_v(k,j,i) >= 1.00001 .or. & 1009 ql_vp(k,j,i)/ql_v(k,j,i) <= 0.99999 ) THEN 1010 WRITE( message_string, * ) 'LWC is not conserved during',& 1011 ' collision! ', & 1012 'LWC after condensation: ', & 1013 ql_v(k,j,i), & 1014 ' LWC after collision: ', & 1015 ql_vp(k,j,i) 1016 CALL message( 'advec_particles', '', 2, 2, -1, 6, 1 ) 1017 ENDIF 1018 ENDIF 1019 1020 ENDDO 1021 ENDDO 1022 ENDDO 1023 1024 ENDIF 963 ENDIF ! collision handling 1025 964 1026 965 CALL cpu_log( log_point_s(43), 'advec_part_coll', 'stop' ) 1027 966 1028 ENDIF 967 ENDIF ! cloud droplet handling 1029 968 1030 969 … … 2272 2211 particles(n)%speed_x = particles(n)%speed_x * exp_term + & 2273 2212 u_int * ( 1.0 - exp_term ) 2274 2213 particles(n)%speed_y = particles(n)%speed_y * exp_term + & 2275 2214 v_int * ( 1.0 - exp_term ) 2276 2277 2278 2215 particles(n)%speed_z = particles(n)%speed_z * exp_term + & 2216 ( w_int - ( 1.0 - dens_ratio ) * g / exp_arg ) & 2217 * ( 1.0 - exp_term ) 2279 2218 ENDIF 2280 2219 -
palm/trunk/SOURCE/check_parameters.f90
r824 r825 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! check for c urvature_solution_effects6 ! check for collision_kernel and curvature_solution_effects 7 7 ! 8 8 ! Former revisions: … … 709 709 ENDIF 710 710 711 ! 712 !-- Collision kernels: 713 SELECT CASE ( TRIM( collision_kernel ) ) 714 715 CASE ( 'hall' ) 716 hall_kernel = .TRUE. 717 718 CASE ( 'palm' ) 719 palm_kernel = .TRUE. 720 721 CASE ( 'wang' ) 722 wang_kernel = .TRUE. 723 724 CASE ( 'none' ) 725 726 727 CASE DEFAULT 728 message_string = 'unknown collision kernel: collision_kernel = "' // & 729 TRIM( collision_kernel ) // '"' 730 CALL message( 'check_parameters', 'PA0350', 1, 2, 0, 6, 0 ) 731 732 END SELECT 733 734 711 735 IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. & 712 736 TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN -
palm/trunk/SOURCE/data_output_ptseries.f90
r824 r825 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! mean particle radius added as output quantity.6 ! mean/minimum/maximum particle radius added as output quantity. 7 7 ! particle attributes speed_x|y|z_sgs renamed rvar1|2|3. 8 8 ! … … 35 35 !------------------------------------------------------------------------------! 36 36 37 USE cloud_parameters 37 38 USE control_parameters 38 39 USE cpulog … … 48 49 INTEGER :: i, inum, j, n 49 50 50 REAL, DIMENSION(0:number_of_particle_groups,30) :: pts_value, pts_value_l 51 51 REAL, DIMENSION(:,:), ALLOCATABLE :: pts_value, pts_value_l 52 52 53 53 … … 70 70 ENDIF 71 71 72 ALLOCATE( pts_value(0:number_of_particle_groups,dopts_num), & 73 pts_value_l(0:number_of_particle_groups,dopts_num) ) 74 72 75 pts_value_l = 0.0 76 pts_value_l(:,16) = 9999999.9 ! for calculation of minimum radius 73 77 74 78 ! … … 88 92 pts_value_l(0,7) = pts_value_l(0,7) + particles(n)%speed_y ! mean v 89 93 pts_value_l(0,8) = pts_value_l(0,8) + particles(n)%speed_z ! mean w 90 pts_value_l(0,9) = pts_value_l(0,9) + particles(n)%rvar1 ! mean sgsu 91 pts_value_l(0,10) = pts_value_l(0,10) + particles(n)%rvar2 ! mean sgsv 92 pts_value_l(0,11) = pts_value_l(0,11) + particles(n)%rvar3 ! mean sgsw 94 IF ( .NOT. curvature_solution_effects ) THEN 95 pts_value_l(0,9) = pts_value_l(0,9) + particles(n)%rvar1 ! mean sgsu 96 pts_value_l(0,10) = pts_value_l(0,10) + particles(n)%rvar2 ! mean sgsv 97 pts_value_l(0,11) = pts_value_l(0,11) + particles(n)%rvar3 ! mean sgsw 98 ENDIF 93 99 IF ( particles(n)%speed_z > 0.0 ) THEN 94 100 pts_value_l(0,12) = pts_value_l(0,12) + 1.0 ! # of upward moving prts … … 99 105 particles(n)%speed_z ! mean w down 100 106 ENDIF 101 pts_value_l(0,15) = pts_value_l(0,15) + particles(n)%radius 102 pts_value_l(0,16) = number_of_particles 103 pts_value_l(0,17) = number_of_particles 107 pts_value_l(0,15) = pts_value_l(0,15) + particles(n)%radius ! mean rad 108 pts_value_l(0,16) = MIN( pts_value_l(0,16), particles(n)%radius ) ! minrad 109 pts_value_l(0,17) = MAX( pts_value_l(0,17), particles(n)%radius ) ! maxrad 110 pts_value_l(0,18) = number_of_particles 111 pts_value_l(0,19) = number_of_particles 104 112 105 113 ! … … 119 127 pts_value_l(j,7) = pts_value_l(j,7) + particles(n)%speed_y 120 128 pts_value_l(j,8) = pts_value_l(j,8) + particles(n)%speed_z 121 pts_value_l(j,9) = pts_value_l(j,9) + particles(n)%rvar1 122 pts_value_l(j,10) = pts_value_l(j,10) + particles(n)%rvar2 123 pts_value_l(j,11) = pts_value_l(j,11) + particles(n)%rvar3 129 IF ( .NOT. curvature_solution_effects ) THEN 130 pts_value_l(j,9) = pts_value_l(j,9) + particles(n)%rvar1 131 pts_value_l(j,10) = pts_value_l(j,10) + particles(n)%rvar2 132 pts_value_l(j,11) = pts_value_l(j,11) + particles(n)%rvar3 133 ENDIF 124 134 IF ( particles(n)%speed_z > 0.0 ) THEN 125 135 pts_value_l(j,12) = pts_value_l(j,12) + 1.0 … … 129 139 ENDIF 130 140 pts_value_l(j,15) = pts_value_l(j,15) + particles(n)%radius 131 pts_value_l(j,16) = pts_value_l(j,16) + 1.0 132 pts_value_l(j,17) = pts_value_l(j,17) + 1.0 141 pts_value_l(j,16) = MIN( pts_value(j,16), particles(n)%radius ) 142 pts_value_l(j,17) = MAX( pts_value(j,17), particles(n)%radius ) 143 pts_value_l(j,18) = pts_value_l(j,18) + 1.0 144 pts_value_l(j,19) = pts_value_l(j,19) + 1.0 133 145 134 146 ENDIF … … 146 158 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 147 159 CALL MPI_ALLREDUCE( pts_value_l(0,16), pts_value(0,16), inum, MPI_REAL, & 160 MPI_MIN, comm2d, ierr ) 161 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 162 CALL MPI_ALLREDUCE( pts_value_l(0,17), pts_value(0,17), inum, MPI_REAL, & 148 163 MPI_MAX, comm2d, ierr ) 149 164 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 150 CALL MPI_ALLREDUCE( pts_value_l(0,17), pts_value(0,17), inum, MPI_REAL, & 165 CALL MPI_ALLREDUCE( pts_value_l(0,18), pts_value(0,18), inum, MPI_REAL, & 166 MPI_MAX, comm2d, ierr ) 167 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 168 CALL MPI_ALLREDUCE( pts_value_l(0,19), pts_value(0,19), inum, MPI_REAL, & 151 169 MPI_MIN, comm2d, ierr ) 152 170 #else 153 pts_value(:,1:1 7) = pts_value_l(:,1:17)171 pts_value(:,1:19) = pts_value_l(:,1:19) 154 172 #endif 155 173 156 174 ! 157 !-- Normalize the above calculated quantities with the total number of158 !-- particles175 !-- Normalize the above calculated quantities (except min/max values) with the 176 !-- total number of particles 159 177 IF ( number_of_particle_groups > 1 ) THEN 160 178 inum = number_of_particle_groups … … 186 204 DO n = 1, number_of_particles 187 205 188 pts_value_l(0, 18) = pts_value_l(0,18) + ( particles(n)%x - &206 pts_value_l(0,20) = pts_value_l(0,20) + ( particles(n)%x - & 189 207 particles(n)%origin_x - pts_value(0,2) )**2 ! x*2 190 pts_value_l(0, 19) = pts_value_l(0,19) + ( particles(n)%y - &208 pts_value_l(0,21) = pts_value_l(0,21) + ( particles(n)%y - & 191 209 particles(n)%origin_y - pts_value(0,3) )**2 ! y*2 192 pts_value_l(0,2 0) = pts_value_l(0,20) + ( particles(n)%z - &210 pts_value_l(0,22) = pts_value_l(0,22) + ( particles(n)%z - & 193 211 particles(n)%origin_z - pts_value(0,4) )**2 ! z*2 194 pts_value_l(0,21) = pts_value_l(0,21) + ( particles(n)%speed_x - & 195 pts_value(0,6) )**2 ! u*2 196 pts_value_l(0,22) = pts_value_l(0,22) + ( particles(n)%speed_y - & 197 pts_value(0,7) )**2 ! v*2 198 pts_value_l(0,23) = pts_value_l(0,23) + ( particles(n)%speed_z - & 199 pts_value(0,8) )**2 ! w*2 200 pts_value_l(0,24) = pts_value_l(0,24) + ( particles(n)%rvar1 - & 201 pts_value(0,9) )**2 ! u"2 202 pts_value_l(0,25) = pts_value_l(0,25) + ( particles(n)%rvar2 - & 203 pts_value(0,10) )**2 ! v"2 204 pts_value_l(0,26) = pts_value_l(0,26) + ( particles(n)%rvar3 - & 205 pts_value(0,11) )**2 ! w"2 212 pts_value_l(0,23) = pts_value_l(0,23) + ( particles(n)%speed_x - & 213 pts_value(0,6) )**2 ! u*2 214 pts_value_l(0,24) = pts_value_l(0,24) + ( particles(n)%speed_y - & 215 pts_value(0,7) )**2 ! v*2 216 pts_value_l(0,25) = pts_value_l(0,25) + ( particles(n)%speed_z - & 217 pts_value(0,8) )**2 ! w*2 218 IF ( .NOT. curvature_solution_effects ) THEN 219 pts_value_l(0,26) = pts_value_l(0,26) + ( particles(n)%rvar1 - & 220 pts_value(0,9) )**2 ! u"2 221 pts_value_l(0,27) = pts_value_l(0,27) + ( particles(n)%rvar2 - & 222 pts_value(0,10) )**2 ! v"2 223 pts_value_l(0,28) = pts_value_l(0,28) + ( particles(n)%rvar3 - & 224 pts_value(0,11) )**2 ! w"2 225 ENDIF 206 226 ! 207 227 !-- Repeat the same for the respective particle group … … 209 229 j = particles(n)%group 210 230 211 pts_value_l(j, 18) = pts_value_l(j,18) + ( particles(n)%x - &231 pts_value_l(j,20) = pts_value_l(j,20) + ( particles(n)%x - & 212 232 particles(n)%origin_x - pts_value(j,2) )**2 213 pts_value_l(j, 19) = pts_value_l(j,19) + ( particles(n)%y - &233 pts_value_l(j,21) = pts_value_l(j,21) + ( particles(n)%y - & 214 234 particles(n)%origin_y - pts_value(j,3) )**2 215 pts_value_l(j,2 0) = pts_value_l(j,20) + ( particles(n)%z - &235 pts_value_l(j,22) = pts_value_l(j,22) + ( particles(n)%z - & 216 236 particles(n)%origin_z - pts_value(j,4) )**2 217 pts_value_l(j,21) = pts_value_l(j,21) + ( particles(n)%speed_x - & 218 pts_value(j,6) )**2 219 pts_value_l(j,22) = pts_value_l(j,22) + ( particles(n)%speed_y - & 220 pts_value(j,7) )**2 221 pts_value_l(j,23) = pts_value_l(j,23) + ( particles(n)%speed_z - & 222 pts_value(j,8) )**2 223 pts_value_l(j,24) = pts_value_l(j,24) + ( particles(n)%rvar1 - & 224 pts_value(j,9) )**2 225 pts_value_l(j,25) = pts_value_l(j,25) + ( particles(n)%rvar2 - & 226 pts_value(j,10) )**2 227 pts_value_l(j,26) = pts_value_l(j,26) + ( particles(n)%rvar3 - & 228 pts_value(j,11) )**2 237 pts_value_l(j,23) = pts_value_l(j,23) + ( particles(n)%speed_x - & 238 pts_value(j,6) )**2 239 pts_value_l(j,24) = pts_value_l(j,24) + ( particles(n)%speed_y - & 240 pts_value(j,7) )**2 241 pts_value_l(j,25) = pts_value_l(j,25) + ( particles(n)%speed_z - & 242 pts_value(j,8) )**2 243 IF ( .NOT. curvature_solution_effects ) THEN 244 pts_value_l(j,26) = pts_value_l(j,26) + ( particles(n)%rvar1 - & 245 pts_value(j,9) )**2 246 pts_value_l(j,27) = pts_value_l(j,27) + ( particles(n)%rvar2 - & 247 pts_value(j,10) )**2 248 pts_value_l(j,28) = pts_value_l(j,28) + ( particles(n)%rvar3 - & 249 pts_value(j,11) )**2 250 ENDIF 229 251 ENDIF 230 252 231 253 ENDDO 232 254 233 pts_value_l(0,2 7) = ( number_of_particles - pts_value(0,1) / numprocs )**2255 pts_value_l(0,29) = ( number_of_particles - pts_value(0,1) / numprocs )**2 234 256 ! variance of particle numbers 235 257 IF ( number_of_particle_groups > 1 ) THEN 236 258 DO j = 1, number_of_particle_groups 237 pts_value_l(j,2 7) = ( pts_value_l(j,1) - &259 pts_value_l(j,29) = ( pts_value_l(j,1) - & 238 260 pts_value(j,1) / numprocs )**2 239 261 ENDDO … … 246 268 247 269 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 248 CALL MPI_ALLREDUCE( pts_value_l(0, 18), pts_value(0,18), inum*10, MPI_REAL, &270 CALL MPI_ALLREDUCE( pts_value_l(0,20), pts_value(0,20), inum*10, MPI_REAL, & 249 271 MPI_SUM, comm2d, ierr ) 250 272 #else 251 pts_value(:, 18:27) = pts_value_l(:,18:27)273 pts_value(:,20:29) = pts_value_l(:,20:29) 252 274 #endif 253 275 … … 264 286 265 287 IF ( pts_value(j,1) > 0.0 ) THEN 266 pts_value(j, 18:26) = pts_value(j,18:26) / pts_value(j,1)267 ENDIF 268 pts_value(j,2 7) = pts_value(j,27) / numprocs288 pts_value(j,20:28) = pts_value(j,20:28) / pts_value(j,1) 289 ENDIF 290 pts_value(j,29) = pts_value(j,29) / numprocs 269 291 270 292 ENDDO … … 286 308 #endif 287 309 310 DEALLOCATE( pts_value, pts_value_l ) 311 288 312 CALL cpu_log( log_point(36), 'data_output_ptseries','stop', 'nobarrier' ) 289 313 -
palm/trunk/SOURCE/diffusion_e.f90
r791 r825 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! wang_collision_kernel renamed wang_kernel 7 7 ! 8 8 ! Former revisions: … … 164 164 !-- Store dissipation if needed for calculating the sgs particle 165 165 !-- velocities 166 IF ( use_sgs_for_particles .OR. wang_ collision_kernel ) THEN166 IF ( use_sgs_for_particles .OR. wang_kernel ) THEN 167 167 DO j = nys, nyn 168 168 DO k = nzb_s_inner(j,i)+1, nzt … … 254 254 !-- Store dissipation if needed for calculating the sgs particle 255 255 !-- velocities 256 IF ( use_sgs_for_particles .OR. wang_ collision_kernel ) THEN256 IF ( use_sgs_for_particles .OR. wang_kernel ) THEN 257 257 DO j = nys, nyn 258 258 DO k = nzb_s_inner(j,i)+1, nzt … … 268 268 ! 269 269 !-- Boundary condition for dissipation 270 IF ( use_sgs_for_particles .OR. wang_ collision_kernel ) THEN270 IF ( use_sgs_for_particles .OR. wang_kernel ) THEN 271 271 DO i = nxl, nxr 272 272 DO j = nys, nyn … … 372 372 ! 373 373 !-- Store dissipation if needed for calculating the sgs particle velocities 374 IF ( use_sgs_for_particles .OR. wang_ collision_kernel ) THEN374 IF ( use_sgs_for_particles .OR. wang_kernel ) THEN 375 375 DO k = nzb_s_inner(j,i)+1, nzt 376 376 diss(k,j,i) = dissipation(k) -
palm/trunk/SOURCE/header.f90
r824 r825 1478 1478 WRITE ( io, 433 ) 1479 1479 IF ( curvature_solution_effects ) WRITE ( io, 434 ) 1480 IF ( collision_kernel /= 'none' ) THEN 1481 WRITE ( io, 435 ) TRIM( collision_kernel ) 1482 ELSE 1483 WRITE ( io, 436 ) 1484 ENDIF 1480 1485 ENDIF 1481 1486 … … 1955 1960 434 FORMAT (' Curvature and solution effecs are considered for growth of', & 1956 1961 ' droplets < 1.0E-6 m') 1962 435 FORMAT (' Droplet collision is handled by ',A,'-kernel') 1963 436 FORMAT (' Droplet collision is switched off') 1957 1964 450 FORMAT (//' LES / Turbulence quantities:'/ & 1958 1965 ' ---------------------------'/) -
palm/trunk/SOURCE/init_3d_model.f90
r791 r825 7 7 ! Current revisions: 8 8 ! ------------------ 9 ! 9 ! wang_collision_kernel renamed wang_kernel 10 10 ! 11 11 ! Former revisions: … … 334 334 !-- 3D-array for storing the dissipation, needed for calculating the sgs 335 335 !-- particle velocities 336 IF ( use_sgs_for_particles .OR. wang_ collision_kernel ) THEN336 IF ( use_sgs_for_particles .OR. wang_kernel ) THEN 337 337 ALLOCATE ( diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 338 338 ELSE -
palm/trunk/SOURCE/lpm_collision_kernels.f90
r824 r825 1 MODULE wang_kernel_mod1 MODULE lpm_collision_kernels_mod 2 2 3 3 !------------------------------------------------------------------------------! 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! routine renamed from wang_kernel to lpm_collision_kernels, 7 ! turbulence_effects on collision replaced by wang_kernel 7 8 ! 8 9 ! Former revisions: … … 106 107 ALLOCATE( ec(pstart:pend,pstart:pend), winf(pstart:pend) ) 107 108 108 IF ( turbulence_effects_on_collision) THEN109 IF ( wang_kernel ) THEN 109 110 110 111 ALLOCATE( gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend) ) … … 756 757 END SUBROUTINE turb_enhance_eff 757 758 758 END MODULE wang_kernel_mod759 END MODULE lpm_collision_kernels_mod -
palm/trunk/SOURCE/modules.f90
r824 r825 6 6 ! +bfactor, curvature_solution_effects, eps_ros, molecular_weight_of_water, 7 7 ! vanthoff, -b_cond in cloud_parameters, 8 ! dopts_num increased to 2 7, particle attributes speed_x|y|z_sgs renamed8 ! dopts_num increased to 29, particle attributes speed_x|y|z_sgs renamed 9 9 ! rvar1|2|3 10 ! wang_collision_kernel and turbulence_effects_on_collision replaced by 11 ! collision_kernel, hall_kernel, palm_kernel, wang_kernel 10 12 ! 11 13 ! Former revisions: … … 1047 1049 #endif 1048 1050 1049 INTEGER, PARAMETER :: dopr_norm_num = 7, dopts_num = 2 7, dots_max = 1001051 INTEGER, PARAMETER :: dopr_norm_num = 7, dopts_num = 29, dots_max = 100 1050 1052 1051 1053 INTEGER :: dots_num = 23 … … 1062 1064 (/ 'tnpt ', 'x_ ', 'y_ ', 'z_ ', 'z_abs ', 'u ', & 1063 1065 'v ', 'w ', 'u" ', 'v" ', 'w" ', 'npt_up ', & 1064 'w_up ', 'w_down ', 'radius ', ' npt_max', 'npt_min', 'x*2', &1065 ' y*2 ', 'z*2 ', 'u*2 ', 'v*2 ', 'w*2 ', 'u"2 ', &1066 ' v"2 ', 'w"2 ', 'npt*2 ' /)1066 'w_up ', 'w_down ', 'radius ', 'r_min ', 'r_max ', 'npt_max', & 1067 'npt_min', 'x*2 ', 'y*2 ', 'z*2 ', 'u*2 ', 'v*2 ', & 1068 'w*2 ', 'u"2 ', 'v"2 ', 'w"2 ', 'npt*2 ' /) 1067 1069 1068 1070 CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_unit = & 1069 1071 (/ 'number ', 'm ', 'm ', 'm ', 'm ', 'm/s ', & 1070 1072 'm/s ', 'm/s ', 'm/s ', 'm/s ', 'm/s ', 'number ', & 1071 'm/s ', 'm/s ', 'm ', ' number ', 'number ', 'm2', &1072 ' m2 ', 'm2 ', 'm2/s2 ', 'm2/s2', 'm2/s2 ', 'm2/s2 ', &1073 'm2/s2 ', 'm2/s2 ', ' number2' /)1073 'm/s ', 'm/s ', 'm ', 'm ', 'm ', 'number ', & 1074 'number ', 'm2 ', 'm2 ', 'm2 ', 'm2/s2 ', 'm2/s2 ', & 1075 'm2/s2 ', 'm2/s2 ', 'm2/s2 ', 'm2/s2 ', 'number2' /) 1074 1076 1075 1077 CHARACTER (LEN=7), DIMENSION(dots_max) :: dots_label = & … … 1184 1186 1185 1187 CHARACTER (LEN=15) :: bc_par_lr = 'cyclic', bc_par_ns = 'cyclic', & 1186 bc_par_b = 'reflect', bc_par_t = 'absorb' 1188 bc_par_b = 'reflect', bc_par_t = 'absorb', & 1189 collision_kernel = 'none' 1187 1190 1188 1191 #if defined( __parallel ) … … 1206 1209 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: prt_count, prt_start_index 1207 1210 1208 LOGICAL :: turbulence_effects_on_collision = .FALSE.,&1211 LOGICAL :: hall_kernel = .FALSE., palm_kernel = .FALSE., & 1209 1212 particle_advection = .FALSE., random_start_position = .FALSE., & 1210 1213 read_particles_from_restartfile = .TRUE., & 1211 1214 uniform_particles = .TRUE., use_particle_tails = .FALSE., & 1212 1215 use_sgs_for_particles = .FALSE., & 1213 wang_collision_kernel = .FALSE., & 1214 write_particle_statistics = .FALSE. 1216 wang_kernel = .FALSE., write_particle_statistics = .FALSE. 1215 1217 1216 1218 LOGICAL, DIMENSION(max_number_of_particle_groups) :: & -
palm/trunk/SOURCE/package_parin.f90
r791 r825 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! wang_collision_kernel and turbulence_effects_on_collision in particles_par 7 ! replaced by collision_kernel 7 8 ! 8 9 ! Former revisions: … … 71 72 72 73 NAMELIST /particles_par/ bc_par_b, bc_par_lr, bc_par_ns, bc_par_t, & 74 collision_kernel, & 73 75 density_ratio, radius, dt_dopts, & 74 76 dt_min_part, dt_prel, dt_sort_particles, & … … 87 89 read_particles_from_restartfile, & 88 90 skip_particles_for_tail, & 89 turbulence_effects_on_collision, &90 91 use_particle_tails, use_sgs_for_particles, & 91 92 vertical_particle_advection, & 92 wang_collision_kernel, &93 93 write_particle_statistics 94 94 NAMELIST /spectra_par/ averaging_interval_sp, comp_spectra_level, & -
palm/trunk/SOURCE/time_integration.f90
r791 r825 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! wang_collision_kernel renamed wang_kernel 7 7 ! 8 8 ! Former revisions: … … 218 218 CALL exchange_horiz( ql_vp, nbgp ) 219 219 ENDIF 220 IF ( wang_ collision_kernel ) CALL exchange_horiz( diss, nbgp )220 IF ( wang_kernel ) CALL exchange_horiz( diss, nbgp ) 221 221 222 222 CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'stop' )
Note: See TracChangeset
for help on using the changeset viewer.
