Changeset 1071
- Timestamp:
- Nov 29, 2012 4:54:55 PM (12 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 3 edited
-
lpm_collision_kernels.f90 (modified) (3 diffs)
-
lpm_droplet_collision.f90 (modified) (8 diffs)
-
lpm_droplet_condensation.f90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lpm_collision_kernels.f90
r1037 r1071 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Bugfix: collision efficiencies for Hall kernel should not be < 1.0E-20 23 23 ! 24 24 ! Former revisions: … … 607 607 608 608 !------------------------------------------------------------------------------! 609 ! Calculation of collision effic encies for the Hall kernel609 ! Calculation of collision efficiencies for the Hall kernel 610 610 !------------------------------------------------------------------------------! 611 611 SUBROUTINE effic … … 728 728 ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq) 729 729 ec(j,i) = MIN( ek, 1.0 ) 730 ENDIF 730 ENDIF 731 732 IF ( ec(j,i) < 1.0E-20 ) ec(j,i) = 0.0 731 733 732 734 ec(i,j) = ec(j,i) 733 IF ( ec(i,j) < 1.0E-20 ) ec(i,j) = 0.0734 735 735 736 ENDDO -
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 -
palm/trunk/SOURCE/lpm_droplet_condensation.f90
r1037 r1071 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Ventilation effect for evaporation of large droplets included 23 ! Check for unreasonable results included in calculation of Rosenbrock method 24 ! since physically unlikely results were observed and for the same 25 ! reason the first internal time step in Rosenbrock method should be < 1.0E02 in 26 ! case of evaporation 27 ! Unnecessary calculation of ql_int removed 28 ! Unnecessary calculations in Rosenbrock method (d2rdt2, drdt_m, dt_ros_last) 29 ! removed 30 ! Bugfix: factor in calculation of surface tension changed from 0.00155 to 31 ! 0.000155 23 32 ! 24 33 ! Former revisions: … … 55 64 IMPLICIT NONE 56 65 57 INTEGER :: i, internal_timestep_count, j, jtry, k, n 66 INTEGER :: i, internal_timestep_count, j, jtry, k, n, ros_count 58 67 59 68 INTEGER, PARAMETER :: maxtry = 40 60 69 70 LOGICAL :: repeat 71 61 72 REAL :: aa, afactor, arg, bb, cc, dd, ddenom, delta_r, drdt, drdt_ini, & 62 drdt_m, dt_ros, dt_ros_last, dt_ros_next, dt_ros_sum, & 63 dt_ros_sum_ini, d2rdtdr, d2rdt2, errmax, err_ros, g1, g2, g3, g4, & 64 e_a, e_s, gg, new_r, p_int, pt_int, pt_int_l, pt_int_u, q_int, & 65 q_int_l, q_int_u, ql_int, ql_int_l, ql_int_u, r_ros, r_ros_ini, & 66 sigma, t_int, x, y 73 dt_ros, dt_ros_next, dt_ros_sum, dt_ros_sum_ini, d2rdtdr, errmax, & 74 err_ros, g1, g2, g3, g4, e_a, e_s, gg, new_r, p_int, pt_int, & 75 pt_int_l, pt_int_u, q_int, q_int_l, q_int_u, ql_int, ql_int_l, & 76 ql_int_u, r_ros, r_ros_ini, sigma, t_int, x, y, re_p 67 77 68 78 ! 69 79 !-- Parameters for Rosenbrock method 70 80 REAL, PARAMETER :: a21 = 2.0, a31 = 48.0/25.0, a32 = 6.0/25.0, & 71 a2x = 1.0, a3x = 3.0/5.0, b1 = 19.0/9.0, b2 = 0.5, & 72 b3 = 25.0/108.0, b4 = 125.0/108.0, c21 = -8.0, & 73 c31 = 372.0/25.0, c32 = 12.0/5.0, & 74 c41 = -112.0/125.0, c42 = -54.0/125.0, & 75 c43 = -2.0/5.0, c1x = 0.5, c2x= -3.0/2.0, & 76 c3x = 121.0/50.0, c4x = 29.0/250.0, & 81 b1 = 19.0/9.0, b2 = 0.5, b3 = 25.0/108.0, & 82 b4 = 125.0/108.0, c21 = -8.0, c31 = 372.0/25.0, & 83 c32 = 12.0/5.0, c41 = -112.0/125.0, & 84 c42 = -54.0/125.0, c43 = -2.0/5.0, & 77 85 errcon = 0.1296, e1 = 17.0/54.0, e2 = 7.0/36.0, & 78 86 e3 = 0.0, e4 = 125.0/108.0, gam = 0.5, grow = 1.5, & … … 121 129 ( q_int_u - q_int_l ) 122 130 123 ql_int_l = ( ( gg - aa ) * ql(k,j,i) + ( gg - bb ) * ql(k,j,i+1) &124 + ( gg - cc ) * ql(k,j+1,i) + ( gg - dd ) * ql(k,j+1,i+1) &125 ) / ( 3.0 * gg )126 127 ql_int_u = ( ( gg-aa ) * ql(k+1,j,i) + ( gg-bb ) * ql(k+1,j,i+1) &128 + ( gg-cc ) * ql(k+1,j+1,i) + ( gg-dd ) * ql(k+1,j+1,i+1) &129 ) / ( 3.0 * gg )130 131 ql_int = ql_int_l + ( particles(n)%z - zu(k) ) / dz * &132 ( ql_int_u - ql_int_l )133 134 131 ! 135 132 !-- Calculate real temperature and saturation vapor pressure … … 151 148 ! 152 149 !-- Change in radius by condensation/evaporation 153 IF ( particles(n)%radius >= 1.0E-6 .OR. & 154 .NOT. curvature_solution_effects ) THEN 155 ! 156 !-- Approximation for large radii, where curvature and solution 157 !-- effects can be neglected 150 IF ( particles(n)%radius >= 4.0E-5 .AND. e_a/e_s < 1.0 ) THEN 151 ! 152 !-- Approximation for large radii, where curvature and solution effects 153 !-- can be neglected but ventilation effect has to be included in case of 154 !-- evaporation. 155 !-- First calculate the droplet's Reynolds number 156 re_p = 2.0 * particles(n)%radius * ABS( particles(n)%speed_z ) / & 157 molecular_viscosity 158 ! 159 !-- Ventilation coefficient after Rogers and Yau, 1989 160 IF ( re_p > 2.5 ) THEN 161 afactor = 0.78 + 0.28 * SQRT( re_p ) 162 ELSE 163 afactor = 1.0 + 0.09 * re_p 164 ENDIF 165 166 arg = particles(n)%radius**2 + 2.0 * dt_3d * afactor * & 167 ( e_a / e_s - 1.0 ) / & 168 ( ( l_d_rv / t_int - 1.0 ) * l_v * rho_l / t_int / & 169 thermal_conductivity_l + & 170 rho_l * r_v * t_int / diff_coeff_l / e_s ) 171 172 new_r = SQRT( arg ) 173 174 ELSEIF ( particles(n)%radius >= 1.0E-6 .OR. & 175 .NOT. curvature_solution_effects ) THEN 176 ! 177 !-- Approximation for larger radii in case that curvature and solution 178 !-- effects are neglected and ventilation effects does not play a role 158 179 arg = particles(n)%radius**2 + 2.0 * dt_3d * & 159 180 ( e_a / e_s - 1.0 ) / & … … 178 199 !-- For larger radii the simple analytic method (see ELSE) gives 179 200 !-- almost the same results. 180 ! 181 !-- Surface tension after (Straka, 2009) 182 sigma = 0.0761 - 0.00155 * ( t_int - 273.15 ) 183 184 r_ros = particles(n)%radius 185 dt_ros_sum = 0.0 ! internal integrated time (s) 186 internal_timestep_count = 0 187 188 ddenom = 1.0 / ( rho_l * r_v * t_int / ( e_s * diff_coeff_l ) + & 189 ( l_v / ( r_v * t_int ) - 1.0 ) * & 190 rho_l * l_v / ( thermal_conductivity_l * t_int )& 191 ) 192 193 afactor = 2.0 * sigma / ( rho_l * r_v * t_int ) 194 195 IF ( particles(n)%rvar3 == -9999999.9 ) THEN 196 ! 197 !-- First particle timestep. Derivative has to be calculated. 198 drdt_m = ddenom / r_ros * ( e_a / e_s - 1.0 - & 199 afactor / r_ros + & 201 202 ros_count = 0 203 repeat = .TRUE. 204 ! 205 !-- Carry out the Rosenbrock algorithm. In case of unreasonable results 206 !-- the switch "repeat" will be set true and the algorithm will be carried 207 !-- out again with the internal time step set to its initial (small) value. 208 !-- Unreasonable results may occur if the external conditions, especially the 209 !-- supersaturation, has significantly changed compared to the last PALM 210 !-- timestep. 211 DO WHILE ( repeat ) 212 213 repeat = .FALSE. 214 ! 215 !-- Surface tension after (Straka, 2009) 216 sigma = 0.0761 - 0.000155 * ( t_int - 273.15 ) 217 218 r_ros = particles(n)%radius 219 dt_ros_sum = 0.0 ! internal integrated time (s) 220 internal_timestep_count = 0 221 222 ddenom = 1.0 / ( rho_l * r_v * t_int / ( e_s * diff_coeff_l ) + & 223 ( l_v / ( r_v * t_int ) - 1.0 ) * & 224 rho_l * l_v / ( thermal_conductivity_l * t_int )& 225 ) 226 227 afactor = 2.0 * sigma / ( rho_l * r_v * t_int ) 228 229 ! 230 !-- Take internal time step values from the end of last PALM time step 231 dt_ros_next = particles(n)%rvar1 232 233 ! 234 !-- Internal time step should not be > 1.0E-2 in case of evaporation 235 !-- because larger values may lead to secondary solutions which are 236 !-- physically unlikely 237 IF ( dt_ros_next > 1.0E-2 .AND. e_a/e_s < 1.0 ) THEN 238 dt_ros_next = 1.0E-3 239 ENDIF 240 ! 241 !-- If calculation of Rosenbrock method is repeated due to unreasonalble 242 !-- results during previous try the initial internal time step has to be 243 !-- reduced 244 IF ( ros_count > 1 ) THEN 245 dt_ros_next = dt_ros_next - ( 0.2 * dt_ros_next ) 246 ELSEIF ( ros_count > 5 ) THEN 247 ! 248 !-- Prevent creation of infinite loop 249 message_string = 'ros_count > 5 in Rosenbrock method' 250 CALL message( 'lpm_droplet_condensation', 'PA0018', 2, 2, & 251 0, 6, 0 ) 252 ENDIF 253 254 ! 255 !-- Internal time step must not be larger than PALM time step 256 dt_ros = MIN( dt_ros_next, dt_3d ) 257 ! 258 !-- Integrate growth equation in time unless PALM time step is reached 259 DO WHILE ( dt_ros_sum < dt_3d ) 260 261 internal_timestep_count = internal_timestep_count + 1 262 263 ! 264 !-- Derivative at starting value 265 drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - afactor / r_ros + & 200 266 bfactor / r_ros**3 ) 201 ELSE 202 ! 203 !-- Take value from last PALM timestep 204 drdt_m = particles(n)%rvar3 205 ENDIF 206 ! 207 !-- Take internal timestep values from the end of last PALM timestep 208 dt_ros_last = particles(n)%rvar1 209 dt_ros_next = particles(n)%rvar2 210 ! 211 !-- Internal timestep must not be larger than PALM timestep 212 dt_ros = MIN( dt_ros_next, dt_3d ) 213 ! 214 !-- Integrate growth equation in time unless PALM timestep is reached 215 DO WHILE ( dt_ros_sum < dt_3d ) 216 217 internal_timestep_count = internal_timestep_count + 1 218 219 ! 220 !-- Derivative at starting value 221 drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - afactor / r_ros + & 222 bfactor / r_ros**3 ) 223 drdt_ini = drdt 224 dt_ros_sum_ini = dt_ros_sum 225 r_ros_ini = r_ros 226 227 ! 228 !-- Calculate time derivative of dr/dt 229 d2rdt2 = ( drdt - drdt_m ) / dt_ros_last 230 231 ! 232 !-- Calculate radial derivative of dr/dt 233 d2rdtdr = ddenom * ( ( 1.0 - e_a/e_s ) / r_ros**2 + & 234 2.0 * afactor / r_ros**3 - & 235 4.0 * bfactor / r_ros**5 ) 236 ! 237 !-- Adjust stepsize unless required accuracy is reached 238 DO jtry = 1, maxtry+1 239 240 IF ( jtry == maxtry+1 ) THEN 241 message_string = 'maxtry > 40 in Rosenbrock method' 242 CALL message( 'lpm_droplet_condensation', 'PA0347', 2, 2, & 243 0, 6, 0 ) 267 drdt_ini = drdt 268 dt_ros_sum_ini = dt_ros_sum 269 r_ros_ini = r_ros 270 271 ! 272 !-- Calculate radial derivative of dr/dt 273 d2rdtdr = ddenom * ( ( 1.0 - e_a/e_s ) / r_ros**2 + & 274 2.0 * afactor / r_ros**3 - & 275 4.0 * bfactor / r_ros**5 ) 276 ! 277 !-- Adjust stepsize unless required accuracy is reached 278 DO jtry = 1, maxtry+1 279 280 IF ( jtry == maxtry+1 ) THEN 281 message_string = 'maxtry > 40 in Rosenbrock method' 282 CALL message( 'lpm_droplet_condensation', 'PA0347', 2, 2, & 283 0, 6, 0 ) 284 ENDIF 285 286 aa = 1.0 / ( gam * dt_ros ) - d2rdtdr 287 g1 = drdt_ini / aa 288 r_ros = r_ros_ini + a21 * g1 289 drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - & 290 afactor / r_ros + & 291 bfactor / r_ros**3 ) 292 293 g2 = ( drdt + c21 * g1 / dt_ros )& 294 / aa 295 r_ros = r_ros_ini + a31 * g1 + a32 * g2 296 drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - & 297 afactor / r_ros + & 298 bfactor / r_ros**3 ) 299 300 g3 = ( drdt + & 301 ( c31 * g1 + c32 * g2 ) / dt_ros ) / aa 302 g4 = ( drdt + & 303 ( c41 * g1 + c42 * g2 + c43 * g3 ) / dt_ros ) / aa 304 r_ros = r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + b4 * g4 305 306 dt_ros_sum = dt_ros_sum_ini + dt_ros 307 308 IF ( dt_ros_sum == dt_ros_sum_ini ) THEN 309 message_string = 'zero stepsize in Rosenbrock method' 310 CALL message( 'lpm_droplet_condensation', 'PA0348', 2, 2, & 311 0, 6, 0 ) 312 ENDIF 313 ! 314 !-- Calculate error 315 err_ros = e1*g1 + e2*g2 + e3*g3 + e4*g4 316 errmax = 0.0 317 errmax = MAX( errmax, ABS( err_ros / r_ros_ini ) ) / eps_ros 318 ! 319 !-- Leave loop if accuracy is sufficient, otherwise try again 320 !-- with a reduced stepsize 321 IF ( errmax <= 1.0 ) THEN 322 EXIT 323 ELSE 324 dt_ros = SIGN( MAX( ABS( 0.9 * dt_ros * errmax**pshrnk ), & 325 shrnk * ABS( dt_ros ) ), dt_ros ) 326 ENDIF 327 328 ENDDO ! loop for stepsize adjustment 329 330 ! 331 !-- Calculate next internal time step 332 IF ( errmax > errcon ) THEN 333 dt_ros_next = 0.9 * dt_ros * errmax**pgrow 334 ELSE 335 dt_ros_next = grow * dt_ros 244 336 ENDIF 245 337 246 aa = 1.0 / ( gam * dt_ros ) - d2rdtdr 247 g1 = ( drdt_ini + dt_ros * c1x * d2rdt2 ) / aa 248 r_ros = r_ros_ini + a21 * g1 249 drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - & 250 afactor / r_ros + & 251 bfactor / r_ros**3 ) 252 253 g2 = ( drdt + dt_ros * c2x * d2rdt2 + c21 * g1 / dt_ros )& 254 / aa 255 r_ros = r_ros_ini + a31 * g1 + a32 * g2 256 drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - & 257 afactor / r_ros + & 258 bfactor / r_ros**3 ) 259 260 g3 = ( drdt + dt_ros * c3x * d2rdt2 + & 261 ( c31 * g1 + c32 * g2 ) / dt_ros ) / aa 262 g4 = ( drdt + dt_ros * c4x * d2rdt2 + & 263 ( c41 * g1 + c42 * g2 + c43 * g3 ) / dt_ros ) / aa 264 r_ros = r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + b4 * g4 265 266 dt_ros_sum = dt_ros_sum_ini + dt_ros 267 268 IF ( dt_ros_sum == dt_ros_sum_ini ) THEN 269 message_string = 'zero stepsize in Rosenbrock method' 270 CALL message( 'lpm_droplet_condensation', 'PA0348', 2, 2, & 271 0, 6, 0 ) 338 ! 339 !-- Estimated time step is reduced if the PALM time step is exceeded 340 IF ( ( dt_ros_next + dt_ros_sum ) >= dt_3d ) THEN 341 dt_ros = dt_3d - dt_ros_sum 342 ELSE 343 dt_ros = dt_ros_next 272 344 ENDIF 273 ! 274 !-- Calculate error 275 err_ros = e1*g1 + e2*g2 + e3*g3 + e4*g4 276 errmax = 0.0 277 errmax = MAX( errmax, ABS( err_ros / r_ros_ini ) ) / eps_ros 278 ! 279 !-- Leave loop if accuracy is sufficient, otherwise try again 280 !-- with a reduced stepsize 281 IF ( errmax <= 1.0 ) THEN 282 EXIT 283 ELSE 284 dt_ros = SIGN( MAX( ABS( 0.9 * dt_ros * errmax**pshrnk ), & 285 shrnk * ABS( dt_ros ) ), dt_ros ) 286 ENDIF 287 288 ENDDO ! loop for stepsize adjustment 289 290 ! 291 !-- Calculate next internal timestep 292 IF ( errmax > errcon ) THEN 293 dt_ros_next = 0.9 * dt_ros * errmax**pgrow 294 ELSE 295 dt_ros_next = grow * dt_ros 345 346 ENDDO 347 ! 348 !-- Store internal time step value for next PALM step 349 particles(n)%rvar1 = dt_ros_next 350 351 new_r = r_ros 352 ! 353 !-- Radius should not fall below 1E-8 because Rosenbrock method may 354 !-- lead to errors otherwise 355 new_r = MAX( new_r, 1.0E-8 ) 356 ! 357 !-- Check if calculated droplet radius change is reasonable since in 358 !-- case of droplet evaporation the Rosenbrock method may lead to 359 !-- secondary solutions which are physically unlikely. 360 !-- Due to the solution effect the droplets may grow for relative 361 !-- humidities below 100%, but change of radius should not be too large. 362 !-- In case of unreasonable droplet growth the Rosenbrock method is 363 !-- recalculated using a smaller initial time step. 364 !-- Limiting values are tested for droplets down to 1.0E-7 365 IF ( new_r - particles(n)%radius >= 3.0E-7 .AND. & 366 e_a/e_s < 0.97 ) THEN 367 ros_count = ros_count + 1 368 repeat = .TRUE. 296 369 ENDIF 297 370 298 ! 299 !-- Estimated timestep is reduced if the PALM time step is exceeded 300 dt_ros_last = dt_ros 301 IF ( ( dt_ros_next + dt_ros_sum ) >= dt_3d ) THEN 302 dt_ros = dt_3d - dt_ros_sum 303 ELSE 304 dt_ros = dt_ros_next 305 ENDIF 306 307 drdt_m = drdt 308 309 ENDDO 310 ! 311 !-- Store derivative and internal timestep values for next PALM step 312 particles(n)%rvar1 = dt_ros_last 313 particles(n)%rvar2 = dt_ros_next 314 particles(n)%rvar3 = drdt_m 315 316 new_r = r_ros 317 ! 318 !-- Radius should not fall below 1E-8 because Rosenbrock method may 319 !-- lead to errors otherwise 320 new_r = MAX( new_r, 1.0E-8 ) 371 ENDDO ! Rosenbrock method 321 372 322 373 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
