! $Id$ ! This module takes care of everything related to tracer particles. ! !** AUTOMATIC CPARAM.INC GENERATION **************************** ! ! Declare (for generation of cparam.inc) the number of f array ! variables and auxiliary variables added by this module ! ! CPARAM logical, parameter :: lparticles=.true. ! CPARAM character (len=20), parameter :: particles_module="tracer" ! ! MAUX CONTRIBUTION 2 ! MPVAR CONTRIBUTION 3 ! MPAUX CONTRIBUTION 3 ! ! PENCILS PROVIDED np; rhop; epsp; rhop_swarm; grhop(3); peh ! !*************************************************************** module Particles ! use Cdata use General, only: keep_compiler_quiet use Messages use Particles_cdata use Particles_map use Particles_sub use Particles_mpicomm ! implicit none ! include 'particles.h' include 'particles_common.h' ! real :: xp0=0.0, yp0=0.0, zp0=0.0, tausp=0.0 real :: nu_epicycle=0.0, nu_epicycle2=0.0 logical :: ldragforce_equi_global_eps=.false. logical :: lquadratic_interpolation=.false. logical :: ltrace_dust=.false. real :: particles_insert_rate=0.,max_particle_insert_time=huge1 character (len=labellen), dimension (ninit) :: initxxp='nothing' character (len=labellen), dimension (ninit) :: insertxxp='nothing' character (len=labellen) :: gravz_profile='zero' logical :: lglobalrandom=.false. ! namelist /particles_init_pars/ & initxxp, xp0, yp0, zp0, bcpx, bcpy, bcpz, eps_dtog, tausp, & ldragforce_equi_global_eps, lquadratic_interpolation, & lparticlemesh_cic, lparticlemesh_tsc, ltrace_dust, & gravz_profile, nu_epicycle, lglobalrandom ! namelist /particles_run_pars/ & bcpx, bcpy, bcpz, lquadratic_interpolation, & lparticlemesh_cic, lparticlemesh_tsc, ltrace_dust, & lcheck_exact_frontier,particles_insert_rate, & linsert_particles_continuously,dsnap_par, & max_particle_insert_time,insertxxp ! integer :: idiag_xpm=0, idiag_ypm=0, idiag_zpm=0 integer :: idiag_xp2m=0, idiag_yp2m=0, idiag_zp2m=0 integer :: idiag_nparmax=0, idiag_npmax=0, idiag_npmin=0 integer :: idiag_npmx=0, idiag_rhopmx=0, idiag_epspmx=0 integer :: idiag_npmz=0, idiag_rhopmz=0, idiag_epspmz=0 integer :: idiag_epsmin=0,idiag_epsmax=0 ! contains !*********************************************************************** subroutine register_particles() ! ! Set up indices for access to the fp and dfp arrays ! ! 29-dec-04/anders: coded ! use FArrayManager ! if (lroot) call svn_id( & "$Id$") ! ! Indices for particle position. ! call append_npvar('ixp',ixp) call append_npvar('iyp',iyp) call append_npvar('izp',izp) ! ! Indices for particle velocity. ! call append_npaux('ivpx',ivpx) call append_npaux('ivpy',ivpy) call append_npaux('ivpz',ivpz) ! ! Set indices for auxiliary variables. ! call farray_register_auxiliary('np',inp) call farray_register_auxiliary('rhop',irhop) ! endsubroutine register_particles !*********************************************************************** subroutine initialize_particles(f,fp) ! ! Perform any post-parameter-read initialization i.e. calculate derived ! parameters. ! ! 29-dec-04/anders: coded ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mpar_loc,mparray), intent (in) :: fp ! real :: rhom ! if (rhop_swarm==0.0 .or. mp_swarm==0.0) then ! For stratification, take into account gas present outside the simulation box. if ((lgravz .and. lgravz_gas) .or. gravz_profile=='linear' ) then rhom=sqrt(2*pi)*1.0*1.0/Lz ! rhom = Sigma/Lz, Sigma=sqrt(2*pi)*H*rho1 else rhom=1.0 endif if (rhop_swarm==0.0) & rhop_swarm = eps_dtog*rhom/(real(npar)/(nxgrid*nygrid*nzgrid)) if (mp_swarm==0.0) & mp_swarm = eps_dtog*rhom*box_volume/(real(npar)) if (lroot) print*, 'initialize_particles: '// & 'dust-to-gas ratio eps_dtog=', eps_dtog endif ! if (lroot) then print*, 'initialize_particles: '// & 'mass per constituent particle mpmat=', mpmat print*, 'initialize_particles: '// & 'mass per superparticle mp_swarm =', mp_swarm print*, 'initialize_particles: '// & 'number density per superparticle np_swarm=', np_swarm print*, 'initialize_particles: '// & 'mass density per superparticle rhop_swarm=', rhop_swarm endif ! ! Calculate nu_epicycle**2 for gravity. ! if (nu_epicycle/=0.0) then gravz_profile='linear' nu_epicycle2=nu_epicycle**2 endif ! call keep_compiler_quiet(f) ! endsubroutine initialize_particles !*********************************************************************** subroutine init_particles(f,fp,ineargrid) ! ! Initial positions and velocities of tracer particles. ! ! 29-dec-04/anders: coded ! use Density, only: beta_glnrho_global use EquationOfState, only: gamma, cs20 use General, only: random_number_wrapper use Mpicomm, only: stop_it use Sub ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mpar_loc,mparray) :: fp integer, dimension (mpar_loc,3) :: ineargrid ! real, dimension (nx) :: eps real, dimension (3) :: Lxyz_par, xyz0_par, xyz1_par real :: r, p, cs integer :: j, k logical :: lnothing ! real :: rad,rad_scl,phi,tmp ! intent (inout) :: f intent (out) :: fp ! ! Use either a local random position or a global random position for certain ! initial conditions. The default is a local random position, but the equal ! number of particles per processors means that this is not completely random. ! if (lglobalrandom) then Lxyz_par=Lxyz xyz0_par=xyz0 xyz1_par=xyz1 else Lxyz_par=Lxyz_loc xyz0_par=xyz0_loc xyz1_par=xyz1_loc endif ! ! Initial particle position. ! lnothing=.false. do j=1,ninit ! select case (initxxp(j)) ! case ('nothing') if (lroot .and. .not. lnothing) print*, 'init_particles: nothing' lnothing=.true. ! case ('origin') if (lroot) print*, 'init_particles: All particles at origin' fp(1:npar_loc,ixp:izp)=0. ! case ('constant') if (lroot) & print*, 'init_particles: All particles at x,y,z=', xp0, yp0, zp0 fp(1:npar_loc,ixp)=xp0 fp(1:npar_loc,iyp)=yp0 fp(1:npar_loc,izp)=zp0 ! case ('random') if (lroot) print*, 'init_particles: Random particle positions' do k=1,npar_loc if (nxgrid/=1) call random_number_wrapper(fp(k,ixp)) if (nygrid/=1) call random_number_wrapper(fp(k,iyp)) if (nzgrid/=1) call random_number_wrapper(fp(k,izp)) enddo if (nxgrid/=1) & fp(1:npar_loc,ixp)=xyz0_loc(1)+fp(1:npar_loc,ixp)*Lxyz_loc(1) if (nygrid/=1) & fp(1:npar_loc,iyp)=xyz0_loc(2)+fp(1:npar_loc,iyp)*Lxyz_loc(2) if (nzgrid/=1) & fp(1:npar_loc,izp)=xyz0_loc(3)+fp(1:npar_loc,izp)*Lxyz_loc(3) ! case ('random-xy') if (lroot) print*, 'init_particles: Random particle positions' do k=1,npar_loc if (nxgrid/=1) call random_number_wrapper(fp(k,ixp)) if (nygrid/=1) call random_number_wrapper(fp(k,iyp)) enddo if (nxgrid/=1) & fp(1:npar_loc,ixp)=xyz0_loc(1)+fp(1:npar_loc,ixp)*Lxyz_loc(1) if (nygrid/=1) & fp(1:npar_loc,iyp)=xyz0_loc(2)+fp(1:npar_loc,iyp)*Lxyz_loc(2) fp(1:npar_loc,izp)=zp0 ! case ('random-xz') if (lroot) print*, 'init_particles: Random particle positions' do k=1,npar_loc if (nxgrid/=1) call random_number_wrapper(fp(k,ixp)) if (nzgrid/=1) call random_number_wrapper(fp(k,izp)) enddo if (nxgrid/=1) & fp(1:npar_loc,ixp)=xyz0_loc(1)+fp(1:npar_loc,ixp)*Lxyz_loc(1) if (nzgrid/=1) & fp(1:npar_loc,izp)=xyz0_loc(3)+fp(1:npar_loc,izp)*Lxyz_loc(3) fp(1:npar_loc,iyp)=yp0 ! case ('random-cylindrical','random-cyl') if (lroot) print*, 'init_particles: Random particle '//& 'cylindrical positions with power-law =',dustdensity_powerlaw ! do k=1,npar_loc ! ! Start the particles obeying a power law ! tmp=2-dustdensity_powerlaw call random_number_wrapper(rad_scl) rad_scl = rp_int**tmp + rad_scl*(rp_ext**tmp-rp_int**tmp) rad = rad_scl**(1./tmp) ! ! Random in azimuth ! call random_number_wrapper(phi) ! if (lcartesian_coords) then phi = 2*pi*phi if (nxgrid/=1) fp(k,ixp)=rad*cos(phi) if (nygrid/=1) fp(k,iyp)=rad*sin(phi) elseif (lcylindrical_coords) then phi = xyz0_par(2)+phi*Lxyz_par(2) if (nxgrid/=1) fp(k,ixp)=rad if (nygrid/=1) fp(k,iyp)=phi elseif (lspherical_coords) then call stop_it("init_particles: random-cylindrical not implemented "//& "for spherical coordinates") endif ! if (nzgrid/=1) call random_number_wrapper(fp(k,izp)) if (nzgrid/=1) & fp(k,izp)=xyz0_par(3)+fp(k,izp)*Lxyz_par(3) ! enddo ! case ('gaussian-z') if (lroot) print*, 'init_particles: Gaussian particle positions' do k=1,npar_loc if (nxgrid/=1) call random_number_wrapper(fp(k,ixp)) if (nygrid/=1) call random_number_wrapper(fp(k,iyp)) call random_number_wrapper(r) call random_number_wrapper(p) if (nprocz==2) then if (lfirst_proc_z) fp(k,izp)=-abs(zp0*sqrt(-2*alog(r))*cos(2*pi*p)) if (llast_proc_z) fp(k,izp)=abs(zp0*sqrt(-2*alog(r))*cos(2*pi*p)) else fp(k,izp)=zp0*sqrt(-2*alog(r))*cos(2*pi*p) endif enddo if (nxgrid/=1) & fp(1:npar_loc,ixp)=xyz0_loc(1)+fp(1:npar_loc,ixp)*Lxyz_loc(1) if (nygrid/=1) & fp(1:npar_loc,iyp)=xyz0_loc(2)+fp(1:npar_loc,iyp)*Lxyz_loc(2) ! case ('dragforce_equilibrium') ! ! Equilibrium between drag forces on dust and gas and other forces ! (from Nakagawa, Sekiya, & Hayashi 1986). ! if (lroot) then print*, 'init_particles: drag equilibrium' print*, 'init_particles: beta_glnrho_global=', beta_glnrho_global endif ! if (.not.(lcartesian_coords.and.(all(lequidist)))) & call fatal_error("init_particles","dragforce_equilibrium " //& "initial condition not implemented for polar or " //& "non-equidistant grids") ! Calculate average dust-to-gas ratio in box. if (ldensity_nolog) then eps = rhop_swarm*sum(f(l1:l2,m1:m2,n1:n2,inp))/ & sum(f(l1:l2,m1:m2,n1:n2,irho)) else eps = rhop_swarm*sum(f(l1:l2,m1:m2,n1:n2,inp))/ & sum(exp(f(l1:l2,m1:m2,n1:n2,ilnrho))) endif ! if (lroot) & print*, 'init_particles: average dust-to-gas ratio=', eps(1) ! Set gas velocity field. do m=m1,m2; do n=n1,n2 cs=sqrt(cs20) ! Take either global or local dust-to-gas ratio. if (.not. ldragforce_equi_global_eps) then if (ldensity_nolog) then eps = rhop_swarm*f(l1:l2,m,n,inp)/f(l1:l2,m,n,irho) else eps = rhop_swarm*f(l1:l2,m,n,inp)/exp(f(l1:l2,m,n,ilnrho)) endif endif ! f(l1:l2,m,n,iuy) = f(l1:l2,m,n,iuy) + & beta_glnrho_global(1)/(2*(1.0+eps))*cs ! enddo; enddo ! case default if (lroot) & print*, 'init_particles: No such such value for initxxp: ', & trim(initxxp(j)) call stop_it("") ! endselect ! enddo ! ! Particles are not allowed to be present in non-existing dimensions. ! This would give huge problems with interpolation later. ! if (nxgrid==1) fp(1:npar_loc,ixp)=x(l1) if (nygrid==1) fp(1:npar_loc,iyp)=y(m1) if (nzgrid==1) fp(1:npar_loc,izp)=z(n1) ! ! Redistribute particles among processors (now that positions are determined). ! call boundconds_particles(fp,ipar) ! ! Map particle positions on the grid. ! call map_nearest_grid(fp,ineargrid) call map_xxp_grid(f,fp,ineargrid) ! ! Sort particles so that they can be accessed contiguously in the memory. ! call sort_particles_imn(fp,ineargrid,ipar) ! endsubroutine init_particles !*********************************************************************** subroutine particles_dragforce_stiff(f,fp,ineargrid) ! ! 10-june-11/anders: dummy ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mpar_loc,mparray) :: fp integer, dimension (mpar_loc,3) :: ineargrid ! call keep_compiler_quiet(f) call keep_compiler_quiet(fp) call keep_compiler_quiet(ineargrid) ! endsubroutine particles_dragforce_stiff !*********************************************************************** subroutine pencil_criteria_particles() ! ! All pencils that the Particles module depends on are specified here. ! ! 20-04-06/anders: coded ! lpenc_diagnos(i_np)=.true. ! if (idiag_epsmin/=0.or.idiag_epsmax/=0) & lpenc_diagnos(i_epsp)=.true. ! endsubroutine pencil_criteria_particles !*********************************************************************** subroutine pencil_interdep_particles(lpencil_in) ! ! Interdependency among pencils provided by the Particles module ! is specified here. ! ! 15-feb-06/anders: coded ! logical, dimension(npencils) :: lpencil_in ! if (lpencil_in(i_rhop) .and. irhop==0) then lpencil_in(i_np)=.true. endif ! if (lpencil_in(i_epsp)) then lpencil_in(i_rho)=.true. lpencil_in(i_rhop)=.true. endif ! endsubroutine pencil_interdep_particles !*********************************************************************** subroutine calc_pencils_particles(f,p) ! ! Calculate particle pencils. ! ! 15-feb-06/anders: coded ! real, dimension (mx,my,mz,mfarray) :: f type (pencil_case) :: p ! np if (lpencil(i_np)) p%np=f(l1:l2,m,n,inp) ! rhop if (lpencil(i_rhop)) then if (irhop/=0) then p%rhop=f(l1:l2,m,n,irhop) else p%rhop=rhop_swarm*p%np endif endif ! epsp if (lpencil(i_epsp)) p%epsp=p%rhop/p%rho ! endsubroutine calc_pencils_particles !*********************************************************************** subroutine dxxp_dt(f,df,fp,dfp,ineargrid) ! ! Evolution of tracer particle position. ! ! 02-jan-05/anders: coded ! use Diagnostics, only: max_name ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df real, dimension (mpar_loc,mparray) :: fp real, dimension (mpar_loc,mpvar) :: dfp integer, dimension (mpar_loc,3) :: ineargrid ! intent (in) :: f, fp, ineargrid intent (inout) :: dfp, df ! logical :: lheader, lfirstcall=.true. ! lheader=lfirstcall .and. lroot ! ! Vertical gravity in the short friction time approximation. ! select case (gravz_profile) ! case ('zero') if (lheader) & print*, 'dxxp_dt_pencil: No gravity in z-direction.' ! case ('linear') if (lheader) & print*, 'dxxp_dt_pencil: Linear gravity field in z-direction.' dfp(1:npar_loc,izp)=dfp(1:npar_loc,izp) - & tausp*nu_epicycle2*fp(1:npar_loc,izp) ! endselect ! ! Diagnostic output ! if (ldiagnos) then if (idiag_xpm/=0) call sum_par_name(fp(1:npar_loc,ixp),idiag_xpm) if (idiag_ypm/=0) call sum_par_name(fp(1:npar_loc,iyp),idiag_ypm) if (idiag_zpm/=0) call sum_par_name(fp(1:npar_loc,izp),idiag_zpm) if (idiag_xp2m/=0) call sum_par_name(fp(1:npar_loc,ixp)**2,idiag_xp2m) if (idiag_yp2m/=0) call sum_par_name(fp(1:npar_loc,iyp)**2,idiag_yp2m) if (idiag_zp2m/=0) call sum_par_name(fp(1:npar_loc,izp)**2,idiag_zp2m) if (idiag_nparmax/=0) call max_name(npar_loc,idiag_nparmax) endif ! if (lfirstcall) lfirstcall=.false. ! call keep_compiler_quiet(f) call keep_compiler_quiet(df) call keep_compiler_quiet(ineargrid) ! endsubroutine dxxp_dt !*********************************************************************** subroutine dvvp_dt(f,df,p,fp,dfp,ineargrid) ! ! Evolution of dust particle velocity. ! ! 22-aug-05/anders: dummy ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df real, dimension (mpar_loc,mparray) :: fp real, dimension (mpar_loc,mpvar) :: dfp integer, dimension (mpar_loc,3) :: ineargrid type (pencil_case) :: p ! call keep_compiler_quiet(f) call keep_compiler_quiet(df) call keep_compiler_quiet(fp) call keep_compiler_quiet(dfp) call keep_compiler_quiet(ineargrid) ! endsubroutine dvvp_dt !*********************************************************************** subroutine dxxp_dt_pencil(f,df,fp,dfp,p,ineargrid) ! ! Evolution of tracer particle position (called from main pencil loop). ! ! 25-apr-06/anders: coded ! use Diagnostics ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df type (pencil_case) :: p real, dimension (mpar_loc,mparray) :: fp real, dimension (mpar_loc,mpvar) :: dfp integer, dimension (mpar_loc,3) :: ineargrid ! real, dimension (3) :: uu integer :: k logical :: lheader, lfirstcall=.true. ! intent (in) :: f, ineargrid intent (inout) :: df, dfp, fp ! ! Identify module and boundary conditions. ! lheader=lfirstcall .and. lroot if (lheader) then print*, 'dxxp_dt: Calculate dxxp_dt' print*, 'dxxp_dt: Particles boundary condition bcpx=', bcpx print*, 'dxxp_dt: Particles boundary condition bcpy=', bcpy print*, 'dxxp_dt: Particles boundary condition bcpz=', bcpz print*, 'dxxp_dt: Set rate of change of particle '// & 'position equal to gas velocity.' endif ! ! Interpolate gas velocity to position of particles. ! Then set particle velocity equal to the local gas velocity. ! if (npar_imn(imn)/=0) then do k=k1_imn(imn),k2_imn(imn) if (lparticlemesh_tsc) then if (ltrace_dust) then call interpolate_quadratic_spline(f,iudx(1),iudz(1), & fp(k,ixp:izp),uu,ineargrid(k,:),0,ipar(k)) else call interpolate_quadratic_spline(f,iux,iuz, & fp(k,ixp:izp),uu,ineargrid(k,:),0,ipar(k)) endif else if (ltrace_dust) then call interpolate_linear(f,iudx(1),iudz(1), & fp(k,ixp:izp),uu,ineargrid(k,:),0,ipar(k)) else call interpolate_linear(f,iux,iuz, & fp(k,ixp:izp),uu,ineargrid(k,:),0,ipar(k)) endif endif ! ! Advance of particle position ! if (lcartesian_coords) then if (nxgrid/=1) dfp(k,ixp) = dfp(k,ixp) + uu(1) ; fp(k,ivpx) = uu(1) if (nygrid/=1) dfp(k,iyp) = dfp(k,iyp) + uu(2) ; fp(k,ivpy) = uu(2) if (nzgrid/=1) dfp(k,izp) = dfp(k,izp) + uu(3) ; fp(k,ivpz) = uu(3) elseif (lcylindrical_coords) then if (nxgrid/=1) dfp(k,ixp) = dfp(k,ixp) + uu(1) if (nygrid/=1) dfp(k,iyp) = dfp(k,iyp) + uu(2)/max(fp(k,ixp),tini) if (nzgrid/=1) dfp(k,izp) = dfp(k,izp) + uu(3) elseif (lspherical_coords) then if (nxgrid/=1) dfp(k,ixp) = dfp(k,ixp) + uu(1) if (nygrid/=1) dfp(k,iyp) = dfp(k,iyp) + uu(2)/max(fp(k,ixp),tini) if (nzgrid/=1) dfp(k,izp) = dfp(k,izp) + uu(3)/(max(fp(k,ixp),tini)*sin(fp(k,iyp))) endif ! ! With shear there is an extra term due to the background shear flow. ! if (lshear.and.nygrid/=1)& dfp(k,iyp) = dfp(k,iyp) - qshear*Omega*fp(k,ixp) ! enddo endif ! if (ldiagnos) then if (idiag_npmax/=0) call max_mn_name(p%np,idiag_npmax) if (idiag_npmin/=0) call max_mn_name(-p%np,idiag_npmin,lneg=.true.) if (idiag_epsmax/=0) call max_mn_name(p%epsp,idiag_epsmax) if (idiag_epsmin/=0) call max_mn_name(-p%epsp,idiag_epsmin,lneg=.true.) endif ! if (l1davgfirst) then call yzsum_mn_name_x(p%np,idiag_npmx) call yzsum_mn_name_x(p%rhop,idiag_rhopmx) call yzsum_mn_name_x(p%epsp,idiag_epspmx) call xysum_mn_name_z(p%np,idiag_npmz) call xysum_mn_name_z(p%rhop,idiag_rhopmz) call xysum_mn_name_z(p%epsp,idiag_epspmz) endif ! if (lfirstcall) lfirstcall=.false. ! call keep_compiler_quiet(df) ! endsubroutine dxxp_dt_pencil !*********************************************************************** subroutine dvvp_dt_pencil(f,df,fp,dfp,p,ineargrid) ! ! Evolution of dust particle velocity (called from main pencil loop). ! ! 20-apr-06/anders: dummy ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df real, dimension (mpar_loc,mparray) :: fp real, dimension (mpar_loc,mpvar) :: dfp type (pencil_case) :: p integer, dimension (mpar_loc,3) :: ineargrid ! call keep_compiler_quiet(f) call keep_compiler_quiet(df) call keep_compiler_quiet(fp) call keep_compiler_quiet(dfp) call keep_compiler_quiet(p) call keep_compiler_quiet(ineargrid) ! endsubroutine dvvp_dt_pencil !*********************************************************************** subroutine dxxp_dt_blocks(f,df,fp,dfp,ineargrid) ! ! Evolution of particle position in blocks. ! ! 29-nov-09/anders: dummy ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df real, dimension (mpar_loc,mparray) :: fp real, dimension (mpar_loc,mpvar) :: dfp integer, dimension (mpar_loc,3) :: ineargrid ! call keep_compiler_quiet(f) call keep_compiler_quiet(df) call keep_compiler_quiet(fp) call keep_compiler_quiet(dfp) call keep_compiler_quiet(ineargrid) ! endsubroutine dxxp_dt_blocks !*********************************************************************** subroutine dvvp_dt_blocks(f,df,fp,dfp,ineargrid) ! ! Evolution of particle velocity in blocks. ! ! 29-nov-09/anders: dummy ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df real, dimension (mpar_loc,mparray) :: fp real, dimension (mpar_loc,mpvar) :: dfp integer, dimension (mpar_loc,3) :: ineargrid ! call keep_compiler_quiet(f) call keep_compiler_quiet(df) call keep_compiler_quiet(fp) call keep_compiler_quiet(dfp) call keep_compiler_quiet(ineargrid) ! endsubroutine dvvp_dt_blocks !*********************************************************************** subroutine remove_particles_sink_simple(f,fp,dfp,ineargrid) ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mpar_loc,mparray) :: fp real, dimension (mpar_loc,mpvar) :: dfp integer, dimension (mpar_loc,3) :: ineargrid ! call keep_compiler_quiet(f) call keep_compiler_quiet(fp) call keep_compiler_quiet(dfp) call keep_compiler_quiet(ineargrid) ! endsubroutine remove_particles_sink_simple !*********************************************************************** subroutine create_particles_sink_simple(f,fp,dfp,ineargrid) ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mpar_loc,mparray) :: fp real, dimension (mpar_loc,mpvar) :: dfp integer, dimension (mpar_loc,3) :: ineargrid ! call keep_compiler_quiet(f) call keep_compiler_quiet(fp) call keep_compiler_quiet(dfp) call keep_compiler_quiet(ineargrid) ! endsubroutine create_particles_sink_simple !*********************************************************************** subroutine read_particles_init_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=particles_init_pars, IOSTAT=iostat) ! endsubroutine read_particles_init_pars !*********************************************************************** subroutine write_particles_init_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=particles_init_pars) ! endsubroutine write_particles_init_pars !*********************************************************************** subroutine read_particles_run_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=particles_run_pars, IOSTAT=iostat) ! endsubroutine read_particles_run_pars !*********************************************************************** subroutine write_particles_run_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=particles_run_pars) ! endsubroutine write_particles_run_pars !*********************************************************************** subroutine powersnap_particles(f) ! ! Calculate power spectra of particle variables. ! ! 01-jan-06/anders: dummy ! real, dimension (mx,my,mz,mfarray) :: f ! call keep_compiler_quiet(f) ! endsubroutine powersnap_particles !*********************************************************************** subroutine insert_lost_particles(f,fp,ineargrid) ! ! 14-oct-12/dhruba: dummy ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mpar_loc,mparray) :: fp integer, dimension (mpar_loc,3) :: ineargrid ! intent (inout) :: fp,ineargrid ! call keep_compiler_quiet(f) call keep_compiler_quiet(fp) call keep_compiler_quiet(ineargrid) ! endsubroutine insert_lost_particles !*********************************************************************** subroutine insert_particles(f,fp,ineargrid) ! ! Insert particles continuously (when linsert_particles_continuously == T), ! i.e. in each timestep. If number of particles to be inserted are less ! than unity, accumulate number over several timesteps until the integer value ! is larger than one. Keep the remainder and accumulate this to the next insert. ! ! Works only for particles_dust - add neccessary variable ! declarations in particles_tracers to make it work here. ! ! use General, only: random_number_wrapper ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mpar_loc,mparray) :: fp integer, dimension (mpar_loc,3) :: ineargrid logical :: linsertmore=.true. real :: avg_n_insert,remaining_particles ! integer :: j, k, n_insert, npar_loc_old, iii ! intent (inout) :: fp,ineargrid ! ! Stop call to this routine when maximum number of particles is reached! ! Since root inserts all new particles, make sure ! npar_total + n_insert < mpar ! so that a processor can not exceed its maximum number of particles. ! if (lroot) then avg_n_insert=particles_insert_rate*dt n_insert=int(avg_n_insert + remaining_particles) ! Remaining particles saved for subsequent timestep: remaining_particles=avg_n_insert + remaining_particles - n_insert if ((n_insert+npar_total <= mpar_loc) & .and. (t