Changeset 1890 for palm/trunk/SOURCE
- Timestamp:
- Apr 22, 2016 8:52:11 AM (9 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 2 edited
-
lpm_droplet_condensation.f90 (modified) (11 diffs)
-
lpm_init.f90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lpm_droplet_condensation.f90
r1872 r1890 19 19 ! Current revisions: 20 20 ! ------------------ 21 ! 21 ! Some improvements of the Rosenbrock method. If the Rosenbrock method needs more 22 ! than 40 iterations to find a sufficient time setp, the model is not aborted. 23 ! This might lead to small erros. But they can be assumend as negligible, since 24 ! the maximum timesetp of the Rosenbrock method after 40 iterations will be 25 ! smaller than 10^-16 s. 22 26 ! 23 27 ! Former revisions: … … 131 135 IMPLICIT NONE 132 136 133 INTEGER(iwp) :: i !<134 137 INTEGER(iwp) :: ip !< 135 138 INTEGER(iwp) :: internal_timestep_count !< 136 INTEGER(iwp) :: j !<137 139 INTEGER(iwp) :: jp !< 138 140 INTEGER(iwp) :: jtry !< 139 INTEGER(iwp) :: k !<140 141 INTEGER(iwp) :: kp !< 141 142 INTEGER(iwp) :: n !< 142 143 INTEGER(iwp) :: ros_count !< 143 144 144 INTEGER(iwp), PARAMETER :: maxtry = 100!<145 INTEGER(iwp), PARAMETER :: maxtry = 40 !< 145 146 146 147 LOGICAL :: repeat !< … … 328 329 !-- Internal time step should not be > 1.0E-2 and < 1.0E-6 329 330 IF ( dt_ros_next > 1.0E-2_wp ) THEN 330 dt_ros_next = 1.0E- 3_wp331 dt_ros_next = 1.0E-2_wp 331 332 ELSEIF ( dt_ros_next < 1.0E-6_wp ) THEN 332 333 dt_ros_next = 1.0E-6_wp … … 366 367 drdt_ini = drdt 367 368 dt_ros_sum_ini = dt_ros_sum 368 r_ros_ini = MAX( r_ros, particles(n)%rvar2 )369 r_ros_ini = r_ros 369 370 370 371 ! … … 380 381 381 382 IF ( jtry == maxtry+1 ) THEN 382 message_string = 'maxtry > 100 in Rosenbrock method'383 CALL message( 'lpm_droplet_condensation', 'PA0347', 2, &384 2, 0, 6, 0 )383 message_string = 'maxtry > 40 in Rosenbrock method' 384 CALL message( 'lpm_droplet_condensation', 'PA0347', 0, & 385 1, 0, 6, 0 ) 385 386 ENDIF 386 387 387 388 aa = 1.0_wp / ( gam * dt_ros ) - d2rdtdr 388 389 g1 = drdt_ini / aa 389 r_ros = MAX( r_ros_ini + a21 * g1, particles(n)%rvar2 )390 r_ros = r_ros_ini + a21 * g1 390 391 drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0_wp - & 391 392 afactor / r_ros + & … … 395 396 g2 = ( drdt + c21 * g1 / dt_ros )& 396 397 / aa 397 r_ros = MAX( r_ros_ini + a31 * g1 + a32 * g2, particles(n)%rvar2 )398 r_ros = r_ros_ini + a31 * g1 + a32 * g2 398 399 drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0_wp - & 399 400 afactor / r_ros + & … … 405 406 g4 = ( drdt + & 406 407 ( c41 * g1 + c42 * g2 + c43 * g3 ) / dt_ros ) / aa 407 r_ros = MAX( r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + & 408 b4 * g4, particles(n)%rvar2 ) 408 r_ros = r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + b4 * g4 409 409 410 410 dt_ros_sum = dt_ros_sum_ini + dt_ros … … 426 426 EXIT 427 427 ELSE 428 dt_ros = SIGN( MAX( ABS( 0.9_wp * dt_ros * & 429 errmax**pshrnk ), & 430 shrnk * ABS( dt_ros ) ), dt_ros ) 428 dt_ros = MAX( ABS( 0.9_wp * dt_ros * errmax**pshrnk ), & 429 shrnk * ABS( dt_ros ) ) 431 430 ENDIF 432 431 … … 483 482 !-- volume (this is needed later in lpm_calc_liquid_water_content for 484 483 !-- calculating the release of latent heat) 485 i = ip 486 j = jp 487 k = kp 488 ! only exact if equidistant 489 490 ql_c(k,j,i) = ql_c(k,j,i) + particles(n)%weight_factor * & 484 ql_c(kp,jp,ip) = ql_c(kp,jp,ip) + particles(n)%weight_factor * & 491 485 rho_l * 1.33333333_wp * pi * & 492 486 ( new_r(n)**3 - particles(n)%radius**3 ) / & 493 487 ( rho_surface * dx * dy * dz ) 494 IF ( ql_c(k ,j,i) > 100.0_wp ) THEN495 WRITE( message_string, * ) 'k=',k ,' j=',j,' i=',i, &496 ' ql_c=',ql_c(k ,j,i), ' &part(',n,')%wf=', &488 IF ( ql_c(kp,jp,ip) > 100.0_wp ) THEN 489 WRITE( message_string, * ) 'k=',kp,' j=',jp,' i=',ip, & 490 ' ql_c=',ql_c(kp,jp,ip), ' &part(',n,')%wf=', & 497 491 particles(n)%weight_factor,' delta_r=',delta_r 498 492 CALL message( 'lpm_droplet_condensation', 'PA0143', 2, 2, -1, 6, 1 ) … … 501 495 ! 502 496 !-- Change the droplet radius 503 IF ( ( new_r(n) - particles(n)%radius ) < 0.0_wp .AND. 497 IF ( ( new_r(n) - particles(n)%radius ) < 0.0_wp .AND. & 504 498 new_r(n) < 0.0_wp ) THEN 505 WRITE( message_string, * ) '#1 k=',k ,' j=',j,' i=',i,&499 WRITE( message_string, * ) '#1 k=',kp,' j=',jp,' i=',ip, & 506 500 ' e_s=',e_s, ' e_a=',e_a,' t_int=',t_int, & 507 ' &delta_r=',delta_r, 501 ' &delta_r=',delta_r, & 508 502 ' particle_radius=',particles(n)%radius 509 503 CALL message( 'lpm_droplet_condensation', 'PA0144', 2, 2, -1, 6, 1 ) … … 513 507 !-- Sum up the total volume of liquid water (needed below for 514 508 !-- re-calculating the weighting factors) 515 ql_v(k ,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * new_r(n)**3509 ql_v(kp,jp,ip) = ql_v(kp,jp,ip) + particles(n)%weight_factor * new_r(n)**3 516 510 517 511 particles(n)%radius = new_r(n) -
palm/trunk/SOURCE/lpm_init.f90
r1874 r1890 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Initialization of aerosol equilibrium radius not possible in supersaturated 22 ! environments. Therefore, a maximum supersaturation of -1 % is assumed during 23 ! initialization. 22 24 ! 23 25 ! Former revisions: … … 811 813 REAL(wp) :: weight_sum !< sum of all weighting factors 812 814 813 INTEGER(iwp), DIMENSION( :,:,:), INTENT(IN) :: local_start !<815 INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) :: local_start !< 814 816 815 817 INTEGER(iwp) :: n !< … … 980 982 981 983 ! 982 !-- The formula is only valid for subsaturated environments. In (super-)983 !-- s aturated air, the inital radius is used.984 IF ( e_a / e_s < 1.0_wp ) THEN984 !-- The formula is only valid for subsaturated environments. For 985 !-- supersaturations higher than -1 %, the supersaturation is set to -1%. 986 IF ( e_a / e_s < 0.99_wp ) THEN 985 987 986 988 DO n = local_start(kp,jp,ip), number_of_particles !only new particles … … 995 997 ENDDO 996 998 999 ELSE 1000 1001 DO n = local_start(kp,jp,ip), number_of_particles !only new particles 1002 1003 bfactor = vanthoff * molecular_weight_of_water * & 1004 rho_s * particles(n)%rvar2**3 / & 1005 ( molecular_weight_of_solute * rho_l ) 1006 particles(n)%radius = particles(n)%rvar2 * ( bfactor / & 1007 particles(n)%rvar2**3 )**(1.0_wp/3.0_wp) *& 1008 0.01_wp**(-1.0_wp/3.0_wp) 1009 1010 ENDDO 1011 997 1012 ENDIF 998 1013
Note: See TracChangeset
for help on using the changeset viewer.
