Changeset 1071 for palm/trunk/SOURCE/lpm_droplet_condensation.f90
- Timestamp:
- Nov 29, 2012 4:54:55 PM (11 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/lpm_droplet_condensation.f90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
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.
