Changeset 2312 for palm/trunk/SOURCE/lpm_droplet_condensation.f90
- Timestamp:
- Jul 14, 2017 8:26:51 PM (7 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/lpm_droplet_condensation.f90 (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lpm_droplet_condensation.f90
r2101 r2312 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 23 ! 22 ! 23 ! 24 24 ! Former revisions: 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Rosenbrock scheme improved. Gas-kinetic effect added. 27 28 ! 28 29 ! 2000 2016-08-20 18:09:15Z knoop 29 30 ! Forced header and separation lines into 80 columns 30 ! 31 ! 31 32 ! 1890 2016-04-22 08:52:11Z hoffmann 32 33 ! Some improvements of the Rosenbrock method. If the Rosenbrock method needs more 33 34 ! than 40 iterations to find a sufficient time setp, the model is not aborted. 34 ! This might lead to small erros. But they can be assumend as negligible, since 35 ! the maximum timesetp of the Rosenbrock method after 40 iterations will be 36 ! smaller than 10^-16 s. 37 ! 35 ! This might lead to small erros. But they can be assumend as negligible, since 36 ! the maximum timesetp of the Rosenbrock method after 40 iterations will be 37 ! smaller than 10^-16 s. 38 ! 38 39 ! 1871 2016-04-15 11:46:09Z hoffmann 39 40 ! Initialization of aerosols added. … … 53 54 ! 54 55 ! 1682 2015-10-07 23:56:08Z knoop 55 ! Code annotations made doxygen readable 56 ! 56 ! Code annotations made doxygen readable 57 ! 57 58 ! 1359 2014-04-11 17:15:14Z hoffmann 58 ! New particle structure integrated. 59 ! New particle structure integrated. 59 60 ! Kind definition added to all floating point numbers. 60 61 ! 61 62 ! 1346 2014-03-27 13:18:20Z heinze 62 ! Bugfix: REAL constants provided with KIND-attribute especially in call of 63 ! Bugfix: REAL constants provided with KIND-attribute especially in call of 63 64 ! intrinsic function like MAX, MIN, SIGN 64 65 ! … … 108 109 !------------------------------------------------------------------------------! 109 110 SUBROUTINE lpm_droplet_condensation (ip,jp,kp) 110 111 111 112 112 113 USE arrays_3d, & 113 ONLY: hyp, pt, q, 114 ONLY: hyp, pt, q, ql_c, ql_v 114 115 115 116 USE cloud_parameters, & … … 139 140 use_kernel_tables, vanthoff, wang_kernel 140 141 141 142 142 IMPLICIT NONE 143 143 144 144 INTEGER(iwp) :: ip !< 145 INTEGER(iwp) :: internal_timestep_count !<146 145 INTEGER(iwp) :: jp !< 147 INTEGER(iwp) :: jtry !<148 146 INTEGER(iwp) :: kp !< 149 147 INTEGER(iwp) :: n !< 150 INTEGER(iwp) :: ros_count !< 151 152 INTEGER(iwp), PARAMETER :: maxtry = 40 !< 153 154 LOGICAL :: repeat !< 155 156 REAL(wp) :: aa !< 148 157 149 REAL(wp) :: afactor !< curvature effects 158 150 REAL(wp) :: arg !< … … 161 153 REAL(wp) :: delta_r !< 162 154 REAL(wp) :: diameter !< diameter of cloud droplets 163 REAL(wp) :: diff_coeff _v!< diffusivity for water vapor155 REAL(wp) :: diff_coeff !< diffusivity for water vapor 164 156 REAL(wp) :: drdt !< 165 REAL(wp) :: drdt_ini !<166 157 REAL(wp) :: dt_ros !< 167 REAL(wp) :: dt_ros_next !<168 158 REAL(wp) :: dt_ros_sum !< 169 REAL(wp) :: dt_ros_sum_ini !<170 159 REAL(wp) :: d2rdtdr !< 171 REAL(wp) :: errmax !<172 160 REAL(wp) :: e_a !< current vapor pressure 173 161 REAL(wp) :: e_s !< current saturation vapor pressure 174 REAL(wp) :: err _ros !<175 REAL(wp) :: g1 !<176 REAL(wp) :: g2 !<177 REAL(wp) :: g3 !<178 REAL(wp) :: g4 !<179 REAL(wp) :: r_ros !<180 REAL(wp) :: r _ros_ini !<181 REAL(wp) :: sigma !< 182 REAL(wp) :: thermal_conductivity _v!< thermal conductivity for water162 REAL(wp) :: error !< local truncation error in Rosenbrock 163 REAL(wp) :: k1 !< 164 REAL(wp) :: k2 !< 165 REAL(wp) :: r_err !< First order estimate of Rosenbrock radius 166 REAL(wp) :: r_ros !< Rosenbrock radius 167 REAL(wp) :: r_ros_ini !< initial Rosenbrock radius 168 REAL(wp) :: r0 !< gas-kinetic lengthscale 169 REAL(wp) :: sigma !< surface tension of water 170 REAL(wp) :: thermal_conductivity !< thermal conductivity for water 183 171 REAL(wp) :: t_int !< temperature 184 172 REAL(wp) :: w_s !< terminal velocity of droplets 185 REAL(wp) :: re_p !< 186 187 ! 188 !-- Parameters for Rosenbrock method 189 REAL(wp), PARAMETER :: a21 = 2.0_wp !< 190 REAL(wp), PARAMETER :: a31 = 48.0_wp / 25.0_wp !< 191 REAL(wp), PARAMETER :: a32 = 6.0_wp / 25.0_wp !< 192 REAL(wp), PARAMETER :: b1 = 19.0_wp / 9.0_wp !< 193 REAL(wp), PARAMETER :: b2 = 0.5_wp !< 194 REAL(wp), PARAMETER :: b3 = 25.0_wp / 108.0_wp !< 195 REAL(wp), PARAMETER :: b4 = 125.0_wp / 108.0_wp !< 196 REAL(wp), PARAMETER :: c21 = -8.0_wp !< 197 REAL(wp), PARAMETER :: c31 = 372.0_wp / 25.0_wp !< 198 REAL(wp), PARAMETER :: c32 = 12.0_wp / 5.0_wp !< 199 REAL(wp), PARAMETER :: c41 = -112.0_wp / 125.0_wp !< 200 REAL(wp), PARAMETER :: c42 = -54.0_wp / 125.0_wp !< 201 REAL(wp), PARAMETER :: c43 = -2.0_wp / 5.0_wp !< 202 REAL(wp), PARAMETER :: errcon = 0.1296_wp !< 203 REAL(wp), PARAMETER :: e1 = 17.0_wp / 54.0_wp !< 204 REAL(wp), PARAMETER :: e2 = 7.0_wp / 36.0_wp !< 205 REAL(wp), PARAMETER :: e3 = 0.0_wp !< 206 REAL(wp), PARAMETER :: e4 = 125.0_wp / 108.0_wp !< 207 REAL(wp), PARAMETER :: eps_ros = 1.0E-3_wp !< accuracy of Rosenbrock method 208 REAL(wp), PARAMETER :: gam = 0.5_wp !< 209 REAL(wp), PARAMETER :: grow = 1.5_wp !< 210 REAL(wp), PARAMETER :: pgrow = -0.25_wp !< 211 REAL(wp), PARAMETER :: pshrnk = -1.0_wp /3.0_wp !< 212 REAL(wp), PARAMETER :: shrnk = 0.5_wp !< 213 173 REAL(wp) :: re_p !< particle Reynolds number 174 ! 175 !-- Parameters for Rosenbrock method (see Verwer et al., 1999) 176 REAL(wp), PARAMETER :: prec = 1.0E-3_wp !< precision of Rosenbrock solution 177 REAL(wp), PARAMETER :: q_increase = 1.5_wp !< increase factor in timestep 178 REAL(wp), PARAMETER :: q_decrease = 0.9_wp !< decrease factor in timestep 179 REAL(wp), PARAMETER :: gamma = 0.292893218814_wp !< = 1.0 - 1.0 / SQRT(2.0) 214 180 ! 215 181 !-- Parameters for terminal velocity … … 224 190 REAL(wp), DIMENSION(number_of_particles) :: new_r !< 225 191 226 227 228 192 CALL cpu_log( log_point_s(42), 'lpm_droplet_condens', 'start' ) 229 193 230 194 ! 231 !-- Calculate temperature, saturation vapor pressure and current vapor pressure195 !-- Absolute temperature 232 196 t_int = pt(kp,jp,ip) * ( hyp(kp) / 100000.0_wp )**0.286_wp 233 e_s = 611.0_wp * EXP( l_d_rv * ( 3.6609E-3_wp - 1.0_wp / t_int ) ) 234 e_a = q(kp,jp,ip) * hyp(kp) / ( 0.378_wp * q(kp,jp,ip) + 0.622_wp ) 235 ! 236 !-- Thermal conductivity for water (from Rogers and Yau, Table 7.1), 237 !-- diffusivity for water vapor (after Hall und Pruppacher, 1976) 238 thermal_conductivity_v = 7.94048E-05_wp * t_int + 0.00227011_wp 239 diff_coeff_v = 0.211E-4_wp * ( t_int / 273.15_wp )**1.94_wp * & 240 ( 101325.0_wp / hyp(kp) ) 241 ! 242 !-- Calculate effects of heat conductivity and diffusion of water vapor on the 243 !-- condensation/evaporation process (typically known as 1.0 / (F_k + F_d) ) 244 ddenom = 1.0_wp / ( rho_l * r_v * t_int / ( e_s * diff_coeff_v ) + & 197 ! 198 !-- Saturation vapor pressure (Eq. 10 in Bolton, 1980) 199 e_s = 611.2_wp * EXP( 17.62_wp * ( t_int - 273.15_wp ) / ( t_int - 29.65_wp ) ) 200 ! 201 !-- Current vapor pressure 202 e_a = q(kp,jp,ip) * hyp(kp) / ( q(kp,jp,ip) + 0.622_wp ) 203 ! 204 !-- Thermal conductivity for water (from Rogers and Yau, Table 7.1) 205 thermal_conductivity = 7.94048E-05_wp * t_int + 0.00227011_wp 206 ! 207 !-- Moldecular diffusivity of water vapor in air (Hall und Pruppacher, 1976) 208 diff_coeff = 0.211E-4_wp * ( t_int / 273.15_wp )**1.94_wp * & 209 ( 101325.0_wp / hyp(kp) ) 210 ! 211 !-- Lengthscale for gas-kinetic effects (from Mordy, 1959, p. 23): 212 r0 = diff_coeff / 0.036_wp * SQRT( 2.0_wp * pi / ( r_v * t_int ) ) 213 ! 214 !-- Calculate effects of heat conductivity and diffusion of water vapor on the 215 !-- diffusional growth process (usually known as 1.0 / (F_k + F_d) ) 216 ddenom = 1.0_wp / ( rho_l * r_v * t_int / ( e_s * diff_coeff ) + & 245 217 ( l_v / ( r_v * t_int ) - 1.0_wp ) * rho_l * & 246 l_v / ( thermal_conductivity _v * t_int )&218 l_v / ( thermal_conductivity * t_int ) & 247 219 ) 248 249 220 new_r = 0.0_wp 250 251 221 ! 252 222 !-- Determine ventilation effect on evaporation of large drops … … 264 234 ENDIF 265 235 ! 266 !-- First calculate droplet's Reynolds number236 !-- Calculate droplet's Reynolds number 267 237 re_p = 2.0_wp * particles(n)%radius * w_s / molecular_viscosity 268 238 ! … … 275 245 ELSE 276 246 ! 277 !-- For small droplets or in supersaturated environments, the ventilation 247 !-- For small droplets or in supersaturated environments, the ventilation 278 248 !-- effect does not play a role 279 249 ventilation_effect(n) = 1.0_wp … … 281 251 ENDDO 282 252 283 !284 !-- Use analytic model for condensational growth285 253 IF( .NOT. curvature_solution_effects ) then 254 ! 255 !-- Use analytic model for diffusional growth including gas-kinetic 256 !-- effects (Mordy, 1959) but without the impact of aerosols. 286 257 DO n = 1, number_of_particles 287 arg = particles(n)%radius**2 + 2.0_wp * dt_3d * ddenom *&288 ventilation_effect(n) *&289 ( e_a / e_s - 1.0_wp )290 arg = MAX( arg, 1.0E-16_wp)291 new_r(n) = SQRT( arg ) 258 arg = ( particles(n)%radius + r0 )**2 + 2.0_wp * dt_3d * ddenom * & 259 ventilation_effect(n) * & 260 ( e_a / e_s - 1.0_wp ) 261 arg = MAX( arg, ( 0.01E-6 + r0 )**2 ) 262 new_r(n) = SQRT( arg ) - r0 292 263 ENDDO 293 ENDIF 294 295 ! 296 !-- If selected, use numerical solution of the condensational growth 297 !-- equation (e.g., for studying the activation of aerosols). 298 !-- Curvature and solutions effects are included in growth equation. 299 !-- Change in Radius is calculated with a 4th-order Rosenbrock method 300 !-- for stiff o.d.e's with monitoring local truncation error to adjust 301 !-- stepsize (see Numerical recipes in FORTRAN, 2nd edition, p. 731). 302 DO n = 1, number_of_particles 303 IF ( curvature_solution_effects ) THEN 304 305 ros_count = 0 306 repeat = .TRUE. 307 ! 308 !-- Carry out the Rosenbrock algorithm. In case of unreasonable results 309 !-- the switch "repeat" will be set true and the algorithm will be carried 310 !-- out again with the internal time step set to its initial (small) value. 311 !-- Unreasonable results may occur if the external conditions, especially 312 !-- the supersaturation, has significantly changed compared to the last 313 !-- PALM timestep. 314 DO WHILE ( repeat ) 315 316 repeat = .FALSE. 317 ! 318 !-- Curvature effect (afactor) with surface tension parameterization 319 !-- by Straka (2009) 320 sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp ) 321 afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int ) 322 ! 323 !-- Solute effect (bfactor) 324 bfactor = vanthoff * rho_s * particles(n)%rvar2**3 * & 325 molecular_weight_of_water / & 326 ( rho_l * molecular_weight_of_solute ) 327 328 r_ros = particles(n)%radius 329 dt_ros_sum = 0.0_wp ! internal integrated time (s) 330 internal_timestep_count = 0 331 ! 332 !-- Take internal time step values from the end of last PALM time step 333 dt_ros_next = particles(n)%rvar1 334 335 ! 336 !-- Internal time step should not be > 1.0E-2 and < 1.0E-6 337 IF ( dt_ros_next > 1.0E-2_wp ) THEN 338 dt_ros_next = 1.0E-2_wp 339 ELSEIF ( dt_ros_next < 1.0E-6_wp ) THEN 340 dt_ros_next = 1.0E-6_wp 341 ENDIF 342 343 ! 344 !-- If calculation of Rosenbrock method is repeated due to unreasonalble 345 !-- results during previous try the initial internal time step has to be 346 !-- reduced 347 IF ( ros_count > 1 ) THEN 348 dt_ros_next = dt_ros_next * 0.1_wp 349 ELSEIF ( ros_count > 5 ) THEN 350 ! 351 !-- Prevent creation of infinite loop 352 message_string = 'ros_count > 5 in Rosenbrock method' 353 CALL message( 'lpm_droplet_condensation', 'PA0018', 2, 2, & 354 0, 6, 0 ) 355 ENDIF 356 357 ! 358 !-- Internal time step must not be larger than PALM time step 359 dt_ros = MIN( dt_ros_next, dt_3d ) 360 361 ! 362 !-- Integrate growth equation in time unless PALM time step is reached 363 DO WHILE ( dt_ros_sum < dt_3d ) 364 365 internal_timestep_count = internal_timestep_count + 1 366 367 ! 368 !-- Derivative at starting value 369 drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0_wp - & 264 265 ELSE 266 ! 267 !-- Integrate the diffusional growth including gas-kinetic (Mordy, 1959), 268 !-- as well as curvature and solute effects (e.g., Köhler, 1936). 269 ! 270 !-- Curvature effect (afactor) with surface tension (sigma) by Straka (2009) 271 sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp ) 272 ! 273 !-- Solute effect (afactor) 274 afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int ) 275 276 DO n = 1, number_of_particles 277 ! 278 !-- Solute effect (bfactor) 279 bfactor = vanthoff * rho_s * particles(n)%aux1**3 * & 280 molecular_weight_of_water / ( rho_l * molecular_weight_of_solute ) 281 282 dt_ros = particles(n)%aux2 ! use previously stored Rosenbrock timestep 283 dt_ros_sum = 0.0_wp 284 285 r_ros = particles(n)%radius ! initialize Rosenbrock particle radius 286 r_ros_ini = r_ros 287 ! 288 !-- Integrate growth equation using a 2nd-order Rosenbrock method 289 !-- (see Verwer et al., 1999, Eq. (3.2)). The Rosenbrock method adjusts 290 !-- its with internal timestep to minimize the local truncation error. 291 DO WHILE ( dt_ros_sum < dt_3d ) 292 293 dt_ros = MIN( dt_ros, dt_3d - dt_ros_sum ) 294 295 DO 296 297 drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0 - & 370 298 afactor / r_ros + & 371 299 bfactor / r_ros**3 & 372 ) / r_ros 373 374 drdt_ini = drdt 375 dt_ros_sum_ini = dt_ros_sum 376 r_ros_ini = r_ros 377 378 ! 379 !-- Calculate radial derivative of dr/dt 380 d2rdtdr = ddenom * ventilation_effect(n) * & 381 ( ( 1.0_wp - e_a / e_s ) / r_ros**2 + & 382 2.0_wp * afactor / r_ros**3 - & 383 4.0_wp * bfactor / r_ros**5 & 384 ) 385 ! 386 !-- Adjust stepsize unless required accuracy is reached 387 DO jtry = 1, maxtry+1 388 389 IF ( jtry == maxtry+1 ) THEN 390 message_string = 'maxtry > 40 in Rosenbrock method' 391 CALL message( 'lpm_droplet_condensation', 'PA0347', 0, & 392 1, 0, 6, 0 ) 393 ENDIF 394 395 aa = 1.0_wp / ( gam * dt_ros ) - d2rdtdr 396 g1 = drdt_ini / aa 397 r_ros = r_ros_ini + a21 * g1 398 drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0_wp - & 399 afactor / r_ros + & 400 bfactor / r_ros**3 & 401 ) / r_ros 402 403 g2 = ( drdt + c21 * g1 / dt_ros )& 404 / aa 405 r_ros = r_ros_ini + a31 * g1 + a32 * g2 406 drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0_wp - & 407 afactor / r_ros + & 408 bfactor / r_ros**3 & 409 ) / r_ros 410 411 g3 = ( drdt + & 412 ( c31 * g1 + c32 * g2 ) / dt_ros ) / aa 413 g4 = ( drdt + & 414 ( c41 * g1 + c42 * g2 + c43 * g3 ) / dt_ros ) / aa 415 r_ros = r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + b4 * g4 416 417 dt_ros_sum = dt_ros_sum_ini + dt_ros 418 419 IF ( dt_ros_sum == dt_ros_sum_ini ) THEN 420 message_string = 'zero stepsize in Rosenbrock method' 421 CALL message( 'lpm_droplet_condensation', 'PA0348', 2, & 422 2, 0, 6, 0 ) 423 ENDIF 424 ! 425 !-- Calculate error 426 err_ros = e1 * g1 + e2 * g2 + e3 * g3 + e4 * g4 427 errmax = 0.0_wp 428 errmax = MAX( errmax, ABS( err_ros / r_ros_ini ) ) / eps_ros 429 ! 430 !-- Leave loop if accuracy is sufficient, otherwise try again 431 !-- with a reduced stepsize 432 IF ( errmax <= 1.0_wp ) THEN 433 EXIT 434 ELSE 435 dt_ros = MAX( ABS( 0.9_wp * dt_ros * errmax**pshrnk ), & 436 shrnk * ABS( dt_ros ) ) 437 ENDIF 438 439 ENDDO ! loop for stepsize adjustment 440 441 ! 442 !-- Calculate next internal time step 443 IF ( errmax > errcon ) THEN 444 dt_ros_next = 0.9_wp * dt_ros * errmax**pgrow 300 ) / ( r_ros + r0 ) 301 302 d2rdtdr = -ddenom * ventilation_effect(n) * ( & 303 (e_a / e_s - 1.0) * r_ros**4 - & 304 afactor * r0 * r_ros**2 - & 305 2.0 * afactor * r_ros**3 + & 306 3.0 * bfactor * r0 + & 307 4.0 * bfactor * r_ros & 308 ) & 309 / ( r_ros**4 * ( r_ros + r0 )**2 ) 310 311 k1 = drdt / ( 1.0 - gamma * dt_ros * d2rdtdr ) 312 313 r_ros = MAX(r_ros_ini + k1 * dt_ros, particles(n)%aux1) 314 r_err = r_ros 315 316 drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0 - & 317 afactor / r_ros + & 318 bfactor / r_ros**3 & 319 ) / ( r_ros + r0 ) 320 321 k2 = ( drdt - dt_ros * 2.0 * gamma * d2rdtdr * k1 ) / & 322 ( 1.0 - dt_ros * gamma * d2rdtdr ) 323 324 r_ros = MAX(r_ros_ini + dt_ros * ( 1.5 * k1 + 0.5 * k2), particles(n)%aux1) 325 ! 326 !-- Check error of the solution, and reduce dt_ros if necessary. 327 error = ABS(r_err - r_ros) / r_ros 328 IF ( error .GT. prec ) THEN 329 dt_ros = SQRT( q_decrease * prec / error ) * dt_ros 330 r_ros = r_ros_ini 445 331 ELSE 446 dt_ros_next = grow * dt_ros 332 dt_ros_sum = dt_ros_sum + dt_ros 333 dt_ros = q_increase * dt_ros 334 r_ros_ini = r_ros 335 EXIT 447 336 ENDIF 448 337 449 ! 450 !-- Estimated time step is reduced if the PALM time step is exceeded 451 IF ( ( dt_ros_next + dt_ros_sum ) >= dt_3d ) THEN 452 dt_ros = dt_3d - dt_ros_sum 453 ELSE 454 dt_ros = dt_ros_next 455 ENDIF 456 457 ENDDO 458 ! 459 !-- Store internal time step value for next PALM step 460 particles(n)%rvar1 = dt_ros_next 461 462 ! 463 !-- Radius should not fall below 1E-8 because Rosenbrock method may 464 !-- lead to errors otherwise 465 new_r(n) = MAX( r_ros, particles(n)%rvar2 ) 466 ! 467 !-- Check if calculated droplet radius change is reasonable since in 468 !-- case of droplet evaporation the Rosenbrock method may lead to 469 !-- secondary solutions which are physically unlikely. 470 !-- Due to the solution effect the droplets may grow for relative 471 !-- humidities below 100%, but change of radius should not be too 472 !-- large. In case of unreasonable droplet growth the Rosenbrock 473 !-- method is recalculated using a smaller initial time step. 474 !-- Limiting values are tested for droplets down to 1.0E-7 475 IF ( new_r(n) - particles(n)%radius >= 3.0E-7_wp .AND. & 476 e_a / e_s < 0.97_wp ) THEN 477 ros_count = ros_count + 1 478 repeat = .TRUE. 479 ENDIF 480 481 ENDDO ! Rosenbrock method 482 483 ENDIF 484 485 delta_r = new_r(n) - particles(n)%radius 486 487 ! 488 !-- Sum up the change in volume of liquid water for the respective grid 489 !-- volume (this is needed later in lpm_calc_liquid_water_content for 490 !-- calculating the release of latent heat) 338 END DO 339 340 END DO !Rosenbrock loop 341 ! 342 !-- Store new particle radius 343 new_r(n) = r_ros 344 ! 345 !-- Store internal time step value for next PALM step 346 particles(n)%aux2 = dt_ros 347 348 ENDDO !Particle loop 349 350 ENDIF 351 352 DO n = 1, number_of_particles 353 ! 354 !-- Sum up the change in liquid water for the respective grid 355 !-- box for the computation of the release/depletion of water vapor 356 !-- and heat. 491 357 ql_c(kp,jp,ip) = ql_c(kp,jp,ip) + particles(n)%weight_factor * & 492 358 rho_l * 1.33333333_wp * pi * & 493 359 ( new_r(n)**3 - particles(n)%radius**3 ) / & 494 360 ( rho_surface * dx * dy * dz ) 361 ! 362 !-- Check if the increase in liqid water is not too big. If this is the case, 363 !-- the model timestep might be too long. 495 364 IF ( ql_c(kp,jp,ip) > 100.0_wp ) THEN 496 365 WRITE( message_string, * ) 'k=',kp,' j=',jp,' i=',ip, & … … 499 368 CALL message( 'lpm_droplet_condensation', 'PA0143', 2, 2, -1, 6, 1 ) 500 369 ENDIF 501 502 ! 503 !-- Change the droplet radius504 IF ( ( new_r(n) - particles(n)%radius ) < 0.0_wp .AND. &505 370 ! 371 !-- Check if the change in the droplet radius is not too big. If this is the 372 !-- case, the model timestep might be too long. 373 delta_r = new_r(n) - particles(n)%radius 374 IF ( delta_r < 0.0_wp .AND. new_r(n) < 0.0_wp ) THEN 506 375 WRITE( message_string, * ) '#1 k=',kp,' j=',jp,' i=',ip, & 507 376 ' e_s=',e_s, ' e_a=',e_a,' t_int=',t_int, & … … 510 379 CALL message( 'lpm_droplet_condensation', 'PA0144', 2, 2, -1, 6, 1 ) 511 380 ENDIF 512 513 381 ! 514 382 !-- Sum up the total volume of liquid water (needed below for 515 383 !-- re-calculating the weighting factors) 516 384 ql_v(kp,jp,ip) = ql_v(kp,jp,ip) + particles(n)%weight_factor * new_r(n)**3 517 518 particles(n)%radius = new_r(n)519 520 385 ! 521 386 !-- Determine radius class of the particle needed for collision 522 IF ( ( hall_kernel .OR. wang_kernel ) .AND. use_kernel_tables ) & 523 THEN 387 IF ( use_kernel_tables ) THEN 524 388 particles(n)%class = ( LOG( new_r(n) ) - rclass_lbound ) / & 525 389 ( rclass_ubound - rclass_lbound ) * & … … 528 392 particles(n)%class = MAX( particles(n)%class, 1 ) 529 393 ENDIF 394 ! 395 !-- Store new radius to particle features 396 particles(n)%radius = new_r(n) 530 397 531 398 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
