Changeset 420 for palm/trunk
- Timestamp:
- Jan 13, 2010 3:10:53 PM (15 years ago)
- Location:
- palm/trunk
- Files:
-
- 4 edited
-
DOC/app/chapter_3.8.html (modified) (1 diff)
-
SOURCE/CURRENT_MODIFICATIONS (modified) (4 diffs)
-
SOURCE/advec_particles.f90 (modified) (45 diffs)
-
SOURCE/disturb_field.f90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/DOC/app/chapter_3.8.html
r349 r420 128 128 <P STYLE="line-height: 100%">It's also possible to perform precursor 129 129 runs (one atmospheric and one oceanic) followed by a coupled restart 130 run. In order to achieve this, the parameter <A HREF="chapter_4.1.html#coupling_start_time">coupling_ time_start</A>130 run. In order to achieve this, the parameter <A HREF="chapter_4.1.html#coupling_start_time">coupling_start_time</A> 131 131 must be set according to the <A HREF="../misc/precursor_run_control.pdf">documentation.</A></P> 132 132 <HR> -
palm/trunk/SOURCE/CURRENT_MODIFICATIONS
r412 r420 1 1 New: 2 2 --- 3 humidity=.T. is now usable for runs with topography. wall_humidityflux and 3 humidity=.T. is now usable for runs with topography. wall_humidityflux and 4 4 wall_scalarflux are the corresponding new parin arrays. 5 5 (check_parameters, init_3d_model, parin) 6 6 7 Large scale vertical motion (subsidence/ascent) can be applied to the 8 prognostic equation for the potential temperature. (check_parameters, header, 7 Large scale vertical motion (subsidence/ascent) can be applied to the 8 prognostic equation for the potential temperature. (check_parameters, header, 9 9 Makefile, modules, parin, prognostic_equations, read_var_list, subsidence, 10 10 write_var_list) … … 18 18 var_ts is replaced by dots_max (modules,init_3d_model) 19 19 20 init_3d_model, modules 20 Every cloud droplet has now an own weighting factor and can be deleted due to 21 collisions. Condensation and collision of cloud droplets are adjusted 22 accordingly. (advec_particles) 23 24 Collision efficiency for large cloud droplets has changed according to table of 25 Rogers and Yau. (collision_efficiency) 26 27 advec_particles, collision_efficiency, init_3d_model, modules 21 28 22 29 … … 25 32 Dimension of array stat in cascade change to prevent type problems with 26 33 mpi2 libraries (poisfft_hybrid) 34 35 A loop has been splitted to make runs reproducible on HLRN 36 systems. (disturb_field) 27 37 28 38 Bugfix: exchange of ghost points for prho included (time_integration) … … 35 45 messages from gfortran compiler (modules) 36 46 37 calc_precipitation, modules, poisfft_hybrid, sum_up_3d_data, time_integration 47 Bugfix: calculation of cloud droplet velocity (advec_particles) 48 49 Bugfix: transfer of particles at south/left edge (advec_particles) 50 51 Bugfix: calculation of collision_efficiency (collision_efficiency) 52 53 advec_particles, calc_precipitation, collision_efficiency, disturb_field, modules, poisfft_hybrid, sum_up_3d_data, time_integration 38 54 39 55 -
palm/trunk/SOURCE/advec_particles.f90
r392 r420 5 5 ! ----------------- 6 6 ! TEST: PRINT statements on unit 9 (commented out) 7 ! Own weighting factor for every cloud droplet is implemented and condensation 8 ! and collision of cloud droplets are adjusted accordingly. +delta_v, -s_r3, 9 ! -s_r4 10 ! Initialization of variables for the (sub-) timestep is moved to the beginning 11 ! in order to enable deletion of cloud droplets due to collision processes. 12 ! Collision efficiency for large cloud droplets has changed according to table 13 ! of Rogers and Yau. (collision_efficiency) 14 ! Bugfix: calculation of cloud droplet velocity 15 ! Bugfix: transfer of particles at south/left edge when new position 16 ! y>=(ny+0.5)*dy-1.e-12 or x>=(nx+0.5)*dx-1.e-12, very rare 17 ! Bugfix: calculation of y (collision_efficiency) 7 18 ! 8 19 ! Former revisions: … … 117 128 LOGICAL :: dt_3d_reached, dt_3d_reached_l, prt_position 118 129 119 REAL :: aa, arg, bb, cc, dd, delta_r, de ns_ratio, de_dt, de_dt_min,&120 de_d x_int, de_dx_int_l, de_dx_int_u, de_dy_int, de_dy_int_l,&121 de_dy_int_ u, de_dz_int, de_dz_int_l, de_dz_int_u, diss_int,&122 diss_int _l, diss_int_u, distance, dt_gap, dt_particle,&123 dt_particle _m, d_radius, d_sum, e_a, e_int, e_int_l, e_int_u,&124 e_ mean_int, e_s, exp_arg, exp_term, fs_int, gg,&125 lagr_timescale, mean_r, new_r, p_int, pt_int, pt_int_l,&130 REAL :: aa, arg, bb, cc, dd, delta_r, delta_v, dens_ratio, de_dt, & 131 de_dt_min, de_dx_int, de_dx_int_l, de_dx_int_u, de_dy_int, & 132 de_dy_int_l, de_dy_int_u, de_dz_int, de_dz_int_l, de_dz_int_u, & 133 diss_int, diss_int_l, diss_int_u, distance, dt_gap, & 134 dt_particle, dt_particle_m, d_radius, d_sum, e_a, e_int, & 135 e_int_l, e_int_u, e_mean_int, e_s, exp_arg, exp_term, fs_int, & 136 gg,lagr_timescale, mean_r, new_r, p_int, pt_int, pt_int_l, & 126 137 pt_int_u, q_int, q_int_l, q_int_u, ql_int, ql_int_l, ql_int_u, & 127 random_gauss, sl_r3, sl_r4, s_r3, s_r4, t_int, u_int, u_int_l,&128 u_int_u, vv_int, v_int, v_int_l, v_int_u, w_int, w_int_l, &129 w_int_u,x, y138 random_gauss, sl_r3, sl_r4, t_int, u_int, u_int_l, u_int_u, & 139 vv_int, v_int, v_int_l, v_int_u, w_int, w_int_l, w_int_u, & 140 x, y 130 141 131 142 REAL, DIMENSION(1:30) :: de_dxi, de_dyi, de_dzi, dissi, d_gp_pl, ei … … 180 191 CALL cpu_log( log_point_s(40), 'advec_part_io', 'stop' ) 181 192 ENDIF 193 194 ! 195 !-- Initialize variables for the (sub-) timestep, i.e. for 196 !-- marking those particles to be deleted after the timestep 197 particle_mask = .TRUE. 198 deleted_particles = 0 199 trlp_count_recv = 0 200 trnp_count_recv = 0 201 trrp_count_recv = 0 202 trsp_count_recv = 0 203 trlpt_count_recv = 0 204 trnpt_count_recv = 0 205 trrpt_count_recv = 0 206 trspt_count_recv = 0 207 IF ( use_particle_tails ) THEN 208 tail_mask = .TRUE. 209 ENDIF 210 deleted_tails = 0 182 211 183 212 ! … … 215 244 !-- Reset summation arrays 216 245 ql_c = 0.0; ql_v = 0.0; ql_vp = 0.0 246 delta_r=0.0; delta_v=0.0 217 247 218 248 ! … … 343 373 !-- Sum up the total volume of liquid water (needed below for 344 374 !-- re-calculating the weighting factors) 345 ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * &375 ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * & 346 376 particles(n)%radius**3 347 377 ENDDO … … 351 381 !-- Particle growth by collision 352 382 CALL cpu_log( log_point_s(43), 'advec_part_coll', 'start' ) 353 383 354 384 DO i = nxl, nxr 355 385 DO j = nys, nyn … … 374 404 tmp_particle = particles(psi+is) 375 405 js = is 376 DO WHILE ( particles(psi+js-inc)%radius > &406 DO WHILE ( particles(psi+js-inc)%radius > & 377 407 tmp_particle%radius ) 378 408 particles(psi+js) = particles(psi+js-inc) … … 386 416 ! 387 417 !-- Calculate the mean radius of all those particles which 388 !-- are of smaller or equalsize than the current particle418 !-- are of smaller size than the current particle 389 419 !-- and use this radius for calculating the collision efficiency 390 420 psi = prt_start_index(k,j,i) 391 s_r3 = 0.0 392 s_r4 = 0.0 393 DO n = psi, psi+prt_count(k,j,i)-1 394 ! 395 !-- There may be some particles of size equal to the 396 !-- current particle but with larger index 421 422 DO n = psi+prt_count(k,j,i)-1, psi+1, -1 397 423 sl_r3 = 0.0 398 424 sl_r4 = 0.0 399 DO is = n, psi+prt_count(k,j,i)-2 400 IF ( particles(is+1)%radius == & 401 particles(is)%radius ) THEN 402 sl_r3 = sl_r3 + particles(is+1)%radius**3 403 sl_r4 = sl_r4 + particles(is+1)%radius**4 404 ELSE 405 EXIT 406 ENDIF 407 ENDDO 408 409 IF ( ( s_r3 + sl_r3 ) > 0.0 ) THEN 410 411 mean_r = ( s_r4 + sl_r4 ) / ( s_r3 + sl_r3 ) 412 413 CALL collision_efficiency( mean_r, & 414 particles(n)%radius, & 425 426 DO is = n-1, psi, -1 427 IF ( particles(is)%radius < particles(n)%radius ) THEN 428 sl_r3 = sl_r3 + particles(is)%weight_factor & 429 *( particles(is)%radius**3 ) 430 sl_r4 = sl_r4 + particles(is)%weight_factor & 431 *( particles(is)%radius**4 ) 432 ENDIF 433 ENDDO 434 435 IF ( ( sl_r3 ) > 0.0 ) THEN 436 mean_r = ( sl_r4 ) / ( sl_r3 ) 437 438 CALL collision_efficiency( mean_r, & 439 particles(n)%radius, & 415 440 effective_coll_efficiency ) 416 441 … … 418 443 effective_coll_efficiency = 0.0 419 444 ENDIF 420 421 ! 422 !-- Contribution of the current particle to the next one 423 s_r3 = s_r3 + particles(n)%radius**3 424 s_r4 = s_r4 + particles(n)%radius**4 425 426 IF ( effective_coll_efficiency > 1.0 .OR. & 427 effective_coll_efficiency < 0.0 ) & 445 446 IF ( effective_coll_efficiency > 1.0 .OR. & 447 effective_coll_efficiency < 0.0 ) & 428 448 THEN 429 449 WRITE( message_string, * ) 'collision_efficiency ' , & … … 431 451 CALL message( 'advec_particles', 'PA0145', 2, 2, -1, & 432 452 6, 1 ) 433 END 453 ENDIF 434 454 435 455 ! … … 459 479 ) / ( 3.0 * gg ) 460 480 461 ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz * &481 ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz * & 462 482 ( ql_int_u - ql_int_l ) 463 483 … … 478 498 gg = aa + bb + cc + dd 479 499 480 u_int_l = ( ( gg-aa ) * u(kk,jj,ii) + ( gg-bb ) * &481 u(kk,jj,ii+1) &482 + ( gg-cc ) * u(kk,jj+1,ii) + ( gg-dd ) * &483 u(kk,jj+1,ii+1) &500 u_int_l = ( ( gg-aa ) * u(kk,jj,ii) + ( gg-bb ) * & 501 u(kk,jj,ii+1) & 502 + ( gg-cc ) * u(kk,jj+1,ii) + ( gg-dd ) * & 503 u(kk,jj+1,ii+1) & 484 504 ) / ( 3.0 * gg ) - u_gtrans 485 505 IF ( kk+1 == nzt+1 ) THEN … … 491 511 u(kk+1,jj+1,ii+1) & 492 512 ) / ( 3.0 * gg ) - u_gtrans 493 u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz * &513 u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz * & 494 514 ( u_int_u - u_int_l ) 495 515 ENDIF … … 558 578 559 579 ! 560 !-- Change in radius due to collision 561 delta_r = effective_coll_efficiency * & 562 ql_int * rho_surface / ( 1.0 - ql_int ) * & 563 0.25 / rho_l * & 564 SQRT( ( u_int - particles(n)%speed_x )**2 + & 565 ( v_int - particles(n)%speed_y )**2 + & 566 ( w_int - particles(n)%speed_z )**2 & 567 ) * dt_3d 580 !-- Change in radius due to collision 581 delta_r = effective_coll_efficiency / 3.0 & 582 * pi * sl_r3 * ddx * ddy / dz & 583 * SQRT( ( u_int - particles(n)%speed_x )**2 & 584 + ( v_int - particles(n)%speed_y )**2 & 585 + ( w_int - particles(n)%speed_z )**2 & 586 ) * dt_3d 587 ! 588 !-- Change in volume due to collision 589 delta_v = particles(n)%weight_factor & 590 * ( ( particles(n)%radius + delta_r )**3 & 591 - particles(n)%radius**3 ) 592 593 ! 594 !-- Check if collected particles provide enough LWC for volume 595 !-- change of collector particle 596 IF ( delta_v >= sl_r3 .and. sl_r3 > 0.0 ) THEN 597 598 delta_r = ( ( sl_r3/particles(n)%weight_factor ) & 599 + particles(n)%radius**3 )**( 1./3. ) & 600 - particles(n)%radius 601 602 DO is = n-1, psi, -1 603 IF ( particles(is)%radius < particles(n)%radius ) & 604 THEN 605 particles(is)%weight_factor = 0.0 606 particle_mask(is) = .FALSE. 607 deleted_particles = deleted_particles + 1 608 ENDIF 609 ENDDO 610 611 ELSE IF ( delta_v < sl_r3 .and. sl_r3 > 0.0 ) THEN 612 613 DO is = n-1, psi, -1 614 IF ( particles(is)%radius < particles(n)%radius & 615 .and. sl_r3 > 0.0 ) THEN 616 particles(is)%weight_factor = & 617 ( ( particles(is)%weight_factor & 618 * ( particles(is)%radius**3 ) ) & 619 - ( delta_v & 620 * particles(is)%weight_factor & 621 * ( particles(is)%radius**3 ) & 622 / sl_r3 ) ) & 623 / ( particles(is)%radius**3 ) 624 625 IF (particles(is)%weight_factor < 0.0) THEN 626 WRITE( message_string, * ) 'negative ', & 627 'weighting factor: ', & 628 particles(is)%weight_factor 629 CALL message( 'advec_particles', '', 2, 2, & 630 -1, 6, 1 ) 631 ENDIF 632 ENDIF 633 ENDDO 634 ENDIF 568 635 569 636 particles(n)%radius = particles(n)%radius + delta_r 570 571 ql_vp(k,j,i) = ql_vp(k,j,i) + particles(n)%radius**3 572 637 ql_vp(k,j,i) = ql_vp(k,j,i) + particles(n)%weight_factor & 638 *( particles(n)%radius**3 ) 573 639 ENDDO 574 640 575 ENDIF 576 577 ! 578 !-- Re-calculate the weighting factor (total liquid water content 579 !-- must be conserved during collision) 580 IF ( ql_vp(k,j,i) /= 0.0 ) THEN 581 582 ql_vp(k,j,i) = ql_v(k,j,i) / ql_vp(k,j,i) 583 ! 584 !-- Re-assign this weighting factor to the particles of the 585 !-- current gridbox 641 ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor & 642 *( particles(psi)%radius**3 ) 643 644 ELSE IF ( prt_count(k,j,i) == 1 ) THEN 586 645 psi = prt_start_index(k,j,i) 587 DO n = psi, psi + prt_count(k,j,i)-1 588 particles(n)%weight_factor = ql_vp(k,j,i) 589 ENDDO 590 646 ql_vp(k,j,i) = particles(psi)%weight_factor & 647 *( particles(psi)%radius**3 ) 648 ENDIF 649 ! 650 !-- Check if condensation of LWC was conserved 651 !-- during collision process 652 IF (ql_v(k,j,i) /= 0.0 ) THEN 653 IF (ql_vp(k,j,i)/ql_v(k,j,i) >= 1.00001 .or. & 654 ql_vp(k,j,i)/ql_v(k,j,i) <= 0.99999 ) THEN 655 WRITE( message_string, * ) 'LWC is not conserved during',& 656 ' collision! ', & 657 'LWC after condensation: ', & 658 ql_v(k,j,i), & 659 ' LWC after collision: ', & 660 ql_vp(k,j,i) 661 CALL message( 'advec_particles', '', 2, 2, -1, 6, 1 ) 662 ENDIF 591 663 ENDIF 592 664 … … 726 798 ! 727 799 !-- Final values are obtained by division by the total number of grid 728 !-- points used for the summation. 800 !-- points used for the summation. 729 801 hom(:,1,1,0) = sums(:,1) / ngp_2dh_outer(:,0) ! u 730 802 hom(:,1,2,0) = sums(:,2) / ngp_2dh_outer(:,0) ! v … … 761 833 CALL MPI_ALLREDUCE( sums_l(nzb,32,0), sums(nzb,32), nzt+2-nzb, & 762 834 MPI_REAL, MPI_SUM, comm2d, ierr ) 763 835 764 836 #else 765 837 sums(:,8) = sums_l(:,8,0) … … 779 851 ENDIF 780 852 781 ENDIF 853 ENDIF 782 854 783 855 ! … … 815 887 dt_3d_reached_l = .TRUE. 816 888 817 !818 !-- Initialize variables for the (sub-) timestep, i.e. for marking those819 !-- particles to be deleted after the timestep820 particle_mask = .TRUE.821 deleted_particles = 0822 trlp_count_recv = 0823 trnp_count_recv = 0824 trrp_count_recv = 0825 trsp_count_recv = 0826 trlpt_count_recv = 0827 trnpt_count_recv = 0828 trrpt_count_recv = 0829 trspt_count_recv = 0830 IF ( use_particle_tails ) THEN831 tail_mask = .TRUE.832 ENDIF833 deleted_tails = 0834 835 836 889 DO n = 1, number_of_particles 837 890 ! … … 847 900 k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz & 848 901 + offset_ocean_nzt ! only exact if equidistant 849 902 850 903 ! 851 904 !-- Interpolation of the velocity components in the xy-plane … … 942 995 + offset_ocean_nzt ! only exact if eq.dist 943 996 944 IF ( topography == 'flat' ) THEN 997 IF ( topography == 'flat' ) THEN 945 998 946 999 x = particles(n)%x - i * dx … … 1066 1119 !-- gp_outside_of_building(7): i+1,j,k+1 1067 1120 !-- gp_outside_of_building(8): i+1,j+1,k+1 1068 1121 1069 1122 gp_outside_of_building = 0 1070 1123 location = 0.0 … … 1092 1145 location(num_gp,3) = k * dz - 0.5 * dz 1093 1146 ei(num_gp) = e(k,j+1,i) 1094 dissi(num_gp) = diss(k,j+1,i) 1147 dissi(num_gp) = diss(k,j+1,i) 1095 1148 de_dxi(num_gp) = de_dx(k,j+1,i) 1096 1149 de_dyi(num_gp) = de_dy(k,j+1,i) … … 1105 1158 location(num_gp,3) = (k+1) * dz - 0.5 * dz 1106 1159 ei(num_gp) = e(k+1,j,i) 1107 dissi(num_gp) = diss(k+1,j,i) 1160 dissi(num_gp) = diss(k+1,j,i) 1108 1161 de_dxi(num_gp) = de_dx(k+1,j,i) 1109 1162 de_dyi(num_gp) = de_dy(k+1,j,i) … … 1119 1172 location(num_gp,3) = (k+1) * dz - 0.5 * dz 1120 1173 ei(num_gp) = e(k+1,j+1,i) 1121 dissi(num_gp) = diss(k+1,j+1,i) 1174 dissi(num_gp) = diss(k+1,j+1,i) 1122 1175 de_dxi(num_gp) = de_dx(k+1,j+1,i) 1123 1176 de_dyi(num_gp) = de_dy(k+1,j+1,i) … … 1133 1186 location(num_gp,3) = k * dz - 0.5 * dz 1134 1187 ei(num_gp) = e(k,j,i+1) 1135 dissi(num_gp) = diss(k,j,i+1) 1188 dissi(num_gp) = diss(k,j,i+1) 1136 1189 de_dxi(num_gp) = de_dx(k,j,i+1) 1137 1190 de_dyi(num_gp) = de_dy(k,j,i+1) … … 1147 1200 location(num_gp,3) = k * dz - 0.5 * dz 1148 1201 ei(num_gp) = e(k,j+1,i+1) 1149 dissi(num_gp) = diss(k,j+1,i+1) 1202 dissi(num_gp) = diss(k,j+1,i+1) 1150 1203 de_dxi(num_gp) = de_dx(k,j+1,i+1) 1151 1204 de_dyi(num_gp) = de_dy(k,j+1,i+1) … … 1161 1214 location(num_gp,3) = (k+1) * dz - 0.5 * dz 1162 1215 ei(num_gp) = e(k+1,j,i+1) 1163 dissi(num_gp) = diss(k+1,j,i+1) 1216 dissi(num_gp) = diss(k+1,j,i+1) 1164 1217 de_dxi(num_gp) = de_dx(k+1,j,i+1) 1165 1218 de_dyi(num_gp) = de_dy(k+1,j,i+1) … … 1175 1228 location(num_gp,3) = (k+1) * dz - 0.5 * dz 1176 1229 ei(num_gp) = e(k+1,j+1,i+1) 1177 dissi(num_gp) = diss(k+1,j+1,i+1) 1230 dissi(num_gp) = diss(k+1,j+1,i+1) 1178 1231 de_dxi(num_gp) = de_dx(k+1,j+1,i+1) 1179 1232 de_dyi(num_gp) = de_dy(k+1,j+1,i+1) … … 1182 1235 1183 1236 ! 1184 !-- If all gridpoints are situated outside of a building, then the 1237 !-- If all gridpoints are situated outside of a building, then the 1185 1238 !-- ordinary interpolation scheme can be used. 1186 1239 IF ( num_gp == 8 ) THEN … … 1309 1362 de_dxi(num_gp) = 0.0 1310 1363 de_dyi(num_gp) = de_dy(k,j,i) 1311 de_dzi(num_gp) = de_dz(k,j,i) 1312 ENDIF 1313 1364 de_dzi(num_gp) = de_dz(k,j,i) 1365 ENDIF 1366 1314 1367 IF ( gp_outside_of_building(5) == 1 .AND. & 1315 1368 gp_outside_of_building(1) == 0 ) THEN … … 1367 1420 de_dxi(num_gp) = 0.0 1368 1421 de_dyi(num_gp) = de_dy(k,j+1,i) 1369 de_dzi(num_gp) = de_dz(k,j+1,i) 1370 ENDIF 1371 1422 de_dzi(num_gp) = de_dz(k,j+1,i) 1423 ENDIF 1424 1372 1425 IF ( gp_outside_of_building(6) == 1 .AND. & 1373 1426 gp_outside_of_building(2) == 0 ) THEN … … 1425 1478 de_dxi(num_gp) = 0.0 1426 1479 de_dyi(num_gp) = de_dy(k+1,j,i) 1427 de_dzi(num_gp) = de_dz(k+1,j,i) 1428 ENDIF 1429 1480 de_dzi(num_gp) = de_dz(k+1,j,i) 1481 ENDIF 1482 1430 1483 IF ( gp_outside_of_building(7) == 1 .AND. & 1431 1484 gp_outside_of_building(3) == 0 ) THEN … … 1483 1536 de_dxi(num_gp) = 0.0 1484 1537 de_dyi(num_gp) = de_dy(k+1,j+1,i) 1485 de_dzi(num_gp) = de_dz(k+1,j+1,i) 1486 ENDIF 1487 1538 de_dzi(num_gp) = de_dz(k+1,j+1,i) 1539 ENDIF 1540 1488 1541 IF ( gp_outside_of_building(8) == 1 .AND. & 1489 1542 gp_outside_of_building(4) == 0 ) THEN … … 1600 1653 IF ( num_gp == 1 ) THEN 1601 1654 ! 1602 !-- If only one of the gridpoints is situated outside of the 1603 !-- building, it follows that the values at the particle 1655 !-- If only one of the gridpoints is situated outside of the 1656 !-- building, it follows that the values at the particle 1604 1657 !-- location are the same as the gridpoint values 1605 1658 e_int = ei(num_gp) … … 1612 1665 d_sum = 0.0 1613 1666 ! 1614 !-- Evaluation of the distances between the gridpoints 1615 !-- contributing to the interpolated values, and the particle 1667 !-- Evaluation of the distances between the gridpoints 1668 !-- contributing to the interpolated values, and the particle 1616 1669 !-- location 1617 1670 DO agp = 1, num_gp … … 1839 1892 IF ( cloud_droplets ) THEN 1840 1893 exp_arg = 4.5 * dens_ratio * molecular_viscosity / & 1841 ( particles(n)%radius )**2 /&1894 ( particles(n)%radius )**2 * & 1842 1895 ( 1.0 + 0.15 * ( 2.0 * particles(n)%radius * & 1843 1896 SQRT( ( u_int - particles(n)%speed_x )**2 + & … … 2178 2231 deleted_particles = deleted_particles + 1 2179 2232 2233 IF ( trlp(trlp_count)%x >= (nx + 0.5)* dx - 1.e-12 ) THEN 2234 trlp(trlp_count)%x = trlp(trlp_count)%x - 1.e-10 2235 trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x & 2236 - 1 2237 ENDIF 2238 2180 2239 IF ( use_particle_tails .AND. nn /= 0 ) THEN 2181 2240 trlpt_count = trlpt_count + 1 … … 2617 2676 deleted_particles = deleted_particles + 1 2618 2677 2678 IF ( trsp(trsp_count)%y >= (ny+0.5)* dy - 1.e-12 ) THEN 2679 trsp(trsp_count)%y = trsp(trsp_count)%y - 1.e-10 2680 trsp(trsp_count)%origin_y = & 2681 trsp(trsp_count)%origin_y - 1 2682 ENDIF 2683 2619 2684 IF ( use_particle_tails .AND. nn /= 0 ) THEN 2620 2685 trspt_count = trspt_count + 1 … … 3296 3361 IF ( dt_3d_reached ) EXIT 3297 3362 3363 ! 3364 !-- Initialize variables for the next (sub-) timestep, i.e. for marking those 3365 !-- particles to be deleted after the timestep 3366 particle_mask = .TRUE. 3367 deleted_particles = 0 3368 trlp_count_recv = 0 3369 trnp_count_recv = 0 3370 trrp_count_recv = 0 3371 trsp_count_recv = 0 3372 trlpt_count_recv = 0 3373 trnpt_count_recv = 0 3374 trrpt_count_recv = 0 3375 trspt_count_recv = 0 3376 IF ( use_particle_tails ) THEN 3377 tail_mask = .TRUE. 3378 ENDIF 3379 deleted_tails = 0 3380 3298 3381 ENDDO ! timestep loop 3299 3300 3382 3301 3383 ! … … 3308 3390 ENDIF 3309 3391 3310 3311 !3312 !-- Re-evaluate the weighting factors. After advection, particles within a3313 !-- grid box may have different weighting factors if some have been advected3314 !-- from a neighbouring box. The factors are re-evaluated so that they are3315 !-- the same for all particles of one box. This procedure must conserve the3316 !-- liquid water content within one box.3317 3392 IF ( cloud_droplets ) THEN 3318 3393 … … 3322 3397 3323 3398 ! 3324 !-- Re-calculate the weighting factors and calculate the liquid water content3399 !-- Calculate the liquid water content 3325 3400 DO i = nxl, nxr 3326 3401 DO j = nys, nyn … … 3328 3403 3329 3404 ! 3330 !-- Calculate the total volume of particles in the boxes (ql_vp) as 3331 !-- well as the real volume (ql_v, weighting factor has to be 3332 !-- included) 3405 !-- Calculate the total volume in the boxes (ql_v, weighting factor 3406 !-- has to beincluded) 3333 3407 psi = prt_start_index(k,j,i) 3334 3408 DO n = psi, psi+prt_count(k,j,i)-1 3335 ql_vp(k,j,i) = ql_vp(k,j,i) + particles(n)%radius**3 3336 3337 ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * & 3409 ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * & 3338 3410 particles(n)%radius**3 3339 3411 ENDDO 3340 3412 3341 3413 ! 3342 !-- Re-calculate the weighting factors and calculate the liquid 3343 !-- water content 3344 IF ( ql_vp(k,j,i) /= 0.0 ) THEN 3345 ql_vp(k,j,i) = ql_v(k,j,i) / ql_vp(k,j,i) 3346 ql(k,j,i) = ql(k,j,i) + rho_l * 1.33333333 * pi * & 3347 ql_v(k,j,i) / & 3414 !-- Calculate the liquid water content 3415 IF ( ql_v(k,j,i) /= 0.0 ) THEN 3416 ql(k,j,i) = ql(k,j,i) + rho_l * 1.33333333 * pi * & 3417 ql_v(k,j,i) / & 3348 3418 ( rho_surface * dx * dy * dz ) 3419 3420 IF ( ql(k,j,i) < 0.0 ) THEN 3421 WRITE( message_string, * ) 'LWC out of range: ' , & 3422 ql(k,j,i) 3423 CALL message( 'advec_particles', '', 2, 2, -1, 6, 1 ) 3424 ENDIF 3425 3349 3426 ELSE 3350 3427 ql(k,j,i) = 0.0 3351 3428 ENDIF 3352 3353 !3354 !-- Re-assign the weighting factor to the particles3355 DO n = psi, psi+prt_count(k,j,i)-13356 particles(n)%weight_factor = ql_vp(k,j,i)3357 ENDDO3358 3429 3359 3430 ENDDO … … 3833 3904 e = 0.001 3834 3905 ELSEIF ( mean_rm >= 25.0 ) THEN 3835 IF( j <= 3 ) e = 0.55 3906 IF( j <= 2 ) e = 0.0 3907 IF( j == 3 ) e = 0.47 3836 3908 IF( j == 4 ) e = 0.8 3837 3909 IF( j == 5 ) e = 0.9 3838 3910 IF( j >=6 ) e = 1.0 3839 3911 ELSEIF ( rm >= 3000.0 ) THEN 3840 e = 1.0 3912 IF( i == 1 ) e = 0.02 3913 IF( i == 2 ) e = 0.16 3914 IF( i == 3 ) e = 0.33 3915 IF( i == 4 ) e = 0.55 3916 IF( i == 5 ) e = 0.71 3917 IF( i == 6 ) e = 0.81 3918 IF( i == 7 ) e = 0.90 3919 IF( i >= 8 ) e = 0.94 3841 3920 ELSE 3842 3921 x = mean_rm - collected_r(i) 3843 y = rm - collect ed_r(j)3844 dx = collected_r(i+1) - collected_r(i) 3845 dy = collector_r(j+1) - collector_r(j) 3922 y = rm - collector_r(j) 3923 dx = collected_r(i+1) - collected_r(i) 3924 dy = collector_r(j+1) - collector_r(j) 3846 3925 aa = x**2 + y**2 3847 3926 bb = ( dx - x )**2 + y**2 … … 3854 3933 ENDIF 3855 3934 3856 END SUBROUTINE collision_efficiency 3935 END SUBROUTINE collision_efficiency 3857 3936 3858 3937 … … 3902 3981 WRITE( message_string, * ) ' particle out of range: i=', i, ' j=', & 3903 3982 j, ' k=', k, & 3904 ' &nxl=', nxl, ' nxr=', nxr,&3983 ' nxl=', nxl, ' nxr=', nxr, & 3905 3984 ' nys=', nys, ' nyn=', nyn, & 3906 3985 ' nzb=', nzb, ' nzt=', nzt -
palm/trunk/SOURCE/disturb_field.f90
r77 r420 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! A loop has been splitted to make runs reproducible on HLRN systems 7 7 ! 8 8 ! Former revisions: … … 98 98 ! 99 99 !-- Exchange of ghost points for the random perturbation 100 100 101 CALL exchange_horiz( dist1 ) 101 102 … … 104 105 !-- Neighboured grid points in all three directions are used for the 105 106 !-- filter operation. 106 DO i = nxl, nxr 107 DO j = nys, nyn 107 !-- Loop has been splitted to make runs reproducible on HLRN systems using 108 !-- compiler option -O3 109 DO i = nxl, nxr 110 DO j = nys, nyn 108 111 DO k = disturbance_level_ind_b-1, disturbance_level_ind_t+1 109 dist2(k,j,i) = ( dist1(k,j,i-1) + dist1(k,j,i+1) + dist1(k,j-1,i) & 110 + dist1(k,j+1,i) + dist1(k+1,j,i) + dist1(k-1,j,i) & 112 dist2(k,j,i) = ( dist1(k,j,i-1) + dist1(k,j,i+1) & 113 + dist1(k,j+1,i) + dist1(k+1,j,i) & 114 ) / 12.0 115 ENDDO 116 DO k = disturbance_level_ind_b-1, disturbance_level_ind_t+1 117 dist2(k,j,i) = dist2(k,j,i) + ( dist1(k,j-1,i) + dist1(k-1,j,i) & 111 118 + 6.0 * dist1(k,j,i) & 112 119 ) / 12.0 113 120 ENDDO 114 ENDDO115 ENDDO121 ENDDO 122 ENDDO 116 123 117 124 ! … … 119 126 !-- Afterwards, filter operation and exchange of ghost points are repeated. 120 127 CALL exchange_horiz( dist2 ) 128 121 129 DO i = nxl, nxr 122 130 DO j = nys, nyn … … 129 137 ENDDO 130 138 ENDDO 139 131 140 CALL exchange_horiz( dist1 ) 132 141
Note: See TracChangeset
for help on using the changeset viewer.
