! $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 not_implemented('particles_initialize_modules','particles 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 on 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: npar/mpar_loc>ncpus')
      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) call fatal_error('particles_initialize_modules', &
                               'cannot use np_const together with Particles_number module')
        if (.not.lparticles_radius) then
          if (mpmat==0.0) &
            call fatal_error('particles_initialize_modules','must have mpmat /=0 when setting np_const')
        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 == "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)
!
      call remove_particle_initialize
!
      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_write_rmv()
!
!  Write the list of removed particles and clear the list.
!
!  28-jan-24/ccyang: coded.
!
      call remove_particle_writelog
!
    endsubroutine particles_write_rmv
!***********************************************************************
    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
!
      if (lparticles_drag) then
        call particles_boundconds(f)
        call integrate_drag(f, fp, dt)
      endif
!
    endsubroutine split_update_particles
!***********************************************************************
    subroutine particles_discrete_collisions
!
!  Discrete particle collisions.
!
!  13-nov-09/anders: coded
!
      if ( lparticles_stirring .and. llast ) call particle_stirring(fp,ineargrid)
!
      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)   call particles_collisions_blocks(fp,ineargrid)
          if (lparticles_coagulation)  call particles_coagulation_blocks(fp,ineargrid)
        else
          call sort_particles_imn(fp,ineargrid,ipar)
          if (lparticles_collisions)   call particles_collisions_pencils(fp,ineargrid)
          if (lparticles_coagulation)  call particles_coagulation_pencils(fp,ineargrid)
          if (lparticles_condensation) call particles_condensation_pencils(fp,ineargrid)
        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,fp,ineargrid)
      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)
!
!  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
!
      intent (inout)  :: f
      intent (out) :: df
!
!  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,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)    call calc_pchemistry_pencils(f,fp,p,ineargrid)
     if (lparticles_surfspec)     call calc_psurf_pencils(f,fp,p,ineargrid)
!
!  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_calc_pencil_diags(p)

      type (pencil_case) :: p

      if (lparticles)             call calc_diagnostics_particles(fp,p,ineargrid)
      if (lparticles_selfgravity) call calc_diagnostics_particles_selg(p)
      if (lparticles_surfspec)    call calc_diagnostics_particles_surf(p)
      if (lparticles_radius)      call calc_diagnostics_particles_rad(p)
      if (lparticles_number)      call calc_diagnostics_particles_num(p)
      if (lparticles_chemistry)   call calc_diagnostics_particles_chem(p)

    endsubroutine particles_calc_pencil_diags
!***********************************************************************
    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("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) then
          if (snapname == "PVAR0") then
            call output_blocks("BLOCKS0")
          else
            call output_blocks("blocks.dat")
          endif
        endif
        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
!
!  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) call insert_particles(f,fp,ineargrid)
      if (lpartnucleation) call insert_nucleii(f,fp,ineargrid)
!
    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.

            case ('size')
              k=1
              do while (k<=npar_loc)
                if ( fp(k,iap) < remove_particle_criteria_size ) then
                  call remove_particle(fp,ipar,k,dfp,ineargrid)
!                  print*,"Removed particle:",k,fp(k,iap)
                else
                  k=k+1
                endif
              enddo
              remove_particle_at_time = t+dt
            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
