Ignore:
Timestamp:
Jun 9, 2016 4:25:25 PM (8 years ago)
Author:
suehring
Message:

several bugfixes in particle model and serial mode

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/lpm_init.f90

    r1891 r1929  
    1919! Current revisions:
    2020! -----------------
    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.
    2227!
    2328! Former revisions:
     
    3136!
    3237! 1873 2016-04-18 14:50:06Z maronga
    33 ! Module renamed (removed _mod)
     38! Module renamed (removed _mod
     39
    3440!
    3541! 1871 2016-04-15 11:46:09Z hoffmann
     
    162168                prt_count, psb, psl, psn, psr, pss, pst,                       &
    163169                radius, random_start_position, read_particles_from_restartfile,&
    164                 seed_follows_topography, sort_count,                           &
     170                seed_follows_topography, sgs_wf_part, sort_count,              &
    165171                total_number_of_particles,                                     &
    166172                use_sgs_for_particles,                                         &
     
    190196    PUBLIC lpm_init, lpm_create_particle
    191197
    192 CONTAINS
     198 CONTAINS
    193199
    194200!------------------------------------------------------------------------------!
     
    290296
    291297!
    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.
    293300    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
    294301       ALLOCATE( de_dx(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
    295302                 de_dy(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
    296303                 de_dz(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     304
     305       sgs_wf_part = 1.0_wp / 3.0_wp   
    297306    ENDIF
    298307
     
    334343!
    335344!--    Precalculate LOG(z/z0)
    336        height_p    = 0.0_wp
     345       height_p    = z0_av_global
    337346       DO  k = 1, number_of_sublayers
    338347
     
    341350
    342351       ENDDO
    343 
    344352
    345353    ENDIF
     
    444452                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
    445453                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0, 0, 0, &
    446                                       0, .FALSE., -1)
     454                                      0, .FALSE., -1 )
    447455
    448456       particle_groups = particle_groups_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
     
    477485
    478486       CALL lpm_create_particle (PHASE_INIT)
    479 
    480487!
    481488!--    User modification of initial particles
     
    499506    number_of_particles = prt_count(nzb+1,nys,nxl)
    500507    particles => grid_particles(nzb+1,nys,nxl)%particles(1:number_of_particles)
    501 
    502508!
    503509!-- Formats
     
    520526
    521527    USE particle_attributes,                                                   &
    522         ONLY: monodisperse_aerosols
     528        ONLY: deleted_particles, monodisperse_aerosols
    523529
    524530    IMPLICIT  NONE
    525531
    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
    548555
    549556!
     
    618625                            tmp_particle%tail_id       = 0     !unused
    619626
     627
    620628!
    621629!--                         Determine the grid indices of the particle position
     
    628636!--                            Particle height is given relative to topography
    629637                               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))
    631640                               IF ( kp > nzt )  THEN
    632641                                  pos_x = pos_x + pdx(i)
    633642                                  CYCLE xloop
    634643                               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                               
    635648                            ENDIF
    636649
     
    644657                               ENDIF
    645658                               grid_particles(kp,jp,ip)%particles(local_count(kp,jp,ip)) = tmp_particle
     659
    646660                            ENDIF
    647661                         ENDDO
     
    691705                      ENDIF
    692706                   ENDIF
     707
    693708                ENDDO
    694709             ENDDO
    695710          ENDDO
    696711       ENDIF
     712
    697713    ENDDO
    698714
     
    707723
    708724!
    709 !-- Add random fluctuation to particle positions
     725!-- Add random fluctuation to particle positions.
    710726    IF ( random_start_position )  THEN
    711727       DO  ip = nxl, nxr
     
    715731                IF ( number_of_particles <= 0 )  CYCLE
    716732                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
    719740                   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)
    720743                      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                                   )
    723747                   ENDIF
    724748                   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)
    725751                      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                                   )
    728755                   ENDIF
    729756                   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)
    730759                      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                                   )
    733763                   ENDIF
    734764                ENDDO
    735765!
    736 !--             Identify particles located outside the model domain
     766!--             Identify particles located outside the model domain and reflect
     767!--             or absorb them if necessary.
    737768                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
    738783             ENDDO
    739784          ENDDO
Note: See TracChangeset for help on using the changeset viewer.