! $Id$
!
! This module contains all the main structure needed for particles.
!
module Particles_main
!
use Cdata
use General, only: keep_compiler_quiet
use Messages
use Particles
use Particles_adaptation
use Particles_adsorbed
use Particles_cdata
use Particles_chemistry
use Particles_coagulation
use Particles_condensation
use Particles_collisions
use Particles_density
use Particles_diagnos_dv
use Particles_diagnos_state
use Particles_drag
use Particles_map
use Particles_mass
use Particles_mpicomm
use Particles_number
use Particles_radius
use Particles_grad
use Particles_selfgravity
use Particles_sink
use Particles_spin
use Particles_stalker
use Particles_stirring
use Particles_sub
use Particles_surfspec
use Particles_temperature
use Particles_lyapunov
use Particles_caustics
use Particles_tetrad
!
implicit none
!
include 'particles_main.h'
!
real, dimension(mpar_loc,mparray) :: fp = 0.0
real, dimension(mpar_loc,mpvar) :: dfp = 0.0
integer, dimension(mpar_loc,3) :: ineargrid = 0
!
contains
!***********************************************************************
subroutine particles_register_modules
!
! Register particle modules.
!
! 07-jan-05/anders: coded
!
use Special, only: register_particles_special
!
integer :: ipvar
!
call register_particles
call register_particles_lyapunov
call register_particles_radius
call register_particles_grad
call register_particles_spin
call register_particles_number
call register_particles_density
call register_particles_stirring
call register_particles_selfgrav
call register_particles_sink
call register_particles_TT
call register_particles_mass
call register_particles_drag
call register_particles_chem
call register_particles_ads
call register_particles_surfspec
call register_pars_diagnos_state
call register_particles_special (npvar)
!
! Print summary of variable names.
!
! [PAB]: We now have a "particle_index.pro"; can we remove the "pvarname.dat" file?
if (lroot) then
open(3,file=trim(datadir)//'/pvarname.dat',status='replace')
do ipvar=1,mparray
write(3,"(i4,2x,a)") ipvar, pvarname(ipvar)
enddo
close(3)
endif
!
endsubroutine particles_register_modules
!***********************************************************************
subroutine particles_rprint_list(lreset)
!
! Read names of diagnostic particle variables to print out during run.
!
! 07-jan-05/anders: coded
!
logical :: lreset
!
call rprint_particles (lreset,LWRITE=lroot)
call rprint_particles_radius (lreset,LWRITE=lroot)
call rprint_particles_grad (lreset,LWRITE=lroot)
call rprint_particles_lyapunov (lreset,LWRITE=lroot)
call rprint_particles_sink (lreset,LWRITE=lroot)
call rprint_particles_spin (lreset,LWRITE=lroot)
call rprint_particles_number (lreset,LWRITE=lroot)
call rprint_particles_density (lreset,LWRITE=lroot)
call rprint_particles_selfgrav (lreset,LWRITE=lroot)
call rprint_particles_TT (lreset,LWRITE=lroot)
call rprint_particles_mass (lreset,LWRITE=lroot)
call rprint_particles_ads (lreset,LWRITE=lroot)
call rprint_particles_surf (lreset,LWRITE=lroot)
call rprint_particles_chem (lreset,LWRITE=lroot)
call rprint_particles_coagulation (lreset,LWRITE=lroot)
call rprint_particles_condensation (lreset,LWRITE=lroot)
call rprint_particles_collisions (lreset,LWRITE=lroot)
call rprint_particles_diagnos_dv (lreset,LWRITE=lroot)
call rprint_particles_diagnos_state(lreset,LWRITE=lroot)
!
! check for those quantities for which we want video slices
!
if (lwrite_slices) then
where(cnamev=='rhop'.or.cnamev=='np'.or.cnamev=='vvp') cformv='DEFINED'
endif
!
endsubroutine particles_rprint_list
!***********************************************************************
subroutine particles_initialize_modules(f)
!
! Initialize particle modules.
!
! 07-jan-05/anders: coded
!
real, dimension (mx,my,mz,mfarray) :: f
!
if (lyinyang) &
call fatal_error('particles_initialize_modules','Particles not implemented on Yin-Yang grid')
!
! Check if there is enough total space allocated for particles.
!
if (npar/mpar_loc>ncpus) then
if (lroot) then
print*, 'particles_initialize_modules: '// &
'total number of particle slots available at the processors '// &
'is smaller than the number of particles!'
print*, 'particles_initialize_modules: npar/ncpus=', npar/ncpus
print*, 'particles_initialize_modules: mpar_loc-ncpus*npar_mig=', &
mpar_loc-ncpus*npar_mig
endif
call fatal_error('particles_initialize_modules','not enough space')
endif
!
! Set mass and number density of individual particles inside each
! superparticle. The mass of a superparticle is defined through
!
! * mpmat : mass of constituent particles
! * np_swarm : number density of constituent particles
! * rhop_swarm : mass density of superparticle
!
! One can either input these quantities by hand or set the wanted number
! density and mass density per grid cell (np_const or rhop_const).
!
if (rhop_const /= 0.0) then
rhop_swarm = rhop_const / (real(npar) / nwgrid)
if (all(lequidist).and.lcartesian_coords) then
mp_swarm = rhop_swarm * dVol(1)
else
call fatal_error('particles_initialize_modules',&
'dVol only defined for equidistant Cartesian coordinates')
endif
if (mpmat /= 0.0 .or. np_swarm /= 0.0) then
if (lparticles_radius) &
call fatal_error('particles_initialize_modules', 'may not set mpmat or np_swarm when setting rhop_const')
if (mpmat /= 0.0 .and. np_swarm /= 0.0) &
call fatal_error('particles_initialize_modules', 'must set only mpmat or np_swarm when using rhop_const')
if (mpmat == 0.0) mpmat = rhop_swarm / np_swarm
if (np_swarm == 0.0) np_swarm = rhop_swarm / mpmat
endif
elseif (np_const/=0.0) then
if (lparticles_number) then
if (lroot) print*, 'particles_initialize_modules: '// &
'can not use np_const together with Particles_number module'
call fatal_error('particles_initialize_modules','')
endif
if (.not.lparticles_radius) then
if (mpmat==0.0) then
if (lroot) print*, 'particles_initialize_modules: '// &
'must have mpmat non zero when setting np_const'
call fatal_error('particles_initialize_modules','')
endif
endif
np_swarm=np_const/(real(npar)/nwgrid)
rhop_swarm=np_swarm*mpmat
else
if (rhop_swarm==0.0) rhop_swarm=mpmat*np_swarm
endif
!
if (lparticles_radius .and. rhopmat>0.0 .and. &
(np_swarm>0.0 .or. lparticles_number .or. particles_module .eq. "lagrangian")) lignore_rhop_swarm=.true.
!
! Initialize individual modules.
!
call initialize_particles_mpicomm (f)
call initialize_particles (f,fp)
call initialize_particles_map
call initialize_particles_adaptation (f)
call initialize_particles_density (f)
call initialize_particles_number (f)
call initialize_particles_radius (f,fp)
call initialize_particles_grad (f)
call initialize_particles_selfgrav (f)
call initialize_particles_sink (f)
call initialize_particles_spin (f)
call initialize_particles_stalker (f)
call initialize_particles_TT (f)
call initialize_particles_mass (f)
call initialize_particles_drag
call initialize_particles_ads (f)
call initialize_particles_surf (f)
call initialize_particles_coag (f)
call initialize_particles_cond (f)
call initialize_particles_collisions (f)
call initialize_pars_diagnos_state (f)
call initialize_particles_diagnos_dv (f)
!
if (lparticles_blocks) then
if (lrun) then
if (lroot) print*, 'particles_initialize_modules: reblocking particles'
call boundconds_particles(fp,ipar)
call map_nearest_grid(fp,ineargrid)
call sort_particles_iblock(fp,ineargrid,ipar)
call map_xxp_grid(f,fp,ineargrid)
call load_balance_particles(f,fp,ipar)
else
inearblock=0
endif
endif
!
! Stop if rhop_swarm is zero.
!
if (irhop/=0 .and. rhop_swarm==0.0 .and. (.not.lparticles_density) &
.and. (.not.lignore_rhop_swarm)) then
if (lroot) then
print*, 'particles_initialize_modules: rhop_swarm is zero'
print*, 'particles_initialize_modules: '// &
'np_swarm, mpmat, rhop_swarm=', np_swarm, mpmat, rhop_swarm
endif
call fatal_error('particles_initialize_modules','rhop_swarm is zero')
endif
!
! Make sure all requested interpolation variables are available.
!
call interpolation_consistency_check
!
! Set internal and external radii of particles.
!
if (rp_int==-impossible .and. r_int>epsi) rp_int=r_int
if (rp_ext==-impossible) rp_ext=r_ext
!
! Disable impossible spectra.
!
if (inp==0) np_spec=.false.
if (any(iapn==0)) np_ap_spec=.false.
if (irhop==0) rhop_spec=.false.
endsubroutine particles_initialize_modules
!***********************************************************************
subroutine particles_init(f)
!
! Set up initial condition for particle modules.
!
! 07-jan-05/anders: coded
!
real, dimension(mx,my,mz,mfarray), intent(inout) :: f
!
! Added somewhat awkward possibility to swap radius and number to the end by
! setting lswap_radius_and_number=.true. (used for lucky droplet experiment)
!
if (.not.lswap_radius_and_number) then
if (lparticles_radius) call set_particle_radius(f,fp,1,npar_loc,init=.true.)
endif
if (lparticles_grad) call set_particle_grad(f,fp,1,npar_loc,ineargrid,init=.true.)
if (.not.lswap_radius_and_number) then
if (lparticles_number) call init_particles_number(f,fp)
endif
if (lparticles_density) call init_particles_density(f,fp)
call init_particles(f,fp,ineargrid)
if (lparticles_sink) call init_particles_sink(f,fp)
if (lparticles_spin) call init_particles_spin(f,fp)
if (lparticles_temperature) call init_particles_TT(f,fp)
if (lparticles_mass) call init_particles_mass(f,fp)
if (lparticles_drag) call init_particles_drag(f,fp)
if (lparticles_adsorbed) call init_particles_ads(f,fp)
if (lparticles_surfspec) call init_particles_surf(f,fp,ineargrid)
if (lparticles_diagnos_state) call init_particles_diagnos_state(fp)
if (lparticles_lyapunov) call init_particles_lyapunov(fp)
if (lswap_radius_and_number) then
if (lparticles_radius) call set_particle_radius(f,fp,1,npar_loc,init=.true.)
if (lparticles_number) call init_particles_number(f,fp)
endif
!
endsubroutine particles_init
!***********************************************************************
subroutine particles_finalize
!
! Finalize particle modules.
!
! 01-May-2019/PABourdin: coded
!
if (lparticles_stalker) call finalize_particles_stalker
!
endsubroutine particles_finalize
!***********************************************************************
subroutine read_snapshot_particles
!
call particles_read_snapshot("pvar.dat")
if (lparticles_lyapunov) call init_particles_lyapunov(fp)
if (lparticles_caustics) call reinitialize_caustics(fp)
if (lparticles_tetrad) call reinitialize_tetrad(fp)
!
endsubroutine read_snapshot_particles
!***********************************************************************
subroutine write_dim_particles(datadir)
!
character (len=*) :: datadir
!
call particles_write_pdim(trim(datadir)//'/pdim.dat')
call particles_write_block(trim(datadir)//'/bdim.dat')
!
endsubroutine write_dim_particles
!***********************************************************************
subroutine read_all_particles_init_pars
!
use File_io, only: read_namelist
!
call read_namelist(read_particles_init_pars ,'particles' ,lparticles)
call read_namelist(read_particles_rad_init_pars ,'particles_radius' ,lparticles_radius)
call read_namelist(read_particles_cond_init_pars ,'particles_cond' ,lparticles_condensation)
call read_namelist(read_particles_spin_init_pars ,'particles_spin' ,lparticles_spin)
call read_namelist(read_particles_sink_init_pars ,'particles_sink' ,lparticles_sink)
call read_namelist(read_particles_num_init_pars ,'particles_number' ,lparticles_number)
call read_namelist(read_particles_dens_init_pars ,'particles_dens' ,lparticles_density)
call read_namelist(read_particles_selfg_init_pars,'particles_selfgrav',lparticles_selfgravity)
call read_namelist(read_particles_mass_init_pars ,'particles_mass' ,lparticles_mass)
call read_namelist(read_particles_drag_init_pars ,'particles_drag' ,lparticles_drag)
call read_namelist(read_particles_TT_init_pars ,'particles_TT' ,lparticles_temperature)
call read_namelist(read_particles_ads_init_pars ,'particles_ads' ,lparticles_adsorbed)
call read_namelist(read_particles_surf_init_pars ,'particles_surf' ,lparticles_surfspec)
call read_namelist(read_particles_chem_init_pars ,'particles_chem' ,lparticles_chemistry)
call read_namelist(read_pstalker_init_pars ,'particles_stalker' ,lparticles_stalker)
call read_namelist(read_plyapunov_init_pars ,'particles_lyapunov',lparticles_lyapunov)
!
endsubroutine read_all_particles_init_pars
!***********************************************************************
subroutine read_all_particles_run_pars
!
use File_io, only: read_namelist
!
call read_namelist(read_particles_run_pars ,'particles' ,lparticles)
call read_namelist(read_particles_adapt_run_pars ,'particles_adapt' ,lparticles_adaptation)
call read_namelist(read_particles_rad_run_pars ,'particles_radius' ,lparticles_radius)
call read_namelist(read_particles_spin_run_pars ,'particles_spin' ,lparticles_spin)
call read_namelist(read_particles_sink_run_pars ,'particles_sink' ,lparticles_sink)
call read_namelist(read_particles_num_run_pars ,'particles_number' ,lparticles_number)
call read_namelist(read_particles_selfg_run_pars ,'particles_selfgrav' ,lparticles_selfgravity)
call read_namelist(read_particles_coag_run_pars ,'particles_coag' ,lparticles_coagulation)
call read_namelist(read_particles_cond_run_pars ,'particles_cond' ,lparticles_condensation)
call read_namelist(read_particles_coll_run_pars ,'particles_coll' ,lparticles_collisions)
call read_namelist(read_particles_stir_run_pars ,'particles_stirring' ,lparticles_stirring)
call read_namelist(read_pstalker_run_pars ,'particles_stalker' ,lparticles_stalker)
call read_namelist(read_pars_diagnos_dv_run_pars ,'particles_diagnos_dv' ,lparticles_diagnos_dv)
call read_namelist(read_pars_diag_state_run_pars ,'particles_diagnos_state',lparticles_diagnos_state)
call read_namelist(read_particles_mass_run_pars ,'particles_mass' ,lparticles_mass)
call read_namelist(read_particles_drag_run_pars ,'particles_drag' ,lparticles_drag)
call read_namelist(read_particles_TT_run_pars ,'particles_TT' ,lparticles_temperature)
call read_namelist(read_particles_ads_run_pars ,'particles_ads' ,lparticles_adsorbed)
call read_namelist(read_particles_surf_run_pars ,'particles_surf' ,lparticles_surfspec)
call read_namelist(read_particles_chem_run_pars ,'particles_chem' ,lparticles_chemistry)
call read_namelist(read_plyapunov_run_pars ,'particles_lyapunov' ,lparticles_lyapunov)
!
endsubroutine read_all_particles_run_pars
!***********************************************************************
subroutine particles_read_snapshot(filename)
!
! Read particle snapshot from file.
!
! 07-jan-05/anders: coded
!
character (len=*) :: filename
!
call input_particles(filename,fp,ipar)
!
endsubroutine particles_read_snapshot
!***********************************************************************
subroutine write_snapshot_particles(f,enum,snapnum)
!
use General, only: itoa
!
real, dimension (mx,my,mz,mfarray) :: f
logical :: enum
integer, optional :: snapnum
!
if (present(snapnum)) then
call particles_write_snapshot('PVAR'//itoa(snapnum),f,enum=.false.,FLIST='pvarN.list')
!
elseif (enum) then
call particles_write_snapshot('PVAR',f,ENUM=.true.,FLIST='pvarN.list')
!
else
call particles_write_snapshot('pvar.dat',f,enum=.false.)
!
endif
!
endsubroutine write_snapshot_particles
!***********************************************************************
subroutine particles_write_snapshot(chsnap,f,enum,flist)
!
! Write particle snapshot to file.
!
! 07-jan-05/anders: coded
!
character (len=*) :: chsnap
real, dimension (mx,my,mz,mfarray) :: f
logical :: enum
character (len=*), optional :: flist
!
logical :: lsnap
!
if (present(flist)) then
call wsnap_particles(chsnap,f,fp,enum,lsnap,dsnap_par_minor,dsnap_par,ipar,mparray,flist)
else
call wsnap_particles(chsnap,f,fp,enum,lsnap,dsnap_par_minor,dsnap_par,ipar,mparray)
endif
!
endsubroutine particles_write_snapshot
!***********************************************************************
subroutine particles_write_dsnapshot(chsnap,f)
!
! Write particle derivative snapshot to file.
!
! 07-jan-05/anders: coded
!
character (len=*) :: chsnap
real, dimension (mx,my,mz,mfarray) :: f
!
call wsnap_particles(chsnap,f,dfp,.false.,.false.,0.0,0.0,ipar,mpvar,nobound=.true.)
!
endsubroutine particles_write_dsnapshot
!***********************************************************************
subroutine particles_write_pdim(filename)
!
! Write npar and mpvar to file.
!
! 09-jan-05/anders: coded
! 10-oct-14/jonas,nils: added the number of particle auxiliaries
!
character (len=*) :: filename
!
open(1,file=filename)
write(1,'(4i9)') npar, mpvar, npar_stalk, mpaux
close(1)
!
endsubroutine particles_write_pdim
!***********************************************************************
subroutine particles_write_block(filename)
!
! Write block domain decomposition parameters to file.
!
! 05-nov-09/anders: coded
!
character (len=*) :: filename
!
if (lparticles_blocks) then
open(1,file=filename)
write(1,'(4i9)') nbrickx, nbricky, nbrickz, nblockmax
write(1,'(4i9)') mxb, myb, mzb, nghostb
write(1,'(3i9)') nxb, nyb, nzb
write(1,'(6i9)') l1b, l2b, m1b, m2b, n1b, n2b
close(1)
endif
!
endsubroutine particles_write_block
!***********************************************************************
subroutine particles_timestep_first(f)
!
! Setup dfp in the beginning of each itsub.
!
! 07-jan-05/anders: coded
!
real, dimension (mx,my,mz,mfarray) :: f
!
call keep_compiler_quiet(f)
!
if (lfirst) then
dfp(1:npar_loc,:)=0.0
else
dfp(1:npar_loc,:)=alpha_ts(itsub)*dfp(1:npar_loc,:)
endif
!
! Insert new and remove old particles continuously during the run
!
if (lfirst) then
call particles_insert_continuously(f)
call particles_remove_continuously(f)
endif
!
endsubroutine particles_timestep_first
!***********************************************************************
subroutine particles_timestep_second(f)
!
! Time evolution of particle variables.
!
! 07-jan-05/anders: coded
!
real, dimension (mx,my,mz,mfarray) :: f
integer :: k
!
! Use zero particle radius to identify tracer particles.
!
if (lparticles_radius) then
do k=1,npar_loc
if (fp(k,iap)==0.0) dfp(k,ivpx:ivpz)=0.0
enddo
endif
!
! Evolve particle state.
!
if (.not.lpointmasses) &
fp(1:npar_loc,1:mpvar) = fp(1:npar_loc,1:mpvar) + dt_beta_ts(itsub)*dfp(1:npar_loc,1:mpvar)
!
! Discrete particle collisions. Must be done at the end of the time-step.
! This call also sorts the particles into mn
!
call particles_discrete_collisions
!
! Adapt the number of particles in each grid cell to a desired number
!
if (lparticles_adaptation .and. llast) then
call sort_particles_imn(fp,ineargrid,ipar,f=f)
call map_xxp_grid(f,fp,ineargrid)
call map_vvp_grid(f,fp,ineargrid)
call particles_adaptation_pencils(f,fp,dfp,ipar,ineargrid)
endif
!
! Insert "old" (removed) particles.
!
call insert_particles_now(f)
!
endsubroutine particles_timestep_second
!***********************************************************************
subroutine split_update_particles(f, dt)
!
! Wrapper for operator split terms for particle dynamics.
!
! 08-may-16/ccyang: coded.
!
real, dimension(mx,my,mz,mfarray), intent(inout) :: f
real, intent(in) :: dt
!
drag: if (lparticles_drag) then
call particles_boundconds(f)
call integrate_drag(f, fp, dt)
endif drag
!
endsubroutine split_update_particles
!***********************************************************************
subroutine particles_discrete_collisions
!
! Discrete particle collisions.
!
! 13-nov-09/anders: coded
!
if ( lparticles_stirring .and. llast ) then
call particle_stirring(fp,ineargrid)
endif
!
if ( (lparticles_collisions .or. lparticles_coagulation .or. &
lparticles_condensation) .and. llast ) then
!
call boundconds_particles(fp,ipar)
call map_nearest_grid(fp,ineargrid)
!
if (lparticles_blocks) then
call sort_particles_iblock(fp,ineargrid,ipar)
if (lparticles_collisions) then
call particles_collisions_blocks(fp,ineargrid)
endif
if (lparticles_coagulation) then
call particles_coagulation_blocks(fp,ineargrid)
endif
else
call sort_particles_imn(fp,ineargrid,ipar)
if (lparticles_collisions) then
call particles_collisions_pencils(fp,ineargrid)
endif
if (lparticles_coagulation) then
call particles_coagulation_pencils(fp,ineargrid)
endif
if (lparticles_condensation) then
call particles_condensation_pencils(fp,ineargrid)
endif
endif
!
endif
!
endsubroutine particles_discrete_collisions
!***********************************************************************
subroutine particles_load_balance(f)
!
! Redistribute particles among the processors for better load balancing.
!
! 04-nov-09/anders: coded
!
real, dimension (mx,my,mz,mfarray) :: f
!
if (lparticles_blocks .and. mod(it,it1_loadbalance)==0) then
call particles_boundconds(f)
call load_balance_particles(f,fp,ipar)
endif
!
endsubroutine particles_load_balance
!***********************************************************************
subroutine particles_boundconds(f)
!
! Particle boundary conditions and parallel communication.
!
! 16-feb-06/anders: coded
!
real, dimension (mx,my,mz,mfarray) :: f
!
! First apply boundary conditions to the newly updated particle positions.
!
call boundconds_particles(fp,ipar,dfp=dfp)
!
! Remove particles that are too close to sink particles or sink points.
! WARNING: ineargrid and the mapped particle density have not been updated
! yet, and the sink particle subroutine must not rely on those arrays.
!
if (.not.lpencil_check_at_work) then
if (lparticles) &
call remove_particles_sink_simple(f,fp,dfp,ineargrid)
if (lparticles_sink) call remove_particles_sink(f,fp,dfp,ineargrid)
endif
!
! Find nearest grid point for each particle.
!
call map_nearest_grid(fp,ineargrid)
!
! Sort particles so that they can be accessed contiguously in the memory.
!
if (lparticles_blocks) then
call sort_particles_iblock(fp,ineargrid,ipar,dfp=dfp)
else
call sort_particles_imn(fp,ineargrid,ipar,dfp=dfp,f=f)
endif
!
! Map the particle positions and velocities on the grid.
!
call map_xxp_grid(f,fp,ineargrid)
call map_vvp_grid(f,fp,ineargrid)
!
endsubroutine particles_boundconds
!***********************************************************************
subroutine particles_calc_selfpotential(f,rhs_poisson,lcontinued)
!
! Calculate the potential of the dust particles (wrapper).
!
! 13-jun-06/anders: coded
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (nx,ny,nz) :: rhs_poisson
logical :: lcontinued
!
call calc_selfpotential_particles(f,rhs_poisson,lcontinued)
call calc_selfpot_sinkparticles(f,rhs_poisson,fp,ineargrid)
!
endsubroutine particles_calc_selfpotential
!***********************************************************************
subroutine particles_before_boundary(f)
!
! Calculate particle-related properties before boundary conditions are
! set.
!
! 07-feb-09/anders: coded
!
real, dimension (mx,my,mz,mfarray) :: f
!
if (lparticles) then
call particles_dragforce_stiff(f,fp,ineargrid)
call periodic_boundcond_on_aux(f)
endif
if (lparticles_caustics) call reset_caustics(fp)
!
endsubroutine particles_before_boundary
!***********************************************************************
subroutine particles_special_bfre_bdary(f)
!
! Fetch fp array to special module.
!
! 01-mar-08/wlad: coded
!
use Special, only: special_particles_bfre_bdary
!
real, dimension (mx,my,mz,mfarray) :: f
!
call special_particles_bfre_bdary(f,fp,ineargrid)
!
endsubroutine particles_special_bfre_bdary
!***********************************************************************
subroutine particles_special_after_dtsub(f, dtsub)
!
! Send fp to Special for processing in the end of a sub-time-step.
!
! 28-aug-18/ccyang: coded
!
use Special, only: special_particles_after_dtsub
!
real, dimension(mx,my,mz,mfarray), intent(in) :: f
real, intent(in) :: dtsub
!
call special_particles_after_dtsub(f, dtsub, fp, dfp, ineargrid)
!
endsubroutine particles_special_after_dtsub
!***********************************************************************
subroutine particles_pencil_criteria
!
! Request pencils for particles.
!
! 20-apr-06/anders: coded
!
if (lparticles) call pencil_criteria_particles
if (lparticles_radius) call pencil_criteria_par_radius
if (lparticles_spin) call pencil_criteria_par_spin
if (lparticles_number) call pencil_criteria_par_number
if (lparticles_density) call pencil_criteria_par_density
if (lparticles_selfgravity) call pencil_criteria_par_selfgrav
if (lparticles_temperature) call pencil_criteria_par_TT
if (lparticles_mass) call pencil_criteria_par_mass
if (lparticles_adsorbed) call pencil_criteria_par_ads
if (lparticles_chemistry) call pencil_criteria_par_chem
!
endsubroutine particles_pencil_criteria
!***********************************************************************
subroutine particles_pencil_interdep(lpencil_in)
!
! Calculate particle pencils.
!
! 15-feb-06/anders: coded
!
logical, dimension(npencils) :: lpencil_in
!
if (lparticles) call pencil_interdep_particles(lpencil_in)
if (lparticles_selfgravity) call pencil_interdep_par_selfgrav(lpencil_in)
!
endsubroutine particles_pencil_interdep
!***********************************************************************
subroutine particles_calc_pencils(f,p)
!
! Calculate particle pencils.
!
! 14-feb-06/anders: coded
!
real, dimension (mx,my,mz,mfarray) :: f
type (pencil_case) :: p
!
if (lparticles) call calc_pencils_particles(f,p)
if (lparticles_lyapunov) call calc_pencils_par_lyapunov(f,p)
if (lparticles_selfgravity) call calc_pencils_par_selfgrav(f,p)
if (lparticles_chemistry) call calc_pencils_par_chem(f,p)
!
endsubroutine particles_calc_pencils
!***********************************************************************
subroutine particles_pde(f,df,p)
!
! Dynamical evolution of particle variables.
!
! 07-jan-05/anders: coded
!
use Mpicomm
use Special, only: special_calc_particles
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mx,my,mz,mvar) :: df
type (pencil_case) :: p
!
intent (inout) :: f
intent (out) :: df
intent (in) :: p
!
! Write information about local particle environment to file.
!
if (lfirst .and. (.not.lpencil_check_at_work)) &
call particles_stalker_sub(f,fp,ineargrid)
!
! Dynamical equations.
!
if (lparticles) call dxxp_dt(f,df,fp,dfp,ineargrid)
if (lparticles) call dvvp_dt(f,df,p,fp,dfp,ineargrid)
if (lparticles_lyapunov) call dlyapunov_dt(f,df,fp,dfp,ineargrid)
if (lparticles_radius) call dap_dt(f,df,fp,dfp,ineargrid)
if (lparticles_spin) call dps_dt(f,df,fp,dfp,ineargrid)
if (lparticles_mass) call dpmass_dt(f,df,fp,dfp,ineargrid)
if (lparticles_temperature) call dpTT_dt(f,df,fp,dfp,ineargrid)
if (lparticles_adsorbed) call dpads_dt(f,df,fp,dfp,ineargrid)
if (lparticles_surfspec) call dpsurf_dt(f,df,fp,dfp,ineargrid)
if (lparticles_number) call dnpswarm_dt(f,df,fp,dfp,ineargrid)
if (lparticles_density) call drhopswarm_dt(f,df,fp,dfp,ineargrid)
if (lparticles_selfgravity) call dvvp_dt_selfgrav(f,df,fp,dfp,ineargrid)
!
! Interface for your personal subroutines calls
!
if (lspecial) call special_calc_particles(f,df,fp,dfp,ineargrid)
!
! Create new sink particles or sink points.
!
if (.not.lpencil_check_at_work) then
if (lparticles) &
call create_particles_sink_simple(f,fp,dfp,ineargrid)
if (lparticles_sink) call create_particles_sink(f,fp,dfp,ineargrid)
endif
!
! Correct for curvilinear geometry.
!
if (lparticles) call correct_curvilinear
!
! Output particle size distribution to file.
!
if (lparticles_radius .and. loutput_psize_dist .and. ldiagnos .and. &
.not. lpencil_check_at_work) call output_particle_size_dist(fp)
!
endsubroutine particles_pde
!***********************************************************************
subroutine particles_pde_pencil(f,df,p)
!
! Dynamical evolution of particle variables in pencils.
!
! 20-apr-06/anders: coded
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mx,my,mz,mvar) :: df
type (pencil_case) :: p
!
intent (in) :: p
intent (inout) :: f, df
!
! Create shepherd/neighbour list of required.
!
call timing('particles_pde_pencil','entered',mnloop=.true.)
if (lshepherd_neighbour) &
call shepherd_neighbour_pencil(fp,ineargrid,kshepherd,kneighbour)
!
! Interpolate required quantities using the predefined policies. Variables
! are found in interp.
! (Clean-up should be performed at end of this subroutine!)
!
call interpolate_quantities(f,fp,p,ineargrid)
!
! If reactive particles are enabled, needed quantities are calculated
!
if (lparticles_chemistry) then
call calc_pchemistry_pencils(f,fp,p,ineargrid)
endif
if (lparticles_surfspec) then
call calc_psurf_pencils(f,fp,p,ineargrid)
endif
!
! Dynamical equations.
!
if (lparticles) call dxxp_dt_pencil(f,df,fp,dfp,p,ineargrid)
if (lparticles) call dvvp_dt_pencil(f,df,fp,dfp,p,ineargrid)
if (lparticles_lyapunov) call dlyapunov_dt_pencil(f,df,fp,dfp,p,ineargrid)
if (lparticles_radius) call dap_dt_pencil(f,df,fp,dfp,p,ineargrid)
if (lparticles_spin) call dps_dt_pencil(f,df,fp,dfp,p,ineargrid)
if (lparticles_mass) call dpmass_dt_pencil(f,df,fp,dfp,p,ineargrid)
if (lparticles_temperature) call dpTT_dt_pencil(f,df,fp,dfp,p,ineargrid)
if (lparticles_adsorbed) call dpads_dt_pencil(f,df,fp,dfp,p,ineargrid)
if (lparticles_surfspec) call dpsurf_dt_pencil(f,df,fp,dfp,p,ineargrid)
if (lparticles_number) call dnpswarm_dt_pencil(f,df,fp,dfp,p,ineargrid)
if (lparticles_density) call drhopswarm_dt_pencil(f,df,fp,dfp,p,ineargrid)
if (lparticles_selfgravity) &
call dvvp_dt_selfgrav_pencil(f,df,fp,dfp,p,ineargrid)
! if (lparticles_polymer) &
! call dRR_dt_pencil(f,df,fp,dfp,ineargrid)
!
! Time-step contribution from discrete particle collisions.
!
if (lparticles_collisions) &
call particles_collisions_timestep(fp,ineargrid)
if (lparticles_coagulation) &
call particles_coagulation_timestep(fp,ineargrid)
!
call cleanup_chemistry_pencils
call cleanup_surf_pencils
call cleanup_interpolated_quantities
call timing('particles_pde_pencil','finished',mnloop=.true.)
!
endsubroutine particles_pde_pencil
!***********************************************************************
subroutine particles_pde_blocks(f,df)
!
! Dynamical evolution of particle variables in blocks.
!
! 29-nov-09/anders: coded
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mx,my,mz,mvar) :: df
!
intent (inout) :: f, df
!
if (lparticles_blocks) then
!
! Fill adopted blocks with gas density and gas velocity field, for calculating
! drag forces.
!
if (lfill_blocks_density) then
if (ldensity_nolog) then
call fill_blocks_with_bricks(f,fb,mfarray,irho)
else
call fill_blocks_with_bricks(f,fb,mfarray,ilnrho)
endif
endif
!
if (lfill_blocks_velocity) then
call fill_blocks_with_bricks(f,fb,mfarray,iux)
call fill_blocks_with_bricks(f,fb,mfarray,iuy)
call fill_blocks_with_bricks(f,fb,mfarray,iuz)
endif
!
! Fill adopted blocks with gravitational acceleration.
!
if (lfill_blocks_gpotself) then
call fill_blocks_with_bricks(f,fb,mfarray,igpotselfx)
call fill_blocks_with_bricks(f,fb,mfarray,igpotselfy)
call fill_blocks_with_bricks(f,fb,mfarray,igpotselfz)
endif
!
! Zero the block contribution to time evolution of gas velocity.
!
if (lhydro) dfb(:,:,:,iux:iuz,0:nblock_loc-1)=0.0
!
! Dynamical equations.
!
if (lparticles) call dxxp_dt_blocks(f,df,fp,dfp,ineargrid)
if (lparticles) call dvvp_dt_blocks(f,df,fp,dfp,ineargrid)
!
! Add block contribution to df to the main grid.
!
if (lfill_bricks_velocity) then
call fill_bricks_with_blocks(df,dfb,mvar,iux)
call fill_bricks_with_blocks(df,dfb,mvar,iuy)
call fill_bricks_with_blocks(df,dfb,mvar,iuz)
endif
!
endif
!
endsubroutine particles_pde_blocks
!***********************************************************************
subroutine correct_curvilinear
!
! Curvilinear corrections to acceleration only.
! Corrections to velocity were already taken into account
! in the dxx_dt of particles_dust.f90
!
! In the case that the N-body code is used, the update in polar grids
! in done by transforming the variables first to Cartesian, to achieve a
! better conservation of the Jacobi constant, and this code is not called.
!
! 15-sep-07/wlad: coded
!
real :: rpcyl1,rp1,lat,costhp,sinthp,sin1thp,cotthp
integer :: k
!
! The subroutine does not need to be called if either pointmasses
! or particle tracers is used. Pointmasses correct for curvilinear
! coordinates already, and particle tracers do not have particle
! velocities as they follow the gas. The npvar>3 means that only
! ip[xyz] exist. A logical lparticles_tracers is missing.
!
if ((.not.lpointmasses).and.(npvar>3)) then
do k=1,npar_loc
!
! Correct acceleration.
!
if (lcylindrical_coords) then
!
rpcyl1 = 1./max(fp(k,ixp),tini)
!
dfp(k,ivpx) = dfp(k,ivpx) + rpcyl1*fp(k,ivpy)**2
dfp(k,ivpy) = dfp(k,ivpy) - rpcyl1*fp(k,ivpx)*fp(k,ivpy)
!
elseif (lspherical_coords) then
!
rp1 = 1./max(fp(k,ixp),tini)
if (luse_latitude) then
lat=pi/2-fp(k,iyp)
costhp=sin(lat)
else
costhp=cos(fp(k,iyp))
endif
sinthp=sin(fp(k,iyp))
if (abs(sinthp)>tini) then
sin1thp=1./sinthp
else
sin1thp=0.
endif
cotthp=costhp*sin1thp
!
dfp(k,ivpx) = dfp(k,ivpx) + rp1*(fp(k,ivpy)**2 + fp(k,ivpz)**2)
dfp(k,ivpy) = dfp(k,ivpy) - rp1*(fp(k,ivpy)*fp(k,ivpx) - fp(k,ivpz)**2*cotthp)
dfp(k,ivpz) = dfp(k,ivpz) - rp1*(fp(k,ivpz)*fp(k,ivpx) + fp(k,ivpz)*fp(k,ivpy)*cotthp)
!
endif
enddo
endif
!
endsubroutine correct_curvilinear
!***********************************************************************
subroutine write_all_particles_init_pars(unit)
!
! Write particle start parameters to file.
!
integer, intent (in) :: unit
!
call write_particles_init_pars(unit)
if (lparticles_radius) call write_particles_rad_init_pars(unit)
if (lparticles_condensation)call write_particles_cond_init_pars(unit)
if (lparticles_spin) call write_particles_spin_init_pars(unit)
if (lparticles_sink) call write_particles_sink_init_pars(unit)
if (lparticles_number) call write_particles_num_init_pars(unit)
if (lparticles_density) call write_particles_dens_init_pars(unit)
if (lparticles_selfgravity) call write_particles_selfg_init_pars(unit)
if (lparticles_stalker) call write_pstalker_init_pars(unit)
if (lparticles_mass) call write_particles_mass_init_pars(unit)
if (lparticles_drag) call write_particles_drag_init_pars(unit)
if (lparticles_temperature) call write_particles_TT_init_pars(unit)
if (lparticles_adsorbed) call write_particles_ads_init_pars(unit)
if (lparticles_surfspec) call write_particles_surf_init_pars(unit)
if (lparticles_chemistry) call write_particles_chem_init_pars(unit)
if (lparticles_lyapunov) call write_plyapunov_init_pars(unit)
!
endsubroutine write_all_particles_init_pars
!***********************************************************************
subroutine write_all_particles_run_pars(unit)
!
! Write particle run parameters to file.
!
integer, intent (in) :: unit
!
if (lparticles) call write_particles_run_pars(unit)
if (lparticles_radius) call write_particles_rad_run_pars(unit)
if (lparticles_spin) call write_particles_spin_run_pars(unit)
if (lparticles_sink) call write_particles_sink_run_pars(unit)
if (lparticles_number) call write_particles_num_run_pars(unit)
if (lparticles_selfgravity) call write_particles_selfg_run_pars(unit)
if (lparticles_coagulation) call write_particles_coag_run_pars(unit)
if (lparticles_condensation) call write_particles_coag_run_pars(unit)
if (lparticles_collisions) call write_particles_coll_run_pars(unit)
if (lparticles_stirring) call write_particles_stir_run_pars(unit)
if (lparticles_stalker) call write_pstalker_run_pars(unit)
if (lparticles_diagnos_dv) call write_pars_diagnos_dv_run_pars(unit)
if (lparticles_diagnos_state) call write_pars_diag_state_run_pars(unit)
if (lparticles_mass) call write_particles_mass_run_pars(unit)
if (lparticles_drag) call write_particles_drag_run_pars(unit)
if (lparticles_temperature) call write_particles_TT_run_pars(unit)
if (lparticles_adsorbed) call write_particles_ads_run_pars(unit)
if (lparticles_surfspec) call write_particles_surf_run_pars(unit)
if (lparticles_chemistry) call write_particles_chem_run_pars(unit)
if (lparticles_lyapunov) call write_plyapunov_run_pars(unit)
!
endsubroutine write_all_particles_run_pars
!***********************************************************************
subroutine wsnap_particles(snapbase,f,fp,enum,lsnap,dsnap_par_minor,dsnap_par,ipar,varsize,flist,nobound)
!
! Write particle snapshot file, labelled consecutively if enum==.true.
! Otherwise just write a snapshot without label (used e.g. for pvar.dat)
!
! 29-dec-04/anders: adapted from wsnap
! 04-oct-08/ccyang: use a separate log file for minor snapshots
! 26-nov-08/ccyang: add independent sequence for particle snapshots
!
use General, only: itoa, safe_character_assign
use IO, only: log_filename_to_file
use Particles_mpicomm, only: output_blocks
use Sub
!
integer :: varsize
character (len=*) :: snapbase, flist
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mpar_loc,varsize) :: fp
logical :: enum, lsnap, nobound
real :: dsnap_par_minor, dsnap_par
integer, dimension (mpar_loc) :: ipar
!
logical, save :: lfirst_call=.true.
integer, save :: nsnap, nsnap_minor, nsnap_par
real, save :: tsnap, tsnap_minor, tsnap_par
character (len=fnlen), save :: fmajor, fminor, fpar
logical :: lsnap_minor=.false., lsnap_par=.false.
character (len=fnlen) :: snapname
character (len=intlen) :: nsnap_ch,nsnap_minor_ch,nsnap_par_ch
!
optional :: flist, nobound
!
! Output snapshot with label in 'tsnap' time intervals
! file keeps the information about number and time of last snapshot
!
if (enum) then
!
! At first call, need to initialize tsnap.
! tsnap calculated in read_snaptime, but only available to root processor.
!
if (lfirst_call) then
call safe_character_assign(fmajor,trim(datadir)//'/tsnap.dat')
call read_snaptime(fmajor,tsnap,nsnap,dsnap,t)
if (dsnap_par_minor>0.0) then
call safe_character_assign(fminor,trim(datadir)//'/tsnap_minor.dat')
call read_snaptime(fminor,tsnap_minor,nsnap_minor,dsnap_par_minor,t)
endif
if (dsnap_par>0.0) then
call safe_character_assign(fpar,trim(datadir)//'/tsnap_par.dat')
call read_snaptime(fpar,tsnap_par,nsnap_par,dsnap_par,t)
endif
lfirst_call=.false.
endif
!
! Output independent sequence of particle snapshots.
!
if (dsnap_par>0.0) then
call update_snaptime(fpar,tsnap_par,nsnap_par,dsnap_par,t,lsnap_par,nsnap_par_ch)
if (lsnap_par) then
snapname=trim(snapbase)//'_'//nsnap_par_ch
call particles_boundconds(f)
call output_particles(snapname,fp,ipar)
if (ip<=10 .and. lroot) &
print*,'wsnap_particles: written snapshot ', snapname
if (present(flist)) call log_filename_to_file(snapname,flist)
endif
endif
!
! Possible to output minor particle snapshots (e.g. for a movie).
!
if (dsnap_par_minor>0.0) then
call update_snaptime(fminor,tsnap_minor,nsnap_minor,dsnap_par_minor,t,lsnap_minor,nsnap_minor_ch)
if (lsnap_minor) then
snapname=trim(snapbase)//trim(itoa(nsnap-1))//'.'//nsnap_minor_ch
call particles_boundconds(f)
call output_particles(snapname,fp,ipar)
if (ip<=10 .and. lroot) &
print*,'wsnap_particles: written snapshot ', snapname
if (present(flist)) call log_filename_to_file(snapname,flist)
endif
endif
!
! Regular data snapshots must come synchronized with the fluid snapshots.
!
call update_snaptime(fmajor,tsnap,nsnap,dsnap,t,lsnap,nsnap_ch,nowrite=.true.)
if (lsnap) then
snapname=trim(snapbase)//nsnap_ch
call particles_boundconds(f)
call output_particles(snapname,fp,ipar)
if (lparticles_blocks) &
call output_blocks(trim(directory_dist)//'/BLOCKS'//nsnap_ch)
if (ip<=10 .and. lroot) &
print*,'wsnap_particles: written snapshot ', snapname
if (present(flist)) call log_filename_to_file(snapname,flist)
nsnap_minor=1
endif
!
else
!
! Write snapshot without label
!
snapname=snapbase
if (present(nobound)) then
if (.not. nobound) call particles_boundconds(f)
else
call particles_boundconds(f)
endif
call output_particles(snapname,fp,ipar)
if (lparticles_blocks) &
call output_blocks(trim(directory_dist)//'/blocks.dat')
if (ip<=10 .and. lroot) &
print*,'wsnap_particles: written snapshot ', snapname
if (present(flist)) call log_filename_to_file(snapname,flist)
endif
!
endsubroutine wsnap_particles
!***********************************************************************
subroutine particles_powersnap(f)
!
! Calculate power spectra of particle variables.
!
! 01-jan-06/anders: coded
!
real, dimension (mx,my,mz,mfarray) :: f
!
call powersnap_particles(f)
!
endsubroutine particles_powersnap
!***********************************************************************
subroutine get_slices_particles(f,slices)
!
! Write slices for animation of Particle variables.
!
use Slices_methods, only: assign_slices_vec, assign_slices_scal, process_slices
!
real, dimension (mx,my,mz,mfarray) :: f
type (slice_data) :: slices
integer :: l
!
! Loop over slices
!
select case (trim(slices%name))
!
! Particle number density
!
case ('np'); call assign_slices_scal(slices,f,inp)
!
! Particle mass density
!
case ('rhop')
if (irhop/=0) then
call assign_slices_scal(slices,f,irhop)
else
call assign_slices_scal(slices,f,inp)
if (lcartesian_coords.and.(all(lequidist))) then
call process_slices(slices,rhop_swarm) ! multiply with rhop_swarm
else !MR: both implementations identical!!!
call process_slices(slices,rhop_swarm)
endif
endif
!
! Particle velocity field
!
! One needs to set lcalc_uup = .true. in &particles_init_pars/ and add in src/cparam.local:
!
! MAUX CONTRIBUTION 3
! COMMUNICATED AUXILIARIES 3
!
case ('vvp'); call assign_slices_vec(slices,f,iupx)
!
endselect
!
endsubroutine get_slices_particles
!***********************************************************************
subroutine insert_particles_now(f)
!
! Insert particles which has been removed from the system as the
! system has evolved.
!
! 2012-oct-19/dhruba: coded
!
real, dimension (mx,my,mz,mfarray) :: f
!
if (linsert_particle) call insert_lost_particles(f,fp,ineargrid)
!
endsubroutine insert_particles_now
!***********************************************************************
subroutine particles_stochastic
if (lparticles_lyapunov) call particles_stochastic_lyapunov(fp)
endsubroutine particles_stochastic
!***********************************************************************
subroutine particles_insert_continuously(f)
!
! Insert particles continuously, i.e. add particles in
! the beginning of a time step.
!
! sep-09/kragset: coded
!
real, dimension (mx,my,mz,mfarray) :: f
!
if (linsert_particles_continuously) then
call insert_particles(f,fp,ineargrid)
endif
!
endsubroutine particles_insert_continuously
!***********************************************************************
subroutine particles_remove_continuously(f)
!
! Remove particles continuously, i.e. remove particles in
! the beginning of a time step under a certain criteria.
!
! dec-16/aschreib: coded
!
real, dimension (mx,my,mz,mfarray) :: f
!
real :: rp, epsd2g_tmp
integer :: k, ix0, iy0, iz0
!
if (remove_particle_at_time > 0) then
!
! Should this be moved to particles_dust.f90 and _blocks.f90 respectivly? (signed who?)
!
if (remove_particle_at_time-dt < t) then
if (t < remove_particle_at_time+2*dt) then
select case (remove_particle_criteria)
case ('all')
k=1
do while (k <= npar_loc)
call remove_particle(fp,ipar,k,dfp,ineargrid)
k=k+1
enddo
remove_particle_at_time = -1.
case ('none')
remove_particle_at_time = -1.
case ('sphere')
k=1
do while (k <= npar_loc)
rp = sqrt(fp(k,ixp)**2 + fp(k,iyp)**2 + fp(k,izp)**2)
if ( rp > remove_particle_criteria_size) then
call remove_particle(fp,ipar,k,dfp,ineargrid)
else
k=k+1
endif
enddo
remove_particle_at_time = -1.
case ('density-threshold')
k=1
do while (k <= npar_loc)
! find closest cell
ix0=ineargrid(k,1); iy0=ineargrid(k,2); iz0=ineargrid(k,3)
! dust-to-gas ratio of closest cell
epsd2g_tmp = f(ix0,iy0,iz0,irhop)/f(ix0,iy0,iz0,irho)
! if dust-to-gas-ratio of closest cell is bigger than threshold, remove
if (epsd2g_tmp > remove_particle_criteria_edtog) then
call remove_particle(fp,ipar,k,dfp,ineargrid)
else
k=k+1
endif
enddo
remove_particle_at_time = -1.
case ('xycylinder')
k=1
do while (k<=npar_loc)
rp = sqrt(fp(k,ixp)**2 + fp(k,iyp)**2)
if ( rp > remove_particle_criteria_size ) then
call remove_particle(fp,ipar,k,dfp,ineargrid)
else
k=k+1
endif
enddo
remove_particle_at_time = -1.
endselect
else
remove_particle_at_time = -1.
endif
endif
endif
!
endsubroutine particles_remove_continuously
!***********************************************************************
subroutine particles_cleanup
!
! call particles_final_clean_up
if (lparticles_chemistry) then
call particles_chemistry_clean_up
call particles_surfspec_clean_up
call particles_adsorbed_clean_up
endif
!
endsubroutine particles_cleanup
!***********************************************************************
subroutine fetch_nparloc(nparloc_aux)
integer, intent(out) :: nparloc_aux
nparloc_aux=npar_loc
endsubroutine fetch_nparloc
!***********************************************************************
subroutine append_particle_index(label,ilabel)
!
character (len=*), intent(in) :: label
integer, intent(out) :: ilabel
!
call append_npvar(label,ilabel)
!
endsubroutine append_particle_index
!***********************************************************************
subroutine fetch_fp_array(fp_aux,dfp_aux,ixw,iyw,izw,ivxw,ivyw,ivzw)
!
real, dimension(mpar_loc,mparray), intent(out) :: fp_aux
real, dimension(mpar_loc,mpvar), intent(out) :: dfp_aux
integer, intent(out) :: ixw,iyw,izw,ivxw,ivyw,ivzw
!
fp_aux = fp
dfp_aux = dfp
ixw=ixp ; iyw=iyp; izw=izp
ivxw=ivpx ; ivyw=ivpy; ivzw=ivpz
!
endsubroutine fetch_fp_array
!***********************************************************************
subroutine return_fp_array(fp_aux,dfp_aux,flag)
!
real, dimension(mpar_loc,mparray), intent(in) :: fp_aux
real, dimension(mpar_loc,mpvar), intent(in) :: dfp_aux
logical, dimension(mpar_loc), intent(in), optional :: flag
!
integer :: k
!
fp = fp_aux
dfp = dfp_aux
!
if (present(flag)) then
do k=npar_loc,1,-1
if (flag(k)) call remove_particle(fp,ipar,k,dfp,ineargrid)
enddo
endif
!
endsubroutine return_fp_array
!***********************************************************************
endmodule Particles_main