Changeset 1359 for palm/trunk/SOURCE/lpm_droplet_condensation.f90
- Timestamp:
- Apr 11, 2014 5:15:14 PM (10 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/lpm_droplet_condensation.f90 (modified) (22 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lpm_droplet_condensation.f90
r1347 r1359 1 SUBROUTINE lpm_droplet_condensation 1 SUBROUTINE lpm_droplet_condensation (ip,jp,kp) 2 2 3 3 !--------------------------------------------------------------------------------! … … 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 23 ! 22 ! New particle structure integrated. 23 ! Kind definition added to all floating point numbers. 24 ! 24 25 ! Former revisions: 25 26 ! ----------------- … … 100 101 101 102 USE particle_attributes, & 102 ONLY: hall_kernel, number_of_particles, offset_ocean_nzt,&103 offset_ocean_nzt _m1, particles, radius_classes,&104 use_kernel_tables, wang_kernel103 ONLY: block_offset, grid_particles, hall_kernel, number_of_particles, & 104 offset_ocean_nzt, offset_ocean_nzt_m1, particles, & 105 radius_classes, use_kernel_tables, wang_kernel 105 106 106 107 … … 108 109 109 110 INTEGER(iwp) :: i !: 111 INTEGER(iwp) :: ip !: 110 112 INTEGER(iwp) :: internal_timestep_count !: 111 113 INTEGER(iwp) :: j !: 114 INTEGER(iwp) :: jp !: 112 115 INTEGER(iwp) :: jtry !: 113 116 INTEGER(iwp) :: k !: 117 INTEGER(iwp) :: kp !: 114 118 INTEGER(iwp) :: n !: 119 INTEGER(iwp) :: nb !: 115 120 INTEGER(iwp) :: ros_count !: 116 121 117 INTEGER(iwp), PARAMETER :: maxtry = 40 !: 118 119 LOGICAL :: repeat !: 122 INTEGER(iwp), PARAMETER :: maxtry = 40 !: 123 124 INTEGER(iwp), DIMENSION(0:7) :: end_index !: 125 INTEGER(iwp), DIMENSION(0:7) :: start_index !: 126 127 LOGICAL :: repeat !: 128 129 LOGICAL, DIMENSION(number_of_particles) :: flag_1 !: 120 130 121 131 REAL(wp) :: aa !: … … 140 150 REAL(wp) :: g3 !: 141 151 REAL(wp) :: g4 !: 142 REAL(wp) :: e_a !:143 REAL(wp) :: e_s !:144 152 REAL(wp) :: gg !: 145 REAL(wp) :: new_r !:146 REAL(wp) :: p_int !:147 153 REAL(wp) :: pt_int !: 148 154 REAL(wp) :: pt_int_l !: … … 154 160 REAL(wp) :: r_ros_ini !: 155 161 REAL(wp) :: sigma !: 156 REAL(wp) :: t_int !:157 162 REAL(wp) :: x !: 158 163 REAL(wp) :: y !: … … 160 165 161 166 !-- Parameters for Rosenbrock method 162 REAL(wp), PARAMETER :: a21 = 2.0 !: 163 REAL(wp), PARAMETER :: a31 = 48.0/25.0_wp !: 164 REAL(wp), PARAMETER :: a32 = 6.0/25.0_wp !: 165 REAL(wp), PARAMETER :: b1 = 19.0/9.0_wp !: 166 REAL(wp), PARAMETER :: b2 = 0.5 !: 167 REAL(wp), PARAMETER :: b3 = 25.0/108.0_wp !: 168 REAL(wp), PARAMETER :: b4 = 125.0/108.0_wp !: 169 REAL(wp), PARAMETER :: c21 = -8.0 !: 170 REAL(wp), PARAMETER :: c31 = 372.0/25.0_wp !: 171 REAL(wp), PARAMETER :: c32 = 12.0/5.0_wp !: 172 REAL(wp), PARAMETER :: c41 = -112.0/125.0_wp !: 173 REAL(wp), PARAMETER :: c42 = -54.0/125.0_wp !: 174 REAL(wp), PARAMETER :: c43 = -2.0/5.0_wp !: 175 REAL(wp), PARAMETER :: errcon = 0.1296 !: 176 REAL(wp), PARAMETER :: e1 = 17.0/54.0_wp !: 177 REAL(wp), PARAMETER :: e2 = 7.0/36.0_wp !: 178 REAL(wp), PARAMETER :: e3 = 0.0 !: 179 REAL(wp), PARAMETER :: e4 = 125.0/108.0_wp !: 180 REAL(wp), PARAMETER :: gam = 0.5 !: 181 REAL(wp), PARAMETER :: grow = 1.5 !: 182 REAL(wp), PARAMETER :: pgrow = -0.25 !: 183 REAL(wp), PARAMETER :: pshrnk = -1.0/3.0_wp !: 184 REAL(wp), PARAMETER :: shrnk = 0.5 !: 185 167 REAL(wp), PARAMETER :: a21 = 2.0_wp !: 168 REAL(wp), PARAMETER :: a31 = 48.0_wp / 25.0_wp !: 169 REAL(wp), PARAMETER :: a32 = 6.0_wp / 25.0_wp !: 170 REAL(wp), PARAMETER :: b1 = 19.0_wp / 9.0_wp !: 171 REAL(wp), PARAMETER :: b2 = 0.5_wp !: 172 REAL(wp), PARAMETER :: b3 = 25.0_wp / 108.0_wp !: 173 REAL(wp), PARAMETER :: b4 = 125.0_wp / 108.0_wp !: 174 REAL(wp), PARAMETER :: c21 = -8.0_wp !: 175 REAL(wp), PARAMETER :: c31 = 372.0_wp / 25.0_wp !: 176 REAL(wp), PARAMETER :: c32 = 12.0_wp / 5.0_wp !: 177 REAL(wp), PARAMETER :: c41 = -112.0_wp / 125.0_wp !: 178 REAL(wp), PARAMETER :: c42 = -54.0_wp / 125.0_wp !: 179 REAL(wp), PARAMETER :: c43 = -2.0_wp / 5.0_wp !: 180 REAL(wp), PARAMETER :: errcon = 0.1296_wp !: 181 REAL(wp), PARAMETER :: e1 = 17.0_wp / 54.0_wp !: 182 REAL(wp), PARAMETER :: e2 = 7.0_wp / 36.0_wp !: 183 REAL(wp), PARAMETER :: e3 = 0.0_wp !: 184 REAL(wp), PARAMETER :: e4 = 125.0_wp / 108.0_wp !: 185 REAL(wp), PARAMETER :: gam = 0.5_wp !: 186 REAL(wp), PARAMETER :: grow = 1.5_wp !: 187 REAL(wp), PARAMETER :: pgrow = -0.25_wp !: 188 REAL(wp), PARAMETER :: pshrnk = -1.0_wp /3.0_wp !: 189 REAL(wp), PARAMETER :: shrnk = 0.5_wp !: 190 191 REAL(wp), DIMENSION(number_of_particles) :: afactor_v !: 192 REAL(wp), DIMENSION(number_of_particles) :: diff_coeff_v !: 193 REAL(wp), DIMENSION(number_of_particles) :: e_s !: 194 REAL(wp), DIMENSION(number_of_particles) :: e_a !: 195 REAL(wp), DIMENSION(number_of_particles) :: new_r !: 196 REAL(wp), DIMENSION(number_of_particles) :: p_int !: 197 REAL(wp), DIMENSION(number_of_particles) :: thermal_conductivity_v !: 198 REAL(wp), DIMENSION(number_of_particles) :: t_int !: 199 REAL(wp), DIMENSION(number_of_particles) :: xv !: 200 REAL(wp), DIMENSION(number_of_particles) :: yv !: 201 REAL(wp), DIMENSION(number_of_particles) :: zv !: 186 202 187 203 188 204 CALL cpu_log( log_point_s(42), 'lpm_droplet_condens', 'start' ) 189 205 206 start_index = grid_particles(kp,jp,ip)%start_index 207 end_index = grid_particles(kp,jp,ip)%end_index 208 209 xv = particles(1:number_of_particles)%x 210 yv = particles(1:number_of_particles)%y 211 zv = particles(1:number_of_particles)%z 212 213 DO nb = 0,7 214 215 i = ip + block_offset(nb)%i_off 216 j = jp + block_offset(nb)%j_off 217 k = kp + block_offset(nb)%k_off 218 219 DO n = start_index(nb), end_index(nb) 220 ! 221 !-- Interpolate temperature and humidity. 222 x = xv(n) - i * dx 223 y = yv(n) - j * dy 224 aa = x**2 + y**2 225 bb = ( dx - x )**2 + y**2 226 cc = x**2 + ( dy - y )**2 227 dd = ( dx - x )**2 + ( dy - y )**2 228 gg = aa + bb + cc + dd 229 230 pt_int_l = ( ( gg - aa ) * pt(k,j,i) + ( gg - bb ) * pt(k,j,i+1) & 231 + ( gg - cc ) * pt(k,j+1,i) + ( gg - dd ) * pt(k,j+1,i+1) & 232 ) / ( 3.0_wp * gg ) 233 234 pt_int_u = ( ( gg-aa ) * pt(k+1,j,i) + ( gg-bb ) * pt(k+1,j,i+1) & 235 + ( gg-cc ) * pt(k+1,j+1,i) + ( gg-dd ) * pt(k+1,j+1,i+1) & 236 ) / ( 3.0_wp * gg ) 237 238 pt_int = pt_int_l + ( particles(n)%z - zu(k) ) / dz * & 239 ( pt_int_u - pt_int_l ) 240 241 q_int_l = ( ( gg - aa ) * q(k,j,i) + ( gg - bb ) * q(k,j,i+1) & 242 + ( gg - cc ) * q(k,j+1,i) + ( gg - dd ) * q(k,j+1,i+1) & 243 ) / ( 3.0_wp * gg ) 244 245 q_int_u = ( ( gg-aa ) * q(k+1,j,i) + ( gg-bb ) * q(k+1,j,i+1) & 246 + ( gg-cc ) * q(k+1,j+1,i) + ( gg-dd ) * q(k+1,j+1,i+1) & 247 ) / ( 3.0_wp * gg ) 248 249 q_int = q_int_l + ( zv(n) - zu(k) ) / dz * & 250 ( q_int_u - q_int_l ) 251 252 ! 253 !-- Calculate real temperature and saturation vapor pressure 254 p_int(n) = hyp(k) + ( particles(n)%z - zu(k) ) / dz * & 255 ( hyp(k+1)-hyp(k) ) 256 t_int(n) = pt_int * ( p_int(n) / 100000.0_wp )**0.286_wp 257 258 e_s(n) = 611.0_wp * EXP( l_d_rv * ( 3.6609E-3_wp - 1.0_wp / & 259 t_int(n) ) ) 260 261 ! 262 !-- Current vapor pressure 263 e_a(n) = q_int * p_int(n) / ( 0.378_wp * q_int + 0.622_wp ) 264 265 ENDDO 266 ENDDO 267 268 new_r = 0.0_wp 269 flag_1 = .false. 270 190 271 DO n = 1, number_of_particles 191 272 ! 192 !-- Interpolate temperature and humidity. 193 !-- First determine left, south, and bottom index of the arrays. 194 i = particles(n)%x * ddx 195 j = particles(n)%y * ddy 196 k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz & 197 + offset_ocean_nzt ! only exact if equidistant 198 199 x = particles(n)%x - i * dx 200 y = particles(n)%y - j * dy 201 aa = x**2 + y**2 202 bb = ( dx - x )**2 + y**2 203 cc = x**2 + ( dy - y )**2 204 dd = ( dx - x )**2 + ( dy - y )**2 205 gg = aa + bb + cc + dd 206 207 pt_int_l = ( ( gg - aa ) * pt(k,j,i) + ( gg - bb ) * pt(k,j,i+1) & 208 + ( gg - cc ) * pt(k,j+1,i) + ( gg - dd ) * pt(k,j+1,i+1) & 209 ) / ( 3.0 * gg ) 210 211 pt_int_u = ( ( gg-aa ) * pt(k+1,j,i) + ( gg-bb ) * pt(k+1,j,i+1) & 212 + ( gg-cc ) * pt(k+1,j+1,i) + ( gg-dd ) * pt(k+1,j+1,i+1) & 213 ) / ( 3.0 * gg ) 214 215 pt_int = pt_int_l + ( particles(n)%z - zu(k) ) / dz * & 216 ( pt_int_u - pt_int_l ) 217 218 q_int_l = ( ( gg - aa ) * q(k,j,i) + ( gg - bb ) * q(k,j,i+1) & 219 + ( gg - cc ) * q(k,j+1,i) + ( gg - dd ) * q(k,j+1,i+1) & 220 ) / ( 3.0 * gg ) 221 222 q_int_u = ( ( gg-aa ) * q(k+1,j,i) + ( gg-bb ) * q(k+1,j,i+1) & 223 + ( gg-cc ) * q(k+1,j+1,i) + ( gg-dd ) * q(k+1,j+1,i+1) & 224 ) / ( 3.0 * gg ) 225 226 q_int = q_int_l + ( particles(n)%z - zu(k) ) / dz * & 227 ( q_int_u - q_int_l ) 228 229 ! 230 !-- Calculate real temperature and saturation vapor pressure 231 p_int = hyp(k) + ( particles(n)%z - zu(k) ) / dz * ( hyp(k+1)-hyp(k) ) 232 t_int = pt_int * ( p_int / 100000.0 )**0.286 233 234 e_s = 611.0 * EXP( l_d_rv * ( 3.6609E-3 - 1.0 / t_int ) ) 235 236 ! 237 !-- Current vapor pressure 238 e_a = q_int * p_int / ( 0.378 * q_int + 0.622 ) 239 273 !-- Change in radius by condensation/evaporation 274 IF ( particles(n)%radius >= 4.0E-5_wp .AND. & 275 e_a(n)/e_s(n) < 1.0_wp ) THEN 276 ! 277 !-- Approximation for large radii, where curvature and solution effects 278 !-- can be neglected but ventilation effect has to be included in case of 279 !-- evaporation. 280 !-- First calculate the droplet's Reynolds number 281 re_p = 2.0_wp * particles(n)%radius * ABS( particles(n)%speed_z ) / & 282 molecular_viscosity 283 ! 284 !-- Ventilation coefficient (Rogers and Yau, 1989): 285 IF ( re_p > 2.5_wp ) THEN 286 afactor_v(n) = 0.78_wp + 0.28_wp * SQRT( re_p ) 287 ELSE 288 afactor_v(n) = 1.0_wp + 0.09_wp * re_p 289 ENDIF 290 flag_1(n) = .TRUE. 291 ELSEIF ( particles(n)%radius >= 1.0E-6_wp .OR. & 292 .NOT. curvature_solution_effects ) THEN 293 ! 294 !-- Approximation for larger radii in case that curvature and solution 295 !-- effects are neglected and ventilation effects does not play a role 296 afactor_v(n) = 1.0_wp 297 flag_1(n) = .TRUE. 298 ENDIF 299 ENDDO 300 301 DO n = 1, number_of_particles 240 302 ! 241 303 !-- Thermal conductivity for water (from Rogers and Yau, Table 7.1), 242 304 !-- diffusivity for water vapor (after Hall und Pruppacher, 1976) 243 thermal_conductivity_l = 7.94048E-05 * t_int + 0.00227011 244 diff_coeff_l = 0.211E-4 * ( t_int / 273.15 )**1.94 * & 245 ( 101325.0 / p_int) 246 ! 247 !-- Change in radius by condensation/evaporation 248 IF ( particles(n)%radius >= 4.0E-5 .AND. e_a/e_s < 1.0 ) THEN 249 ! 250 !-- Approximation for large radii, where curvature and solution effects 251 !-- can be neglected but ventilation effect has to be included in case of 252 !-- evaporation. 253 !-- First calculate the droplet's Reynolds number 254 re_p = 2.0 * particles(n)%radius * ABS( particles(n)%speed_z ) / & 255 molecular_viscosity 256 ! 257 !-- Ventilation coefficient after Rogers and Yau, 1989 258 IF ( re_p > 2.5 ) THEN 259 afactor = 0.78 + 0.28 * SQRT( re_p ) 260 ELSE 261 afactor = 1.0 + 0.09 * re_p 262 ENDIF 263 264 arg = particles(n)%radius**2 + 2.0 * dt_3d * afactor * & 265 ( e_a / e_s - 1.0 ) / & 266 ( ( l_d_rv / t_int - 1.0 ) * l_v * rho_l / t_int / & 267 thermal_conductivity_l + & 268 rho_l * r_v * t_int / diff_coeff_l / e_s ) 269 270 new_r = SQRT( arg ) 271 272 ELSEIF ( particles(n)%radius >= 1.0E-6 .OR. & 273 .NOT. curvature_solution_effects ) THEN 274 ! 275 !-- Approximation for larger radii in case that curvature and solution 276 !-- effects are neglected and ventilation effects does not play a role 277 arg = particles(n)%radius**2 + 2.0 * dt_3d * & 278 ( e_a / e_s - 1.0 ) / & 279 ( ( l_d_rv / t_int - 1.0 ) * l_v * rho_l / t_int / & 280 thermal_conductivity_l + & 281 rho_l * r_v * t_int / diff_coeff_l / e_s ) 282 IF ( arg < 1.0E-16 ) THEN 283 new_r = 1.0E-8 284 ELSE 285 new_r = SQRT( arg ) 286 ENDIF 287 ENDIF 288 289 IF ( curvature_solution_effects .AND. & 290 ( ( particles(n)%radius < 1.0E-6 ) .OR. ( new_r < 1.0E-6 ) ) ) & 291 THEN 305 thermal_conductivity_v(n) = 7.94048E-05_wp * t_int(n) + 0.00227011_wp 306 diff_coeff_v(n) = 0.211E-4_wp * & 307 ( t_int(n) / 273.15_wp )**1.94_wp * ( 101325.0_wp / p_int(n)) 308 309 IF(flag_1(n)) then 310 arg = particles(n)%radius**2 + 2.0_wp * dt_3d * afactor_v(n) * & 311 ( e_a(n) / e_s(n) - 1.0_wp ) / & 312 ( ( l_d_rv / t_int(n) - 1.0_wp ) * l_v * rho_l / t_int(n) / & 313 thermal_conductivity_v(n) + & 314 rho_l * r_v * t_int(n) / diff_coeff_v(n) / e_s(n) ) 315 316 arg = MAX( arg, 1.0E-16_wp ) 317 new_r(n) = SQRT( arg ) 318 ENDIF 319 ENDDO 320 321 DO n = 1, number_of_particles 322 IF ( curvature_solution_effects .AND. & 323 ( ( particles(n)%radius < 1.0E-6_wp ) .OR. & 324 ( new_r(n) < 1.0E-6_wp ) ) ) THEN 292 325 ! 293 326 !-- Curvature and solutions effects are included in growth equation. … … 304 337 !-- the switch "repeat" will be set true and the algorithm will be carried 305 338 !-- out again with the internal time step set to its initial (small) value. 306 !-- Unreasonable results may occur if the external conditions, especially the307 !-- supersaturation, has significantly changed compared to the last PALM308 !-- timestep.339 !-- Unreasonable results may occur if the external conditions, especially 340 !-- the supersaturation, has significantly changed compared to the last 341 !-- PALM timestep. 309 342 DO WHILE ( repeat ) 310 343 311 344 repeat = .FALSE. 312 345 ! 313 !-- Surface tension after (Straka, 2009)314 sigma = 0.0761 - 0.000155 * ( t_int - 273.15)346 !-- Surface tension (Straka, 2009): 347 sigma = 0.0761_wp - 0.000155_wp * ( t_int(n) - 273.15_wp ) 315 348 316 349 r_ros = particles(n)%radius 317 dt_ros_sum = 0.0 ! internal integrated time (s)350 dt_ros_sum = 0.0_wp ! internal integrated time (s) 318 351 internal_timestep_count = 0 319 352 320 ddenom = 1.0 / ( rho_l * r_v * t_int / ( e_s * diff_coeff_l ) + & 321 ( l_v / ( r_v * t_int ) - 1.0 ) * & 322 rho_l * l_v / ( thermal_conductivity_l * t_int )& 323 ) 324 325 afactor = 2.0 * sigma / ( rho_l * r_v * t_int ) 353 ddenom = 1.0_wp / ( rho_l * r_v * t_int(n) / ( e_s(n) * & 354 diff_coeff_v(n) ) + ( l_v / & 355 ( r_v * t_int(n) ) - 1.0_wp ) * & 356 rho_l * l_v / ( thermal_conductivity_v(n) * & 357 t_int(n) ) & 358 ) 359 360 afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int(n) ) 326 361 327 362 ! … … 333 368 !-- because larger values may lead to secondary solutions which are 334 369 !-- physically unlikely 335 IF ( dt_ros_next > 1.0E-2 .AND. e_a/e_s < 1.0) THEN336 dt_ros_next = 1.0E-3 370 IF ( dt_ros_next > 1.0E-2_wp .AND. e_a(n)/e_s(n) < 1.0_wp ) THEN 371 dt_ros_next = 1.0E-3_wp 337 372 ENDIF 338 373 ! … … 341 376 !-- reduced 342 377 IF ( ros_count > 1 ) THEN 343 dt_ros_next = dt_ros_next - ( 0.2 * dt_ros_next )378 dt_ros_next = dt_ros_next - ( 0.2_wp * dt_ros_next ) 344 379 ELSEIF ( ros_count > 5 ) THEN 345 380 ! … … 361 396 ! 362 397 !-- Derivative at starting value 363 drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - afactor / r_ros +&364 bfactor / r_ros**3 )398 drdt = ddenom / r_ros * ( e_a(n) / e_s(n) - 1.0_wp - afactor / & 399 r_ros + bfactor / r_ros**3 ) 365 400 drdt_ini = drdt 366 401 dt_ros_sum_ini = dt_ros_sum … … 369 404 ! 370 405 !-- Calculate radial derivative of dr/dt 371 d2rdtdr = ddenom * ( ( 1.0 - e_a/e_s ) / r_ros**2 +&372 2.0 * afactor / r_ros**3 -&373 4.0 * bfactor / r_ros**5 )406 d2rdtdr = ddenom * ( ( 1.0_wp - e_a(n)/e_s(n) ) / r_ros**2 + & 407 2.0_wp * afactor / r_ros**3 - & 408 4.0_wp * bfactor / r_ros**5 ) 374 409 ! 375 410 !-- Adjust stepsize unless required accuracy is reached … … 378 413 IF ( jtry == maxtry+1 ) THEN 379 414 message_string = 'maxtry > 40 in Rosenbrock method' 380 CALL message( 'lpm_droplet_condensation', 'PA0347', 2, 2,&381 0, 6, 0 )415 CALL message( 'lpm_droplet_condensation', 'PA0347', 2, & 416 2, 0, 6, 0 ) 382 417 ENDIF 383 418 384 aa = 1.0 / ( gam * dt_ros ) - d2rdtdr419 aa = 1.0_wp / ( gam * dt_ros ) - d2rdtdr 385 420 g1 = drdt_ini / aa 386 421 r_ros = r_ros_ini + a21 * g1 387 drdt = ddenom / r_ros * ( e_a / e_s - 1.0 -&388 afactor / r_ros + &422 drdt = ddenom / r_ros * ( e_a(n) / e_s(n) - 1.0_wp - & 423 afactor / r_ros + & 389 424 bfactor / r_ros**3 ) 390 425 … … 392 427 / aa 393 428 r_ros = r_ros_ini + a31 * g1 + a32 * g2 394 drdt = ddenom / r_ros * ( e_a / e_s - 1.0 -&395 afactor / r_ros + &429 drdt = ddenom / r_ros * ( e_a(n) / e_s(n) - 1.0_wp - & 430 afactor / r_ros + & 396 431 bfactor / r_ros**3 ) 397 432 … … 406 441 IF ( dt_ros_sum == dt_ros_sum_ini ) THEN 407 442 message_string = 'zero stepsize in Rosenbrock method' 408 CALL message( 'lpm_droplet_condensation', 'PA0348', 2, 2,&409 0, 6, 0 )443 CALL message( 'lpm_droplet_condensation', 'PA0348', 2, & 444 2, 0, 6, 0 ) 410 445 ENDIF 411 446 ! 412 447 !-- Calculate error 413 err_ros = e1 *g1 + e2*g2 + e3*g3 + e4*g4414 errmax = 0.0 448 err_ros = e1 * g1 + e2 * g2 + e3 * g3 + e4 * g4 449 errmax = 0.0_wp 415 450 errmax = MAX( errmax, ABS( err_ros / r_ros_ini ) ) / eps_ros 416 451 ! 417 452 !-- Leave loop if accuracy is sufficient, otherwise try again 418 453 !-- with a reduced stepsize 419 IF ( errmax <= 1.0 ) THEN454 IF ( errmax <= 1.0_wp ) THEN 420 455 EXIT 421 456 ELSE 422 dt_ros = SIGN( MAX( ABS( 0.9 * dt_ros * errmax**pshrnk ), & 423 shrnk * ABS( dt_ros ) ), dt_ros ) 457 dt_ros = SIGN( MAX( ABS( 0.9_wp * dt_ros * & 458 errmax**pshrnk ), & 459 shrnk * ABS( dt_ros ) ), dt_ros ) 424 460 ENDIF 425 461 … … 429 465 !-- Calculate next internal time step 430 466 IF ( errmax > errcon ) THEN 431 dt_ros_next = 0.9 * dt_ros * errmax**pgrow467 dt_ros_next = 0.9_wp * dt_ros * errmax**pgrow 432 468 ELSE 433 469 dt_ros_next = grow * dt_ros … … 447 483 particles(n)%rvar1 = dt_ros_next 448 484 449 new_r = r_ros485 new_r(n) = r_ros 450 486 ! 451 487 !-- Radius should not fall below 1E-8 because Rosenbrock method may 452 488 !-- lead to errors otherwise 453 new_r = MAX( new_r, 1.0E-8_wp )489 new_r(n) = MAX( new_r(n), 1.0E-8_wp ) 454 490 ! 455 491 !-- Check if calculated droplet radius change is reasonable since in … … 457 493 !-- secondary solutions which are physically unlikely. 458 494 !-- Due to the solution effect the droplets may grow for relative 459 !-- humidities below 100%, but change of radius should not be too large.460 !-- In case of unreasonable droplet growth the Rosenbrock method is461 !-- recalculated using a smaller initial time step.495 !-- humidities below 100%, but change of radius should not be too 496 !-- large. In case of unreasonable droplet growth the Rosenbrock 497 !-- method is recalculated using a smaller initial time step. 462 498 !-- Limiting values are tested for droplets down to 1.0E-7 463 IF ( new_r - particles(n)%radius >= 3.0E-7.AND. &464 e_a /e_s < 0.97) THEN499 IF ( new_r(n) - particles(n)%radius >= 3.0E-7_wp .AND. & 500 e_a(n)/e_s(n) < 0.97_wp ) THEN 465 501 ros_count = ros_count + 1 466 502 repeat = .TRUE. … … 471 507 ENDIF 472 508 473 delta_r = new_r - particles(n)%radius509 delta_r = new_r(n) - particles(n)%radius 474 510 475 511 ! … … 477 513 !-- volume (this is needed later in lpm_calc_liquid_water_content for 478 514 !-- calculating the release of latent heat) 479 i = ( particles(n)%x + 0.5 * dx ) * ddx480 j = ( particles(n)%y + 0.5 * dy ) * ddy481 k = particles(n)%z / dz + 1 + offset_ocean_nzt_m1515 i = ip 516 j = jp 517 k = kp 482 518 ! only exact if equidistant 483 519 484 ql_c(k,j,i) = ql_c(k,j,i) + particles(n)%weight_factor * &485 rho_l * 1.33333333 * pi *&486 ( new_r **3 - particles(n)%radius**3 ) /&520 ql_c(k,j,i) = ql_c(k,j,i) + particles(n)%weight_factor * & 521 rho_l * 1.33333333_wp * pi * & 522 ( new_r(n)**3 - particles(n)%radius**3 ) / & 487 523 ( rho_surface * dx * dy * dz ) 488 IF ( ql_c(k,j,i) > 100.0 ) THEN524 IF ( ql_c(k,j,i) > 100.0_wp ) THEN 489 525 WRITE( message_string, * ) 'k=',k,' j=',j,' i=',i, & 490 526 ' ql_c=',ql_c(k,j,i), ' &part(',n,')%wf=', & … … 495 531 ! 496 532 !-- Change the droplet radius 497 IF ( ( new_r - particles(n)%radius ) < 0.0 .AND. new_r < 0.0 )&498 THEN499 WRITE( message_string, * ) '#1 k=',k,' j=',j,' i=',i, &500 ' e_s=',e_s , ' e_a=',e_a,' t_int=',t_int,&501 ' &delta_r=',delta_r, &533 IF ( ( new_r(n) - particles(n)%radius ) < 0.0_wp .AND. & 534 new_r(n) < 0.0_wp ) THEN 535 WRITE( message_string, * ) '#1 k=',k,' j=',j,' i=',i, & 536 ' e_s=',e_s(n), ' e_a=',e_a(n),' t_int=',t_int(n), & 537 ' &delta_r=',delta_r, & 502 538 ' particle_radius=',particles(n)%radius 503 539 CALL message( 'lpm_droplet_condensation', 'PA0144', 2, 2, -1, 6, 1 ) … … 507 543 !-- Sum up the total volume of liquid water (needed below for 508 544 !-- re-calculating the weighting factors) 509 ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * new_r **3510 511 particles(n)%radius = new_r 545 ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * new_r(n)**3 546 547 particles(n)%radius = new_r(n) 512 548 513 549 ! 514 550 !-- Determine radius class of the particle needed for collision 515 IF ( ( hall_kernel .OR. wang_kernel ) .AND. use_kernel_tables ) &551 IF ( ( hall_kernel .OR. wang_kernel ) .AND. use_kernel_tables ) & 516 552 THEN 517 particles(n)%class = ( LOG( new_r ) - rclass_lbound ) /&518 ( rclass_ubound - rclass_lbound ) * &553 particles(n)%class = ( LOG( new_r(n) ) - rclass_lbound ) / & 554 ( rclass_ubound - rclass_lbound ) * & 519 555 radius_classes 520 556 particles(n)%class = MIN( particles(n)%class, radius_classes )
Note: See TracChangeset
for help on using the changeset viewer.
