Changeset 824 for palm/trunk/SOURCE/advec_particles.f90
- Timestamp:
- Feb 17, 2012 9:09:57 AM (12 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/advec_particles.f90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_particles.f90
r800 r824 4 4 ! Current revisions: 5 5 ! ------------------ 6 ! 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. 7 9 ! 8 10 ! Former revisions: … … 137 139 IMPLICIT NONE 138 140 139 INTEGER :: agp, deleted_particles, deleted_tails, i, ie, ii, inc, is, j, & 140 jj, js, k, kk, kw, m, n, nc, nd, nn, num_gp, pse, psi, & 141 INTEGER :: agp, deleted_particles, deleted_tails, i, ie, ii, inc, & 142 internal_timestep_count, is, j, & 143 jj, js, jtry, k, kk, kw, m, n, nc, nd, nn, num_gp, pse, psi, & 141 144 tlength, trlp_count, trlp_count_sum, trlp_count_recv, & 142 145 trlp_count_recv_sum, trlpt_count, trlpt_count_recv, & … … 150 153 INTEGER :: gp_outside_of_building(1:8) 151 154 155 INTEGER, PARAMETER :: maxtry = 40 ! for Rosenbrock method 156 152 157 LOGICAL :: dt_3d_reached, dt_3d_reached_l, prt_position 153 158 154 REAL :: aa, arg, bb, cc, dd, delta_r, delta_v, dens_ratio, de_dt, & 155 de_dt_min, de_dx_int, de_dx_int_l, de_dx_int_u, de_dy_int, & 156 de_dy_int_l, de_dy_int_u, de_dz_int, de_dz_int_l, de_dz_int_u, & 157 diss_int, diss_int_l, diss_int_u, distance, dt_gap, & 158 dt_particle, dt_particle_m, d_radius, d_sum, e_a, e_int, & 159 e_int_l, e_int_u, e_mean_int, e_s, exp_arg, exp_term, fs_int, & 160 gg, integral, lagr_timescale, lw_max, mean_r, new_r, p_int, & 161 pt_int, pt_int_l, pt_int_u, q_int, q_int_l, q_int_u, ql_int, & 162 ql_int_l, ql_int_u, random_gauss, sl_r3, sl_r4, t_int, u_int, & 163 u_int_l, u_int_u,vv_int, v_int, v_int_l, v_int_u, w_int, & 164 w_int_l, w_int_u, x, y 159 REAL :: aa, afactor, arg, bb, cc, dd, ddenom, delta_r, delta_v, & 160 dens_ratio, de_dt, de_dt_min, de_dx_int, de_dx_int_l, & 161 de_dx_int_u, de_dy_int, de_dy_int_l, de_dy_int_u, de_dz_int, & 162 de_dz_int_l, de_dz_int_u, diss_int, diss_int_l, diss_int_u, & 163 distance, drdt, drdt_ini, drdt_m, dt_ros, dt_ros_last, & 164 dt_ros_next, dt_ros_sum, dt_ros_sum_ini, d2rdt2, d2rdtdr, & 165 dt_gap, dt_particle, dt_particle_m, d_sum, e_a, e_int, & 166 e_int_l, e_int_u, e_mean_int, e_s, err_ros, errmax, exp_arg, & 167 exp_term, fs_int, gg, g1, g2, g3, g4, integral, & 168 lagr_timescale, lw_max, mean_r, new_r, p_int, pt_int, & 169 pt_int_l, pt_int_u, q_int, q_int_l, q_int_u, ql_int, ql_int_l, & 170 ql_int_u, random_gauss, r_ros, r_ros_ini, sigma, sl_r3, sl_r4, & 171 t_int, u_int, u_int_l, u_int_u,vv_int, v_int, v_int_l, & 172 v_int_u, w_int, w_int_l, w_int_u, x, y 173 ! 174 !-- Parameters for Rosenbrock method 175 REAL, PARAMETER :: a21 = 2.0, a31 = 48.0/25.0, a32 = 6.0/25.0, & 176 a2x = 1.0, a3x = 3.0/5.0, b1 = 19.0/9.0, b2 = 0.5, & 177 b3 = 25.0/108.0, b4 = 125.0/108.0, c21 = -8.0, & 178 c31 = 372.0/25.0, c32 = 12.0/5.0, & 179 c41 = -112.0/125.0, c42 = -54.0/125.0, & 180 c43 = -2.0/5.0, c1x = 0.5, c2x= -3.0/2.0, & 181 c3x = 121.0/50.0, c4x = 29.0/250.0, & 182 errcon = 0.1296, e1 = 17.0/54.0, e2 = 7.0/36.0, & 183 e3 = 0.0, e4 = 125.0/108.0, gam = 0.5, grow = 1.5, & 184 pgrow = -0.25, pshrnk = -1.0/3.0, shrnk = 0.5 165 185 166 186 REAL, DIMENSION(1:30) :: de_dxi, de_dyi, de_dzi, dissi, d_gp_pl, ei … … 337 357 ! 338 358 !-- Change in radius by condensation/evaporation 339 !-- ATTENTION: this is only an approximation for large radii 340 arg = particles(n)%radius**2 + 2.0 * dt_3d * & 341 ( e_a / e_s - 1.0 ) / & 342 ( ( l_d_rv / t_int - 1.0 ) * l_v * rho_l / t_int / & 343 thermal_conductivity_l + & 344 rho_l * r_v * t_int / diff_coeff_l / e_s ) 345 IF ( arg < 1.0E-14 ) THEN 346 new_r = 1.0E-7 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 347 511 ELSE 348 new_r = SQRT( arg ) 512 ! 513 !-- Approximation for large radii, where curvature and solution 514 !-- effects can be neglected 515 arg = particles(n)%radius**2 + 2.0 * dt_3d * & 516 ( e_a / e_s - 1.0 ) / & 517 ( ( l_d_rv / t_int - 1.0 ) * l_v * rho_l / t_int / & 518 thermal_conductivity_l + & 519 rho_l * r_v * t_int / diff_coeff_l / e_s ) 520 IF ( arg < 1.0E-14 ) THEN 521 new_r = 1.0E-7 522 ELSE 523 new_r = SQRT( arg ) 524 ENDIF 525 349 526 ENDIF 350 527 351 528 delta_r = new_r - particles(n)%radius 352 353 ! NOTE: this is the correct formula (indipendent of radius).354 ! nevertheless, it give wrong results for large timesteps355 ! d_radius = 1.0 / particles(n)%radius356 ! delta_r = d_radius * ( e_a / e_s - 1.0 - 3.3E-7 / t_int * d_radius + &357 ! b_cond * d_radius**3 ) / &358 ! ( ( l_d_rv / t_int - 1.0 ) * l_v * rho_l / t_int / &359 ! thermal_conductivity_l + &360 ! rho_l * r_v * t_int / diff_coeff_l / e_s ) * dt_3d361 362 ! new_r = particles(n)%radius + delta_r363 ! IF ( new_r < 1.0E-7 ) new_r = 1.0E-7364 529 365 530 ! … … 389 554 WRITE( message_string, * ) '#1 k=',k,' j=',j,' i=',i, & 390 555 ' e_s=',e_s, ' e_a=',e_a,' t_int=',t_int, & 391 ' &d _radius=',d_radius,' delta_r=',delta_r,&556 ' &delta_r=',delta_r, & 392 557 ' particle_radius=',particles(n)%radius 393 558 CALL message( 'advec_particles', 'PA0144', 2, 2, -1, 6, 1 ) … … 1965 2130 !-- [-5.0*sigma, 5.0*sigma] in order to prevent the SGS velocities 1966 2131 !-- from becoming unrealistically large. 1967 particles(n)% speed_x_sgs = SQRT( 2.0 * sgs_wfu_part * e_int ) *&1968 ( random_gauss( iran_part, 5.0 ) 2132 particles(n)%rvar1 = SQRT( 2.0 * sgs_wfu_part * e_int ) * & 2133 ( random_gauss( iran_part, 5.0 ) & 1969 2134 - 1.0 ) 1970 particles(n)% speed_y_sgs = SQRT( 2.0 * sgs_wfv_part * e_int ) *&1971 ( random_gauss( iran_part, 5.0 ) 2135 particles(n)%rvar2 = SQRT( 2.0 * sgs_wfv_part * e_int ) * & 2136 ( random_gauss( iran_part, 5.0 ) & 1972 2137 - 1.0 ) 1973 particles(n)% speed_z_sgs = SQRT( 2.0 * sgs_wfw_part * e_int ) *&1974 ( random_gauss( iran_part, 5.0 ) 2138 particles(n)%rvar3 = SQRT( 2.0 * sgs_wfw_part * e_int ) * & 2139 ( random_gauss( iran_part, 5.0 ) & 1975 2140 - 1.0 ) 1976 2141 … … 2003 2168 ENDIF 2004 2169 2005 particles(n)% speed_x_sgs = particles(n)%speed_x_sgs -&2006 fs_int * c_0 * diss_int * particles(n)% speed_x_sgs *&2170 particles(n)%rvar1 = particles(n)%rvar1 - & 2171 fs_int * c_0 * diss_int * particles(n)%rvar1 * & 2007 2172 dt_particle / ( 4.0 * sgs_wfu_part * e_int ) + & 2008 2173 ( 2.0 * sgs_wfu_part * de_dt * & 2009 particles(n)% speed_x_sgs /&2174 particles(n)%rvar1 / & 2010 2175 ( 2.0 * sgs_wfu_part * e_int ) + de_dx_int & 2011 2176 ) * dt_particle / 2.0 + & … … 2014 2179 SQRT( dt_particle ) 2015 2180 2016 particles(n)% speed_y_sgs = particles(n)%speed_y_sgs -&2017 fs_int * c_0 * diss_int * particles(n)% speed_y_sgs *&2181 particles(n)%rvar2 = particles(n)%rvar2 - & 2182 fs_int * c_0 * diss_int * particles(n)%rvar2 * & 2018 2183 dt_particle / ( 4.0 * sgs_wfv_part * e_int ) + & 2019 2184 ( 2.0 * sgs_wfv_part * de_dt * & 2020 particles(n)% speed_y_sgs /&2185 particles(n)%rvar2 / & 2021 2186 ( 2.0 * sgs_wfv_part * e_int ) + de_dy_int & 2022 2187 ) * dt_particle / 2.0 + & … … 2025 2190 SQRT( dt_particle ) 2026 2191 2027 particles(n)% speed_z_sgs = particles(n)%speed_z_sgs- &2028 fs_int * c_0 * diss_int * particles(n)% speed_z_sgs* &2192 particles(n)%rvar3 = particles(n)%rvar3 - & 2193 fs_int * c_0 * diss_int * particles(n)%rvar3 * & 2029 2194 dt_particle / ( 4.0 * sgs_wfw_part * e_int ) + & 2030 2195 ( 2.0 * sgs_wfw_part * de_dt * & 2031 particles(n)% speed_z_sgs/ &2196 particles(n)%rvar3 / & 2032 2197 ( 2.0 * sgs_wfw_part * e_int ) + de_dz_int & 2033 2198 ) * dt_particle / 2.0 + & … … 2038 2203 ENDIF 2039 2204 2040 u_int = u_int + particles(n)% speed_x_sgs2041 v_int = v_int + particles(n)% speed_y_sgs2042 w_int = w_int + particles(n)% speed_z_sgs2205 u_int = u_int + particles(n)%rvar1 2206 v_int = v_int + particles(n)%rvar2 2207 w_int = w_int + particles(n)%rvar3 2043 2208 2044 2209 ! … … 3378 3543 particles(n)%speed_z = -particles(n)%speed_z 3379 3544 IF ( use_sgs_for_particles .AND. & 3380 particles(n)% speed_z_sgs> 0.0 ) THEN3381 particles(n)% speed_z_sgs = -particles(n)%speed_z_sgs3545 particles(n)%rvar3 > 0.0 ) THEN 3546 particles(n)%rvar3 = -particles(n)%rvar3 3382 3547 ENDIF 3383 3548 IF ( use_particle_tails .AND. nn /= 0 ) THEN … … 3403 3568 particles(n)%speed_z = -particles(n)%speed_z 3404 3569 IF ( use_sgs_for_particles .AND. & 3405 particles(n)% speed_z_sgs< 0.0 ) THEN3406 particles(n)% speed_z_sgs = -particles(n)%speed_z_sgs3570 particles(n)%rvar3 < 0.0 ) THEN 3571 particles(n)%rvar3 = -particles(n)%rvar3 3407 3572 ENDIF 3408 3573 IF ( use_particle_tails .AND. nn /= 0 ) THEN … … 4187 4352 4188 4353 particles_temp => part_1 4189 part_1 = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &4190 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &4191 0.0, 0, 0, 0, 0 )4354 ! part_1 = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 4355 ! 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 4356 ! 0.0, 0, 0, 0, 0 ) 4192 4357 4193 4358 CASE ( 1 ) 4194 4359 4195 4360 particles_temp => part_2 4196 part_2 = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &4197 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &4198 0.0, 0, 0, 0, 0 )4361 ! part_2 = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 4362 ! 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 4363 ! 0.0, 0, 0, 0, 0 ) 4199 4364 4200 4365 END SELECT
Note: See TracChangeset
for help on using the changeset viewer.
