Changeset 1929 for palm/trunk/SOURCE/lpm_init.f90
- Timestamp:
- Jun 9, 2016 4:25:25 PM (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lpm_init.f90
r1891 r1929 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Bugfix in determining initial particle height and grid index in case of 22 ! seed_follows_topography. 23 ! Bugfix concerning random positions, ensure that particles do not move more 24 ! than one grid length. 25 ! Bugfix logarithmic interpolation. 26 ! Initial setting of sgs_wf_part. 22 27 ! 23 28 ! Former revisions: … … 31 36 ! 32 37 ! 1873 2016-04-18 14:50:06Z maronga 33 ! Module renamed (removed _mod) 38 ! Module renamed (removed _mod 39 34 40 ! 35 41 ! 1871 2016-04-15 11:46:09Z hoffmann … … 162 168 prt_count, psb, psl, psn, psr, pss, pst, & 163 169 radius, random_start_position, read_particles_from_restartfile,& 164 seed_follows_topography, s ort_count,&170 seed_follows_topography, sgs_wf_part, sort_count, & 165 171 total_number_of_particles, & 166 172 use_sgs_for_particles, & … … 190 196 PUBLIC lpm_init, lpm_create_particle 191 197 192 CONTAINS198 CONTAINS 193 199 194 200 !------------------------------------------------------------------------------! … … 290 296 291 297 ! 292 !-- Allocate arrays required for calculating particle SGS velocities 298 !-- Allocate arrays required for calculating particle SGS velocities. 299 !-- Initialize prefactor required for stoachastic Weil equation. 293 300 IF ( use_sgs_for_particles .AND. .NOT. cloud_droplets ) THEN 294 301 ALLOCATE( de_dx(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 295 302 de_dy(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 296 303 de_dz(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 304 305 sgs_wf_part = 1.0_wp / 3.0_wp 297 306 ENDIF 298 307 … … 334 343 ! 335 344 !-- Precalculate LOG(z/z0) 336 height_p = 0.0_wp345 height_p = z0_av_global 337 346 DO k = 1, number_of_sublayers 338 347 … … 341 350 342 351 ENDDO 343 344 352 345 353 ENDIF … … 444 452 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, & 445 453 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0, 0, 0, & 446 0, .FALSE., -1 )454 0, .FALSE., -1 ) 447 455 448 456 particle_groups = particle_groups_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp ) … … 477 485 478 486 CALL lpm_create_particle (PHASE_INIT) 479 480 487 ! 481 488 !-- User modification of initial particles … … 499 506 number_of_particles = prt_count(nzb+1,nys,nxl) 500 507 particles => grid_particles(nzb+1,nys,nxl)%particles(1:number_of_particles) 501 502 508 ! 503 509 !-- Formats … … 520 526 521 527 USE particle_attributes, & 522 ONLY: monodisperse_aerosols528 ONLY: deleted_particles, monodisperse_aerosols 523 529 524 530 IMPLICIT NONE 525 531 526 INTEGER(iwp) :: alloc_size !< 527 INTEGER(iwp) :: i !< 528 INTEGER(iwp) :: ip !< 529 INTEGER(iwp) :: j !< 530 INTEGER(iwp) :: jp !< 531 INTEGER(iwp) :: kp !< 532 INTEGER(iwp) :: loop_stride !< 533 INTEGER(iwp) :: n !< 534 INTEGER(iwp) :: new_size !< 535 536 INTEGER(iwp), INTENT(IN) :: phase !< 537 538 INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: local_count !< 539 INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: local_start !< 540 541 LOGICAL :: first_stride !< 542 543 REAL(wp) :: pos_x !< 544 REAL(wp) :: pos_y !< 545 REAL(wp) :: pos_z !< 546 547 TYPE(particle_type),TARGET :: tmp_particle !< 532 INTEGER(iwp) :: alloc_size !< relative increase of allocated memory for particles 533 INTEGER(iwp) :: i !< loop variable ( particle groups ) 534 INTEGER(iwp) :: ip !< index variable along x 535 INTEGER(iwp) :: j !< loop variable ( particles per point ) 536 INTEGER(iwp) :: jp !< index variable along y 537 INTEGER(iwp) :: kp !< index variable along z 538 INTEGER(iwp) :: loop_stride !< loop variable for initialization 539 INTEGER(iwp) :: n !< loop variable ( number of particles ) 540 INTEGER(iwp) :: new_size !< new size of allocated memory for particles 541 542 INTEGER(iwp), INTENT(IN) :: phase !< mode of inititialization 543 544 INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: local_count !< start address of new particle 545 INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: local_start !< start address of new particle 546 547 LOGICAL :: first_stride !< flag for initialization 548 549 REAL(wp) :: pos_x !< increment for particle position in x 550 REAL(wp) :: pos_y !< increment for particle position in y 551 REAL(wp) :: pos_z !< increment for particle position in z 552 REAL(wp) :: rand_contr !< dummy argument for random position 553 554 TYPE(particle_type),TARGET :: tmp_particle !< temporary particle used for initialization 548 555 549 556 ! … … 618 625 tmp_particle%tail_id = 0 !unused 619 626 627 620 628 ! 621 629 !-- Determine the grid indices of the particle position … … 628 636 !-- Particle height is given relative to topography 629 637 kp = kp + nzb_w_inner(jp,ip) 630 tmp_particle%z = tmp_particle%z + zw(kp) 638 tmp_particle%z = tmp_particle%z + & 639 zw(nzb_w_inner(jp,ip)) 631 640 IF ( kp > nzt ) THEN 632 641 pos_x = pos_x + pdx(i) 633 642 CYCLE xloop 634 643 ENDIF 644 ELSEIF ( .NOT. seed_follows_topography .AND. & 645 tmp_particle%z <= zw(nzb_w_inner(jp,ip)) ) THEN 646 pos_x = pos_x + pdx(i) 647 CYCLE xloop 635 648 ENDIF 636 649 … … 644 657 ENDIF 645 658 grid_particles(kp,jp,ip)%particles(local_count(kp,jp,ip)) = tmp_particle 659 646 660 ENDIF 647 661 ENDDO … … 691 705 ENDIF 692 706 ENDIF 707 693 708 ENDDO 694 709 ENDDO 695 710 ENDDO 696 711 ENDIF 712 697 713 ENDDO 698 714 … … 707 723 708 724 ! 709 !-- Add random fluctuation to particle positions 725 !-- Add random fluctuation to particle positions. 710 726 IF ( random_start_position ) THEN 711 727 DO ip = nxl, nxr … … 715 731 IF ( number_of_particles <= 0 ) CYCLE 716 732 particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles) 717 718 DO n = local_start(kp,jp,ip), number_of_particles !Move only new particles 733 ! 734 !-- Move only new particles. Moreover, limit random fluctuation 735 !-- in order to prevent that particles move more than one grid box, 736 !-- which would lead to problems concerning particle exchange 737 !-- between processors in case pdx/pdy are larger than dx/dy, 738 !-- respectively. 739 DO n = local_start(kp,jp,ip), number_of_particles 719 740 IF ( psl(particles(n)%group) /= psr(particles(n)%group) ) THEN 741 rand_contr = ( random_function( iran_part ) - 0.5_wp ) * & 742 pdx(particles(n)%group) 720 743 particles(n)%x = particles(n)%x + & 721 ( random_function( iran_part ) - 0.5_wp ) * & 722 pdx(particles(n)%group) 744 MERGE( rand_contr, SIGN( dx, rand_contr ), & 745 ABS( rand_contr ) < dx & 746 ) 723 747 ENDIF 724 748 IF ( pss(particles(n)%group) /= psn(particles(n)%group) ) THEN 749 rand_contr = ( random_function( iran_part ) - 0.5_wp ) * & 750 pdy(particles(n)%group) 725 751 particles(n)%y = particles(n)%y + & 726 ( random_function( iran_part ) - 0.5_wp ) * & 727 pdy(particles(n)%group) 752 MERGE( rand_contr, SIGN( dy, rand_contr ), & 753 ABS( rand_contr ) < dy & 754 ) 728 755 ENDIF 729 756 IF ( psb(particles(n)%group) /= pst(particles(n)%group) ) THEN 757 rand_contr = ( random_function( iran_part ) - 0.5_wp ) * & 758 pdz(particles(n)%group) 730 759 particles(n)%z = particles(n)%z + & 731 ( random_function( iran_part ) - 0.5_wp ) * & 732 pdz(particles(n)%group) 760 MERGE( rand_contr, SIGN( dz, rand_contr ), & 761 ABS( rand_contr ) < dz & 762 ) 733 763 ENDIF 734 764 ENDDO 735 765 ! 736 !-- Identify particles located outside the model domain 766 !-- Identify particles located outside the model domain and reflect 767 !-- or absorb them if necessary. 737 768 CALL lpm_boundary_conds( 'bottom/top' ) 769 ! 770 !-- Furthermore, remove particles located in topography. Note, as 771 !-- the particle speed is still zero at this point, wall 772 !-- reflection boundary conditions will not work in this case. 773 particles => & 774 grid_particles(kp,jp,ip)%particles(1:number_of_particles) 775 DO n = local_start(kp,jp,ip), number_of_particles 776 i = ( particles(n)%x + 0.5_wp * dx ) * ddx 777 j = ( particles(n)%y + 0.5_wp * dy ) * ddy 778 IF ( particles(n)%z <= zw(nzb_w_inner(j,i)) ) THEN 779 particles(n)%particle_mask = .FALSE. 780 deleted_particles = deleted_particles + 1 781 ENDIF 782 ENDDO 738 783 ENDDO 739 784 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.