- Timestamp:
- Sep 19, 2012 2:30:36 PM (12 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 7 edited
-
data_output_2d.f90 (modified) (2 diffs)
-
data_output_3d.f90 (modified) (2 diffs)
-
data_output_mask.f90 (modified) (7 diffs)
-
flow_statistics.f90 (modified) (10 diffs)
-
lpm_collision_kernels.f90 (modified) (33 diffs)
-
production_e.f90 (modified) (24 diffs)
-
sum_up_3d_data.f90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/data_output_2d.f90
r979 r1007 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! Bugfix: missing calculation of ql_vp added 7 7 ! 8 8 ! Former revisions: … … 428 428 CASE ( 'ql_vp_xy', 'ql_vp_xz', 'ql_vp_yz' ) 429 429 IF ( av == 0 ) THEN 430 IF ( simulated_time >= particle_advection_start ) THEN 431 DO i = nxl, nxr 432 DO j = nys, nyn 433 DO k = nzb, nzt+1 434 psi = prt_start_index(k,j,i) 435 DO n = psi, psi+prt_count(k,j,i)-1 436 tend(k,j,i) = tend(k,j,i) + & 437 particles(n)%weight_factor / & 438 prt_count(k,j,i) 439 ENDDO 440 ENDDO 441 ENDDO 442 ENDDO 443 CALL exchange_horiz( tend, nbgp ) 444 ELSE 445 tend = 0.0 446 END IF 447 DO i = nxlg, nxrg 448 DO j = nysg, nyng 449 DO k = nzb, nzt+1 450 local_pf(i,j,k) = tend(k,j,i) 451 ENDDO 452 ENDDO 453 ENDDO 454 resorted = .TRUE. 455 ELSE 456 CALL exchange_horiz( ql_vp_av, nbgp ) 430 457 to_be_resorted => ql_vp 431 ELSE432 to_be_resorted => ql_vp_av433 458 ENDIF 434 459 IF ( mode == 'xy' ) level_z = zu -
palm/trunk/SOURCE/data_output_3d.f90
r791 r1007 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! Bugfix: missing calculation of ql_vp added 7 7 ! 8 8 ! Former revisions: … … 295 295 CASE ( 'ql_vp' ) 296 296 IF ( av == 0 ) THEN 297 to_be_resorted => ql_vp 298 ELSE 297 DO i = nxl, nxr 298 DO j = nys, nyn 299 DO k = nzb, nz_do3d 300 psi = prt_start_index(k,j,i) 301 DO n = psi, psi+prt_count(k,j,i)-1 302 tend(k,j,i) = tend(k,j,i) + & 303 particles(n)%weight_factor / & 304 prt_count(k,j,i) 305 ENDDO 306 ENDDO 307 ENDDO 308 ENDDO 309 CALL exchange_horiz( tend, nbgp ) 310 DO i = nxlg, nxrg 311 DO j = nysg, nyng 312 DO k = nzb, nz_do3d 313 local_pf(i,j,k) = tend(k,j,i) 314 ENDDO 315 ENDDO 316 ENDDO 317 resorted = .TRUE. 318 ELSE 319 CALL exchange_horiz( ql_vp_av, nbgp ) 299 320 to_be_resorted => ql_vp_av 300 321 ENDIF -
palm/trunk/SOURCE/data_output_mask.f90
r772 r1007 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! Bugfix: calculation of pr must depend on the particle weighting factor, 7 ! missing calculation of ql_vp added 6 8 ! 7 9 ! Former revisions: … … 159 161 s_r4 = 0.0 160 162 DO n = psi, psi+prt_count(k,j,i)-1 161 s_r3 = s_r3 + particles(n)%radius**3 162 s_r4 = s_r4 + particles(n)%radius**4 163 s_r3 = s_r3 + particles(n)%radius**3 * & 164 particles(n)%weight_factor 165 s_r4 = s_r4 + particles(n)%radius**4 * & 166 particles(n)%weight_factor 163 167 ENDDO 164 168 IF ( s_r3 /= 0.0 ) THEN … … 213 217 to_be_resorted => q_av 214 218 ENDIF 215 219 216 220 CASE ( 'ql' ) 217 221 IF ( av == 0 ) THEN … … 237 241 CASE ( 'ql_vp' ) 238 242 IF ( av == 0 ) THEN 239 to_be_resorted => ql_vp 240 ELSE 243 DO i = nxl, nxr 244 DO j = nys, nyn 245 DO k = nzb, nz_do3d 246 psi = prt_start_index(k,j,i) 247 DO n = psi, psi+prt_count(k,j,i)-1 248 tend(k,j,i) = tend(k,j,i) + & 249 particles(n)%weight_factor / & 250 prt_count(k,j,i) 251 ENDDO 252 ENDDO 253 ENDDO 254 ENDDO 255 CALL exchange_horiz( tend, nbgp ) 256 DO i = 1, mask_size_l(mid,1) 257 DO j = 1, mask_size_l(mid,2) 258 DO k = 1, mask_size_l(mid,3) 259 local_pf(i,j,k) = tend(mask_k(mid,k), & 260 mask_j(mid,j),mask_i(mid,i)) 261 ENDDO 262 ENDDO 263 ENDDO 264 resorted = .TRUE. 265 ELSE 266 CALL exchange_horiz( ql_vp_av, nbgp ) 241 267 to_be_resorted => ql_vp_av 242 268 ENDIF … … 264 290 to_be_resorted => rho_av 265 291 ENDIF 266 292 267 293 CASE ( 's' ) 268 294 IF ( av == 0 ) THEN … … 271 297 to_be_resorted => s_av 272 298 ENDIF 273 299 274 300 CASE ( 'sa' ) 275 301 IF ( av == 0 ) THEN … … 278 304 to_be_resorted => sa_av 279 305 ENDIF 280 306 281 307 CASE ( 'u' ) 282 308 IF ( av == 0 ) THEN -
palm/trunk/SOURCE/flow_statistics.f90
r802 r1007 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! Calculation of buoyancy flux for humidity in case of WS-scheme is now using 7 ! turbulent fluxes of WS-scheme 8 ! Bugfix: Calculation of subgridscale buoyancy flux for humidity and cloud 9 ! droplets at nzb and nzt added 6 10 ! 7 11 ! Former revisions: … … 151 155 'timestep' 152 156 CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 ) 153 157 154 158 ENDIF 155 159 … … 176 180 IF ( ws_scheme_mom .AND. sr == 0 ) THEN 177 181 178 ! 182 ! 179 183 !-- According to the Neumann bc for the horizontal velocity components, 180 184 !-- the corresponding fluxes has to satisfiy the same bc. 181 185 IF ( ocean ) THEN 182 186 sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:) 183 sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:) 187 sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:) 184 188 ENDIF 185 189 186 190 DO i = 0, threads_per_task-1 187 ! 191 ! 188 192 !-- Swap the turbulent quantities evaluated in advec_ws. 189 193 sums_l(:,13,i) = sums_wsus_ws_l(:,i) ! w*u* … … 453 457 ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr) 454 458 ENDIF 455 459 456 460 ! 457 461 !-- Higher moments … … 578 582 * ( q(k+1,j,i) - q(k,j,i) ) & 579 583 * ddzu(k+1) * rmask(j,i,sr) 584 580 585 IF ( cloud_physics ) THEN 581 586 sums_l(k,51,tn) = sums_l(k,51,tn) & … … 618 623 sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + & 619 624 qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 625 sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + ( & 626 ( 1.0 + 0.61 * q(nzb,j,i) ) * & 627 shf(j,i) + 0.61 * pt(nzb,j,i) * & 628 qsws(j,i) ) 629 IF ( cloud_droplets ) THEN 630 sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + ( & 631 ( 1.0 + 0.61 * q(nzb,j,i) - & 632 ql(nzb,j,i) ) * shf(j,i) + & 633 0.61 * pt(nzb,j,i) * qsws(j,i) ) 634 ENDIF 620 635 IF ( cloud_physics ) THEN 621 sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + ( &622 ( 1.0 + 0.61 * q(nzb,j,i) ) * &623 shf(j,i) + 0.61 * pt(nzb,j,i) * &624 qsws(j,i) &625 )626 636 ! 627 637 !-- Formula does not work if ql(nzb) /= 0.0 … … 657 667 sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + & 658 668 qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 669 sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + ( & 670 ( 1.0 + 0.61 * q(nzt,j,i) ) * & 671 tswst(j,i) + 0.61 * pt(nzt,j,i) * & 672 qswst(j,i) ) 673 IF ( cloud_droplets ) THEN 674 sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + ( & 675 ( 1.0 + 0.61 * q(nzt,j,i) - & 676 ql(nzt,j,i) ) * tswst(j,i) + & 677 0.61 * pt(nzt,j,i) * qswst(j,i) ) 678 ENDIF 659 679 IF ( cloud_physics ) THEN 660 sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + ( &661 ( 1.0 + 0.61 * q(nzt,j,i) ) * &662 tswst(j,i) + 0.61 * pt(nzt,j,i) * &663 qsws(j,i) &664 )665 680 ! 666 681 !-- Formula does not work if ql(nzb) /= 0.0 … … 714 729 !-- content 715 730 IF ( humidity ) THEN 716 pts = 0.5 * ( vpt(k,j,i) - hom(k,1,44,sr) + & 717 vpt(k+1,j,i) - hom(k+1,1,44,sr) ) 718 sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * & 719 rmask(j,i,sr) 720 721 IF ( cloud_physics .OR. cloud_droplets ) THEN 722 pts = 0.5 * & 723 ( ( q(k,j,i) - ql(k,j,i) ) - hom(k,1,42,sr) & 724 + ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) ) 725 sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) * & 731 IF ( cloud_physics .OR. cloud_droplets ) THEN 732 pts = 0.5 * ( vpt(k,j,i) - hom(k,1,44,sr) + & 733 vpt(k+1,j,i) - hom(k+1,1,44,sr) ) 734 sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * & 726 735 rmask(j,i,sr) 727 736 sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * & 728 737 rmask(j,i,sr) 729 ENDIF 730 ENDIF 731 738 ELSE 739 IF( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN 740 pts = 0.5 * ( vpt(k,j,i) - hom(k,1,44,sr) + & 741 vpt(k+1,j,i) - hom(k+1,1,44,sr) ) 742 sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * & 743 rmask(j,i,sr) 744 ELSE IF ( ws_scheme_sca .AND. sr == 0 ) THEN 745 sums_l(k,46,tn) = ( 1.0 + 0.61 * hom(k,1,41,sr) ) * & 746 sums_l(k,17,tn) + & 747 0.61 * hom(k,1,4,sr) * sums_l(k,49,tn) 748 END IF 749 END IF 750 ENDIF 732 751 ! 733 752 !-- Passive scalar flux … … 762 781 vst = 0.5 * ( v(k,j,i) - hom(k,1,2,sr) + & 763 782 v(k+1,j,i) - hom(k+1,1,2,sr) ) 764 ! 783 ! 765 784 !-- Momentum flux w*u* 766 785 sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 * & … … 1002 1021 ENDDO 1003 1022 ENDIF 1023 1004 1024 ! 1005 1025 !-- Collect horizontal average in hom. -
palm/trunk/SOURCE/lpm_collision_kernels.f90
r850 r1007 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! converted all units to SI units and replaced some parameters by corresponding 7 ! PALM parameters 8 ! Bugfix: factor in calculation of enhancement factor for collision efficencies 9 ! changed from 10. to 1.0 7 10 ! 8 11 ! Former revisions: … … 64 67 65 68 PUBLIC ckernel, collision_efficiency_rogers, init_kernels, & 66 rclass_lbound, rclass_ubound, recalculate_kernel 69 rclass_lbound, rclass_ubound, recalculate_kernel 67 70 68 71 REAL :: epsilon, eps2, rclass_lbound, rclass_ubound, urms, urms2 … … 126 129 ! ENDIF 127 130 ENDDO 128 ! 129 !-- Collision routines expect radius to be in cm 130 radclass = radclass * 100.0 131 132 ! 133 !-- Set the class bounds for dissipation in interval [0, 1000] cm**2/s**3 131 132 ! 133 !-- Set the class bounds for dissipation in interval [0.0, 0.1] m**2/s**3 134 134 DO i = 1, dissipation_classes 135 epsclass(i) = 1000.0* REAL( i ) / dissipation_classes135 epsclass(i) = 0.1 * REAL( i ) / dissipation_classes 136 136 ! IF ( myid == 0 ) THEN 137 137 ! PRINT*, 'i=', i, ' eps = ', epsclass(i) … … 148 148 149 149 epsilon = epsclass(k) 150 urms = 2 02.0 * ( epsilon / 400.0)**( 1.0 / 3.0 )150 urms = 2.02 * ( epsilon / 0.04 )**( 1.0 / 3.0 ) 151 151 152 152 CALL turbsd … … 183 183 184 184 PRINT*, '*** Hall kernel' 185 WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E4, i = 1,radius_classes ) 185 WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E6, & 186 i = 1,radius_classes ) 186 187 DO j = 1, radius_classes 187 WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j), ( hkernel(i,j), i = 1,radius_classes ) 188 WRITE ( *,'(F4.0,1X,20(F8.4,1X))' ) radclass(j), & 189 ( hkernel(i,j), i = 1,radius_classes ) 188 190 ENDDO 189 191 … … 200 202 201 203 PRINT*, '*** epsilon = ', epsclass(k) 202 WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E4, i = 1,radius_classes ) 204 WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E6, & 205 i = 1,radius_classes ) 203 206 DO j = 1, radius_classes 204 ! WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E4, ( ckernel(i,j,k), i = 1,radius_classes ) 205 WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E4, ( hwratio(i,j), i = 1,radius_classes ) 207 ! WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E6, & 208 ! ( ckernel(i,j,k), i = 1,radius_classes ) 209 WRITE ( *,'(F4.0,1X,20(F8.4,1X))' ) radclass(j)*1.0E6, & 210 ( hwratio(i,j), i = 1,radius_classes ) 206 211 ENDDO 207 212 ENDDO … … 210 215 211 216 DEALLOCATE( ec, ecf, epsclass, gck, hkernel, winf ) 212 213 ckernel = ckernel * 1.0E-6 ! kernel is needed in m**3/s214 217 215 218 ELSEIF( collision_kernel == 'hall' .OR. collision_kernel == 'wang' ) & … … 249 252 250 253 ! 251 !-- Store particle radii on the radclass array. Collision routines 252 !-- expect radii to be in cm. 253 radclass(1:radius_classes) = particles(pstart:pend)%radius * 100.0 254 !-- Store particle radii on the radclass array 255 radclass(1:radius_classes) = particles(pstart:pend)%radius 254 256 255 257 IF ( wang_kernel ) THEN 256 epsilon = diss(k1,j1,i1) * 1.0E4 ! dissipation rate in cm**2/s**-3258 epsilon = diss(k1,j1,i1) ! dissipation rate in m**2/s**3 257 259 ELSE 258 260 epsilon = 0.0 259 261 ENDIF 260 urms = 2 02.0 * ( epsilon / 400.0)**( 0.33333333333 )261 262 IF ( wang_kernel .AND. epsilon > 0.001) THEN262 urms = 2.02 * ( epsilon / 0.04 )**( 0.33333333333 ) 263 264 IF ( wang_kernel .AND. epsilon > 1.0E-7 ) THEN 263 265 ! 264 266 !-- Call routines to calculate efficiencies for the Wang kernel … … 294 296 ENDIF 295 297 296 ckernel = ckernel * 1.0E-6 ! kernel is needed in m**3/s297 298 298 DEALLOCATE( ec, radclass, winf ) 299 299 … … 314 314 USE particle_attributes 315 315 USE arrays_3d 316 USE control_parameters 316 317 317 318 IMPLICIT NONE … … 326 327 v1xysq, v2, v2xysq, wrfin, wrgrav2, wrtur2xy, xx, yy, z 327 328 328 REAL, SAVE :: airdens, airvisc, anu, gravity, waterdens329 330 329 REAL, DIMENSION(1:radius_classes) :: st, tau 331 330 … … 336 335 337 336 first = .FALSE. 338 airvisc = 0.1818 ! dynamic viscosity in mg/cm*s 339 airdens = 1.2250 ! air density in mg/cm**3 340 waterdens = 1000.0 ! water density in mg/cm**3 341 gravity = 980.6650 ! in cm/s**2 342 anu = airvisc/airdens ! kinetic viscosity in cm**2/s 343 344 ENDIF 345 346 lambda = urms * SQRT( 15.0 * anu / epsilon ) ! in cm 347 lambda_re = urms**2 * SQRT( 15.0 / epsilon / anu ) 337 338 ENDIF 339 340 lambda = urms * SQRT( 15.0 * molecular_viscosity / epsilon ) ! in m 341 lambda_re = urms**2 * SQRT( 15.0 / epsilon / molecular_viscosity ) 348 342 tl = urms**2 / epsilon ! in s 349 lf = 0.5 * urms**3 / epsilon ! in cm350 tauk = SQRT( anu / epsilon )! in s351 eta = ( anu**3 / epsilon )**0.25 ! in cm352 vk = eta / tauk ! in cm/s343 lf = 0.5 * urms**3 / epsilon ! in m 344 tauk = SQRT( molecular_viscosity / epsilon ) ! in s 345 eta = ( molecular_viscosity**3 / epsilon )**0.25 ! in m 346 vk = eta / tauk 353 347 354 348 ao = ( 11.0 + 7.0 * lambda_re ) / ( 205.0 + lambda_re ) 355 349 tt = SQRT( 2.0 * lambda_re / ( SQRT( 15.0 ) * ao ) ) * tauk ! in s 356 350 357 CALL fallg ! gives winf in cm/s351 CALL fallg ! gives winf in m/s 358 352 359 353 DO i = 1, radius_classes 360 tau(i) = winf(i) / g ravity! in s354 tau(i) = winf(i) / g ! in s 361 355 st(i) = tau(i) / tauk 362 356 ENDDO … … 368 362 bbb = SQRT( 1.0 - 2.0 * be**2 ) 369 363 d1 = ( 1.0 + bbb ) / ( 2.0 * bbb ) 370 e1 = lf * ( 1.0 + bbb ) * 0.5 ! in cm364 e1 = lf * ( 1.0 + bbb ) * 0.5 ! in m 371 365 d2 = ( 1.0 - bbb ) * 0.5 / bbb 372 e2 = lf * ( 1.0 - bbb ) * 0.5 ! in cm366 e2 = lf * ( 1.0 - bbb ) * 0.5 ! in m 373 367 ccc = SQRT( 1.0 - 2.0 * z**2 ) 374 368 b1 = ( 1.0 + ccc ) * 0.5 / ccc … … 379 373 DO i = 1, radius_classes 380 374 381 v1 = winf(i) ! in cm/s375 v1 = winf(i) ! in m/s 382 376 t1 = tau(i) ! in s 383 377 384 378 DO j = 1, i 385 rrp = radclass(i) + radclass(j) ! radius in cm386 v2 = winf(j) ! in cm/s379 rrp = radclass(i) + radclass(j) 380 v2 = winf(j) ! in m/s 387 381 t2 = tau(j) ! in s 388 382 389 v1xysq = b1 * d1 * phi (c1,e1,v1,t1) - b1 * d2 * phi(c1,e2,v1,t1) &390 - b2 * d1 * phi (c2,e1,v1,t1) + b2 * d2 * phi(c2,e2,v1,t1)391 v1xysq = v1xysq * urms**2 / t1 ! in cm**2/s**2392 vrms1xy = SQRT( v1xysq ) ! in cm/s393 394 v2xysq = b1 * d1 * phi (c1,e1,v2,t2) - b1 * d2 * phi(c1,e2,v2,t2) &395 - b2 * d1 * phi (c2,e1,v2,t2) + b2 * d2 * phi(c2,e2,v2,t2)396 v2xysq = v2xysq * urms**2 / t2 ! in cm**2/s**2397 vrms2xy = SQRT( v2xysq ) ! in cm/s383 v1xysq = b1 * d1 * phi_w(c1,e1,v1,t1) - b1 * d2 * phi_w(c1,e2,v1,t1) & 384 - b2 * d1 * phi_w(c2,e1,v1,t1) + b2 * d2 * phi_w(c2,e2,v1,t1) 385 v1xysq = v1xysq * urms**2 / t1 ! in m**2/s**2 386 vrms1xy = SQRT( v1xysq ) ! in m/s 387 388 v2xysq = b1 * d1 * phi_w(c1,e1,v2,t2) - b1 * d2 * phi_w(c1,e2,v2,t2) & 389 - b2 * d1 * phi_w(c2,e1,v2,t2) + b2 * d2 * phi_w(c2,e2,v2,t2) 390 v2xysq = v2xysq * urms**2 / t2 ! in m**2/s**2 391 vrms2xy = SQRT( v2xysq ) ! in m/s 398 392 399 393 IF ( winf(i) >= winf(j) ) THEN … … 414 408 b2 * d2* zhi(c2,e2,v1,t1,v2,t2) 415 409 fr = d1 * EXP( -rrp / e1 ) - d2 * EXP( -rrp / e2 ) 416 v1v2xy = v1v2xy * fr * urms**2 / tau(i) / tau(j) ! in cm**2/s**2417 wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0 * v1v2xy ! in cm**2/s**2410 v1v2xy = v1v2xy * fr * urms**2 / tau(i) / tau(j) ! in m**2/s**2 411 wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0 * v1v2xy ! in m**2/s**2 418 412 IF ( wrtur2xy < 0.0 ) wrtur2xy = 0.0 419 413 wrgrav2 = pi / 8.0 * ( winf(j) - winf(i) )**2 420 wrfin = SQRT( ( 2.0 / pi ) * ( wrtur2xy + wrgrav2) ) ! in cm/s414 wrfin = SQRT( ( 2.0 / pi ) * ( wrtur2xy + wrgrav2) ) ! in m/s 421 415 422 416 ! … … 433 427 yy = 0.1886 * EXP( 20.306 / lambda_re ) 434 428 435 c1_gr = xx / ( g ravity/ vk * tauk )**yy436 437 ao_gr = ao + ( pi / 8.0) * ( g ravity/ vk * tauk )**2429 c1_gr = xx / ( g / vk * tauk )**yy 430 431 ao_gr = ao + ( pi / 8.0) * ( g / vk * tauk )**2 438 432 fao_gr = 20.115 * SQRT( ao_gr / lambda_re ) 439 433 rc = SQRT( fao_gr * ABS( st(j) - st(i) ) ) * eta ! in cm … … 452 446 453 447 !------------------------------------------------------------------------------! 454 ! phi as a function455 !------------------------------------------------------------------------------! 456 REAL FUNCTION phi ( a, b, vsett, tau0 )448 ! phi_w as a function 449 !------------------------------------------------------------------------------! 450 REAL FUNCTION phi_w( a, b, vsett, tau0 ) 457 451 458 452 IMPLICIT NONE … … 461 455 462 456 aa1 = 1.0 / tau0 + 1.0 / a + vsett / b 463 phi = 1.0 / aa1 - 0.5 * vsett / b / aa1**2 ! in s464 465 END FUNCTION phi 466 467 468 !------------------------------------------------------------------------------! 469 ! z etaas a function457 phi_w = 1.0 / aa1 - 0.5 * vsett / b / aa1**2 ! in s 458 459 END FUNCTION phi_w 460 461 462 !------------------------------------------------------------------------------! 463 ! zhi as a function 470 464 !------------------------------------------------------------------------------! 471 465 REAL FUNCTION zhi( a, b, vsett1, tau1, vsett2, tau2 ) … … 490 484 491 485 !------------------------------------------------------------------------------! 492 ! Calculation of terminal velocity winf 486 ! Calculation of terminal velocity winf following Equations 10-138 to 10-145 487 ! from (Pruppacher and Klett, 1997) 493 488 !------------------------------------------------------------------------------! 494 489 SUBROUTINE fallg … … 498 493 USE particle_attributes 499 494 USE arrays_3d 495 USE control_parameters 500 496 501 497 IMPLICIT NONE … … 505 501 LOGICAL, SAVE :: first = .TRUE. 506 502 507 REAL, SAVE :: cunh, eta, grav, phy, py, rhoa, rhow, sigma, stb, stok, &503 REAL, SAVE :: cunh, eta, phy, py, rho_a, sigma, stb, stok, & 508 504 t0, xlamb 509 505 … … 523 519 -0.542819E-1, 0.238449E-2 /) 524 520 525 eta = 1.818E-4 ! in poise = g/(cm s) 526 xlamb = 6.62E-6 ! in cm 527 rhow = 1.0 ! in g/cm**3 528 rhoa = 1.225E-3 ! in g/cm**3 529 grav = 980.665 ! in cm/s**2 530 cunh = 1.257 * xlamb ! in cm 531 t0 = 273.15 ! in K 532 sigma = 76.1 - 0.155 * ( 293.15 - t0 ) ! in N/m = g/s**2 533 stok = 2.0 * grav * ( rhow - rhoa ) / ( 9.0 * eta ) ! in 1/(cm s) 534 stb = 32.0 * rhoa * ( rhow - rhoa) * grav / (3.0 * eta * eta) 535 ! in 1/cm**3 536 phy = sigma**3 * rhoa**2 / ( eta**4 * grav * ( rhow - rhoa ) ) 521 ! 522 !-- Parameter values for p = 1013,25 hPa and T = 293,15 K 523 eta = 1.818E-5 ! in kg/(m s) 524 xlamb = 6.6E-8 ! in m 525 rho_a = 1.204 ! in kg/m**3 526 cunh = 1.26 * xlamb ! in m 527 sigma = 0.07363 ! in kg/s**2 528 stok = 2.0 * g * ( rho_l - rho_a ) / ( 9.0 * eta ) ! in 1/(m s) 529 stb = 32.0 * rho_a * ( rho_l - rho_a) * g / (3.0 * eta * eta) 530 phy = sigma**3 * rho_a**2 / ( eta**4 * g * ( rho_l - rho_a ) ) 537 531 py = phy**( 1.0 / 6.0 ) 538 532 … … 541 535 DO j = 1, radius_classes 542 536 543 IF ( radclass(j) <= 1.0E- 3) THEN544 545 winf(j) = stok * ( radclass(j)**2 + cunh * radclass(j) ) ! in cm/s546 547 ELSEIF ( radclass(j) > 1.0E- 3 .AND. radclass(j) <= 5.35E-2) THEN537 IF ( radclass(j) <= 1.0E-5 ) THEN 538 539 winf(j) = stok * ( radclass(j)**2 + cunh * radclass(j) ) 540 541 ELSEIF ( radclass(j) > 1.0E-5 .AND. radclass(j) <= 5.35E-4 ) THEN 548 542 549 543 x = LOG( stb * radclass(j)**3 ) … … 553 547 y = y + b(i) * x**(i-1) 554 548 ENDDO 555 556 xrey = ( 1.0 + cunh / radclass(j) ) * EXP( y ) 557 winf(j) = xrey * eta / ( 2.0 * rhoa * radclass(j) ) ! in cm/s 558 559 ELSEIF ( radclass(j) > 5.35E-2 ) THEN 560 561 IF ( radclass(j) > 0.35 ) THEN 562 bond = grav * ( rhow - rhoa ) * 0.35**2 / sigma 549 ! 550 !-- Note: this Eq. is wrong in (Pruppacher and Klett, 1997, p. 418) 551 !-- for correct version see (Beard, 1976) 552 xrey = ( 1.0 + cunh / radclass(j) ) * EXP( y ) 553 554 winf(j) = xrey * eta / ( 2.0 * rho_a * radclass(j) ) 555 556 ELSEIF ( radclass(j) > 5.35E-4 ) THEN 557 558 IF ( radclass(j) > 0.0035 ) THEN 559 bond = g * ( rho_l - rho_a ) * 0.0035**2 / sigma 563 560 ELSE 564 bond = grav * ( rhow - rhoa ) * radclass(j)**2 / sigma561 bond = g * ( rho_l - rho_a ) * radclass(j)**2 / sigma 565 562 ENDIF 566 563 … … 574 571 xrey = py * EXP( y ) 575 572 576 IF ( radclass(j) > 0. 35 ) THEN577 winf(j) = xrey * eta / ( 2.0 * rho a * 0.35 ) ! in cm/s573 IF ( radclass(j) > 0.0035 ) THEN 574 winf(j) = xrey * eta / ( 2.0 * rho_a * 0.0035 ) 578 575 ELSE 579 winf(j) = xrey * eta / ( 2.0 * rho a * radclass(j) ) ! in cm/s576 winf(j) = xrey * eta / ( 2.0 * rho_a * radclass(j) ) 580 577 ENDIF 581 578 … … 668 665 ! 669 666 !-- Calculate the radius class index of particles with respect to array r 667 !-- Radius has to be in µm 670 668 ALLOCATE( ira(1:radius_classes) ) 671 669 DO j = 1, radius_classes 672 particle_radius = radclass(j) * 1.0E 4 ! in microm670 particle_radius = radclass(j) * 1.0E6 673 671 DO k = 1, 15 674 672 IF ( particle_radius < r0(k) ) THEN … … 693 691 IF ( ir < 16 ) THEN 694 692 IF ( ir >= 2 ) THEN 695 pp = ( ( radclass(j) * 1.0E0 4) - r0(ir-1) ) / &693 pp = ( ( radclass(j) * 1.0E06 ) - r0(ir-1) ) / & 696 694 ( r0(ir) - r0(ir-1) ) 697 695 qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) ) … … 754 752 rat = (/ 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 /) 755 753 ! 756 !-- In 100 cm^2/s^3754 !-- for 100 cm**2/s**3 757 755 ecoll_100(:,1) = (/1.74, 1.74, 1.773, 1.49, 1.207, 1.207, 1.0 /) 758 756 ecoll_100(:,2) = (/1.46, 1.46, 1.421, 1.245, 1.069, 1.069, 1.0 /) … … 767 765 ecoll_100(:,11)= (/20.3, 20.3, 14.6 , 8.61, 2.60, 2.60 , 1.0 /) 768 766 ! 769 !-- In 400 cm^2/s^3767 !-- for 400 cm**2/s**3 770 768 ecoll_400(:,1) = (/4.976, 4.976, 3.593, 2.519, 1.445, 1.445, 1.0 /) 771 769 ecoll_400(:,2) = (/2.984, 2.984, 2.181, 1.691, 1.201, 1.201, 1.0 /) … … 784 782 ! 785 783 !-- Calculate the radius class index of particles with respect to array r0 784 !-- Radius has to be in µm 786 785 ALLOCATE( ira(1:radius_classes) ) 787 786 788 787 DO j = 1, radius_classes 789 particle_radius = radclass(j) * 1.0E 4 ! in microm788 particle_radius = radclass(j) * 1.0E6 790 789 DO k = 1, 7 791 790 IF ( particle_radius < r0(k) ) THEN … … 799 798 ! 800 799 !-- Two-dimensional linear interpolation of the collision efficiencies 800 !-- Radius has to be in µm 801 801 DO j = 1, radius_classes 802 802 DO i = 1, j … … 812 812 ENDDO 813 813 814 y1 = 1.0 ! 0 cm2/s3 815 ! 816 !-- 100 cm2/s3, 400 cm2/s3 814 y1 = 0.0001 ! for 0 m**2/s**3 815 817 816 IF ( ir < 8 ) THEN 818 817 IF ( ir >= 2 ) THEN 819 pp = ( radclass(j)*1.0E 4- r0(ir-1) ) / ( r0(ir) - r0(ir-1) )818 pp = ( radclass(j)*1.0E6 - r0(ir-1) ) / ( r0(ir) - r0(ir-1) ) 820 819 qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) ) 821 820 y2 = ( 1.0-pp ) * ( 1.0-qq ) * ecoll_100(ir-1,iq-1) + & 822 821 pp * ( 1.0-qq ) * ecoll_100(ir,iq-1) + & 823 qq * ( 1 0.-pp ) * ecoll_100(ir-1,iq) + &822 qq * ( 1.0-pp ) * ecoll_100(ir-1,iq) + & 824 823 pp * qq * ecoll_100(ir,iq) 825 824 y3 = ( 1.0-pp ) * ( 1.0-qq ) * ecoll_400(ir-1,iq-1) + & … … 838 837 ENDIF 839 838 ! 840 !-- Linear interpolation of dissipation rate in cm2/s3841 IF ( epsilon <= 100.0) THEN842 ecf(j,i) = ( epsilon - 100.0 ) / ( 0.0 - 100.0) * y1 &843 + ( epsilon - 0.0 ) / ( 100.0- 0.0 ) * y2844 ELSEIF ( epsilon <= 600.0) THEN845 ecf(j,i) = ( epsilon - 400.0 ) / ( 100.0 - 400.0) * y2 &846 + ( epsilon - 100.0 ) / ( 400.0 - 100.0) * y3839 !-- Linear interpolation of dissipation rate in m**2/s**3 840 IF ( epsilon <= 0.01 ) THEN 841 ecf(j,i) = ( epsilon - 0.01 ) / ( 0.0 - 0.01 ) * y1 & 842 + ( epsilon - 0.0 ) / ( 0.01 - 0.0 ) * y2 843 ELSEIF ( epsilon <= 0.06 ) THEN 844 ecf(j,i) = ( epsilon - 0.04 ) / ( 0.01 - 0.04 ) * y2 & 845 + ( epsilon - 0.01 ) / ( 0.04 - 0.01 ) * y3 847 846 ELSE 848 ecf(j,i) = ( 600.0 - 400.0 ) / ( 100.0 - 400.0) * y2 &849 + ( 600.0 - 100.0 ) / ( 400.0 - 100.0) * y3847 ecf(j,i) = ( 0.06 - 0.04 ) / ( 0.01 - 0.04 ) * y2 & 848 + ( 0.06 - 0.01 ) / ( 0.04 - 0.01 ) * y3 850 849 ENDIF 851 850 -
palm/trunk/SOURCE/production_e.f90
r941 r1007 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! Bugfix: calculation of buoyancy production has to consider the liquid water 7 ! mixing ratio in case of cloud droplets 7 8 ! 8 9 ! Former revisions: … … 164 165 165 166 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 166 167 167 168 ENDDO 168 169 ENDDO … … 189 190 190 191 IF ( wall_e_y(j,i) /= 0.0 ) THEN 191 ! 192 ! 192 193 !-- Inconsistency removed: as the thermal stratification is 193 194 !-- not taken into account for the evaluation of the wall … … 219 220 220 221 IF ( wall_e_x(j,i) /= 0.0 ) THEN 221 ! 222 ! 222 223 !-- Inconsistency removed: as the thermal stratification is 223 224 !-- not taken into account for the evaluation of the wall … … 275 276 276 277 IF ( wall_e_y(j,i) /= 0.0 ) THEN 277 ! 278 ! 278 279 !-- Inconsistency removed: as the thermal stratification 279 280 !-- is not taken into account for the evaluation of the … … 301 302 302 303 IF ( wall_e_x(j,i) /= 0.0 ) THEN 303 ! 304 ! 304 305 !-- Inconsistency removed: as the thermal stratification 305 306 !-- is not taken into account for the evaluation of the … … 416 417 417 418 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 418 419 419 420 ENDIF 420 421 … … 551 552 DO k = nzb_diff_s_inner(j,i), nzt_diff 552 553 553 IF ( .NOT. cloud_physics )THEN554 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 554 555 k1 = 1.0 + 0.61 * q(k,j,i) 555 556 k2 = 0.61 * pt(k,j,i) 556 ELSE 557 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * & 558 g / vpt(k,j,i) * & 559 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & 560 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & 561 ) * dd2zu(k) 562 ELSE IF ( cloud_physics ) THEN 557 563 IF ( ql(k,j,i) == 0.0 ) THEN 558 564 k1 = 1.0 + 0.61 * q(k,j,i) … … 568 574 k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) 569 575 ENDIF 570 ENDIF 571 572 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & 576 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * & 577 g / vpt(k,j,i) * & 573 578 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & 574 579 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & 575 580 ) * dd2zu(k) 581 ELSE IF ( cloud_droplets ) THEN 582 k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i) 583 k2 = 0.61 * pt(k,j,i) 584 tend(k,j,i) = tend(k,j,i) - & 585 kh(k,j,i) * g / vpt(k,j,i) * & 586 ( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) + & 587 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - & 588 pt(k,j,i) * ( ql(k+1,j,i) - & 589 ql(k-1,j,i) ) ) * dd2zu(k) 590 ENDIF 591 576 592 ENDDO 577 593 … … 584 600 k = nzb_diff_s_inner(j,i)-1 585 601 586 IF ( .NOT. cloud_physics )THEN602 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 587 603 k1 = 1.0 + 0.61 * q(k,j,i) 588 604 k2 = 0.61 * pt(k,j,i) 589 ELSE 605 ELSE IF ( cloud_physics ) THEN 590 606 IF ( ql(k,j,i) == 0.0 ) THEN 591 607 k1 = 1.0 + 0.61 * q(k,j,i) … … 601 617 k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) 602 618 ENDIF 619 ELSE IF ( cloud_droplets ) THEN 620 k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i) 621 k2 = 0.61 * pt(k,j,i) 603 622 ENDIF 604 623 … … 615 634 k = nzt 616 635 617 IF ( .NOT. cloud_physics )THEN636 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 618 637 k1 = 1.0 + 0.61 * q(k,j,i) 619 638 k2 = 0.61 * pt(k,j,i) 620 ELSE 639 ELSE IF ( cloud_physics ) THEN 621 640 IF ( ql(k,j,i) == 0.0 ) THEN 622 641 k1 = 1.0 + 0.61 * q(k,j,i) … … 632 651 k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) 633 652 ENDIF 653 ELSE IF ( cloud_droplets ) THEN 654 k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i) 655 k2 = 0.61 * pt(k,j,i) 634 656 ENDIF 635 657 … … 699 721 700 722 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 701 723 702 724 ENDDO 703 725 … … 721 743 722 744 IF ( wall_e_y(j,i) /= 0.0 ) THEN 723 ! 745 ! 724 746 !-- Inconsistency removed: as the thermal stratification 725 747 !-- is not taken into account for the evaluation of the … … 751 773 752 774 IF ( wall_e_x(j,i) /= 0.0 ) THEN 753 ! 775 ! 754 776 !-- Inconsistency removed: as the thermal stratification 755 777 !-- is not taken into account for the evaluation of the … … 805 827 806 828 IF ( wall_e_y(j,i) /= 0.0 ) THEN 807 ! 829 ! 808 830 !-- Inconsistency removed: as the thermal stratification 809 831 !-- is not taken into account for the evaluation of the … … 831 853 832 854 IF ( wall_e_x(j,i) /= 0.0 ) THEN 833 ! 855 ! 834 856 !-- Inconsistency removed: as the thermal stratification 835 857 !-- is not taken into account for the evaluation of the … … 1041 1063 DO k = nzb_diff_s_inner(j,i), nzt_diff 1042 1064 1043 IF ( .NOT. cloud_physics ) THEN1065 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 1044 1066 k1 = 1.0 + 0.61 * q(k,j,i) 1045 1067 k2 = 0.61 * pt(k,j,i) 1046 ELSE 1068 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & 1069 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & 1070 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & 1071 ) * dd2zu(k) 1072 ELSE IF ( cloud_physics ) THEN 1047 1073 IF ( ql(k,j,i) == 0.0 ) THEN 1048 1074 k1 = 1.0 + 0.61 * q(k,j,i) … … 1058 1084 k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) 1059 1085 ENDIF 1060 ENDIF 1061 1062 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & 1086 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & 1063 1087 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & 1064 1088 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & 1065 1089 ) * dd2zu(k) 1090 ELSE IF ( cloud_droplets ) THEN 1091 k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i) 1092 k2 = 0.61 * pt(k,j,i) 1093 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & 1094 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & 1095 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - & 1096 pt(k,j,i) * ( ql(k+1,j,i) - & 1097 ql(k-1,j,i) ) ) * dd2zu(k) 1098 ENDIF 1066 1099 ENDDO 1067 1100 … … 1069 1102 k = nzb_diff_s_inner(j,i)-1 1070 1103 1071 IF ( .NOT. cloud_physics )THEN1104 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 1072 1105 k1 = 1.0 + 0.61 * q(k,j,i) 1073 1106 k2 = 0.61 * pt(k,j,i) 1074 ELSE 1107 ELSE IF ( cloud_physics ) THEN 1075 1108 IF ( ql(k,j,i) == 0.0 ) THEN 1076 1109 k1 = 1.0 + 0.61 * q(k,j,i) … … 1086 1119 k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) 1087 1120 ENDIF 1121 ELSE IF ( cloud_droplets ) THEN 1122 k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i) 1123 k2 = 0.61 * pt(k,j,i) 1088 1124 ENDIF 1089 1125 … … 1095 1131 k = nzt 1096 1132 1097 IF ( .NOT. cloud_physics )THEN1133 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 1098 1134 k1 = 1.0 + 0.61 * q(k,j,i) 1099 1135 k2 = 0.61 * pt(k,j,i) 1100 ELSE 1136 ELSE IF ( cloud_physics ) THEN 1101 1137 IF ( ql(k,j,i) == 0.0 ) THEN 1102 1138 k1 = 1.0 + 0.61 * q(k,j,i) … … 1112 1148 k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) 1113 1149 ENDIF 1150 ELSE IF ( cloud_droplets ) THEN 1151 k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i) 1152 k2 = 0.61 * pt(k,j,i) 1114 1153 ENDIF 1115 1154 -
palm/trunk/SOURCE/sum_up_3d_data.f90
r979 r1007 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! Bugfix in calculation of ql_vp 7 7 ! 8 8 ! Former revisions: … … 388 388 389 389 CASE ( 'ql_vp' ) 390 DO i = nxlg, nxrg 391 DO j = nysg, nyng 392 DO k = nzb, nzt+1 393 ql_vp_av(k,j,i) = ql_vp_av(k,j,i) + ql_vp(k,j,i) 390 DO i = nxl, nxr 391 DO j = nys, nyn 392 DO k = nzb, nzt+1 393 psi = prt_start_index(k,j,i) 394 DO n = psi, psi+prt_count(k,j,i)-1 395 ql_vp_av(k,j,i) = ql_vp_av(k,j,i) + & 396 particles(n)%weight_factor / & 397 prt_count(k,j,i) 398 ENDDO 394 399 ENDDO 395 400 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
