! $Id$
!
! This module writes information about the local state of the gas at
! the positions of a selected number of 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_stalker=.true.
!
! MSCRATCH CONTRIBUTION 1
!
!***************************************************************
module Particles_stalker
!
use Cdata
use Messages
use Particles_cdata
use Particles_map
use Particles_mpicomm
!
implicit none
!
include 'particles_stalker.h'
!
real :: dstalk=0.1, tstalk=0.0
integer :: iscratch=0, nout=0, nvar_stalk=0
logical :: linterpolate_cic=.false., linterpolate_tsc=.true.
logical :: lstalk_xx=.true., lstalk_vv=.true.
logical :: lstalk_uu=.true., lstalk_guu=.false.
logical :: lstalk_rho=.true., lstalk_grho=.false.
logical :: lstalk_bb=.true., lstalk_ap=.true.
logical :: lstalk_npswarm=.false., lstalk_rhopswarm=.true., lstalk_potself=.true.
logical :: lstalk_aps=.true.
logical :: lstalk_sink_particles=.false.
logical :: lstalk_relvel=.false.
logical :: lstalk_gTT=.false.
!
real, dimension (:), allocatable :: xp, yp, zp, vpx, vpy, vpz, ux, uy, uz
real, dimension (:), allocatable :: rho, drhodx, drhody, drhodz
real, dimension (:), allocatable :: duxdx, duxdy, duxdz, duydx, duydy, duydz, duzdx, duzdy, duzdz
real, dimension (:), allocatable :: gTTx, gTTy, gTTz
real, dimension (:), allocatable :: bx, by, bz, ap, npswarm, rhopswarm
real, dimension (:), allocatable :: potself, aps, relvel
!
namelist /particles_stalker_init_pars/ &
dstalk, linterpolate_cic, linterpolate_tsc, &
lstalk_xx, lstalk_vv, lstalk_uu, lstalk_guu, lstalk_rho, lstalk_grho, &
lstalk_bb, lstalk_ap, lstalk_npswarm, lstalk_rhopswarm, lstalk_potself, &
lstalk_aps, lstalk_sink_particles, lstalk_relvel,lstalk_gTT
!
namelist /particles_stalker_run_pars/ &
dstalk, linterpolate_cic, linterpolate_tsc, lstalk_sink_particles
!
contains
!***********************************************************************
subroutine initialize_particles_stalker(f)
!
! Perform any post-parameter-read initialization i.e. calculate derived
! parameters.
!
! 13-nov-07/anders: coded
!
use General, only: keep_compiler_quiet
use FArrayManager
!
real, dimension (mx,my,mz,mfarray) :: f
!
! Stop if no particles are stalked.
!
if (npar_stalk==0) then
if (lroot) print*, 'initialize_particles_stalker: npar_stalk is zero - set it in cparam.local'
call fatal_error('initialize_particles_stalker','')
endif
!
! Turn off stalking if physics not selected.
!
if (ivpx==0) lstalk_vv=.false.
if (iap==0) lstalk_ap=.false.
if (inpswarm==0) lstalk_npswarm=.false.
if (irhopswarm==0) lstalk_rhopswarm=.false.
if (iaps==0) lstalk_aps=.false.
if (iuu==0) lstalk_uu=.false.
if (iuu==0) lstalk_guu=.false.
if (ilnrho==0) lstalk_rho=.false.
if (ilnrho==0) lstalk_grho=.false.
if (iaa==0) lstalk_bb=.false.
if (ipotself==0) lstalk_potself=.false.
if (iuu==0 .or. ivpx==0) lstalk_relvel=.false.
!
! Need scratch slot in f array to interpolate derived variables.
!
if (lstalk_guu .or. lstalk_grho .or. lstalk_bb) &
call farray_acquire_scratch_area('scratch',iscratch)
!
! Count the number of variables to be stalked.
!
nvar_stalk=0
if (lstalk_xx) nvar_stalk=nvar_stalk+3
if (lstalk_vv) nvar_stalk=nvar_stalk+3
if (lstalk_ap) nvar_stalk=nvar_stalk+1
if (lstalk_npswarm) nvar_stalk=nvar_stalk+1
if (lstalk_rhopswarm) nvar_stalk=nvar_stalk+1
if (lstalk_aps) nvar_stalk=nvar_stalk+1
if (lstalk_uu) nvar_stalk=nvar_stalk+3
if (lstalk_guu) nvar_stalk=nvar_stalk+9
if (lstalk_rho) nvar_stalk=nvar_stalk+1
if (lstalk_grho) nvar_stalk=nvar_stalk+3
if (lstalk_bb) nvar_stalk=nvar_stalk+3
if (lstalk_potself) nvar_stalk=nvar_stalk+1
if (lstalk_relvel) nvar_stalk=nvar_stalk+1
if (lstalk_gTT) nvar_stalk=nvar_stalk+3
!
! Reserve space for particle arrays.
!
if (lstalk_xx) allocate (xp(npar_stalk), yp(npar_stalk), zp(npar_stalk))
if (lstalk_vv) allocate (vpx(npar_stalk), vpy(npar_stalk), vpz(npar_stalk))
if (lstalk_ap) allocate (ap(npar_stalk))
if (lstalk_npswarm) allocate (npswarm(npar_stalk))
if (lstalk_rhopswarm) allocate (rhopswarm(npar_stalk))
if (lstalk_aps) allocate (aps(npar_stalk))
if (lstalk_uu) allocate (ux(npar_stalk), uy(npar_stalk), uz(npar_stalk))
if (lstalk_guu) allocate (duxdx(npar_stalk), duxdy(npar_stalk), duxdz(npar_stalk), duydx(npar_stalk), duydy(npar_stalk), &
duydz(npar_stalk), duzdx(npar_stalk), duzdy(npar_stalk), duzdz(npar_stalk))
if (lstalk_rho) allocate (rho(npar_stalk))
if (lstalk_grho) allocate (drhodx(npar_stalk), drhody(npar_stalk), drhodz(npar_stalk))
if (lstalk_bb) allocate (bx(npar_stalk), by(npar_stalk), bz(npar_stalk))
if (lstalk_potself) allocate (potself(npar_stalk))
if (lstalk_relvel) allocate (relvel(npar_stalk))
if (lstalk_gTT) allocate (gTTx(npar_stalk), gTTy(npar_stalk), gTTz(npar_stalk))
!
! Write information on which variables are stalked to file.
!
if (lroot) then
open(1,file=trim(datadir)//'/particles_stalker_header.dat', &
status='unknown')
if (lstalk_xx) write(1,'(A)',advance='no') 'xp,yp,zp,'
if (lstalk_vv) write(1,'(A)',advance='no') 'vpx,vpy,vpz,'
if (lstalk_ap) write(1,'(A)',advance='no') 'ap,'
if (lstalk_npswarm) write(1,'(A)',advance='no') 'npswarm,'
if (lstalk_rhopswarm) write(1,'(A)',advance='no') 'rhopswarm,'
if (lstalk_aps) write(1,'(A)',advance='no') 'aps,'
if (lstalk_uu) write(1,'(A)',advance='no') 'ux,uy,uz,'
if (lstalk_guu) write(1,'(A)',advance='no') 'duxdx,duxdy,duxdz,'
if (lstalk_guu) write(1,'(A)',advance='no') 'duydx,duydy,duydz,'
if (lstalk_guu) write(1,'(A)',advance='no') 'duzdx,duzdy,duzdz,'
if (lstalk_rho) write(1,'(A)',advance='no') 'rho,'
if (lstalk_grho) write(1,'(A)',advance='no') 'drhodx,drhody,drhodz,'
if (lstalk_bb) write(1,'(A)',advance='no') 'bx,by,bz,'
if (lstalk_potself) write(1,'(A)',advance='no') 'potself,'
if (lstalk_relvel) write(1,'(A)',advance='no') 'relvel,'
if (lstalk_gTT) write(1,'(A)',advance='no') 'dTTdx,dTTdy,dTTdz,'
close (1)
endif
!
! Read time of next stalking from file.
!
if (lrun) then
open(1,file=trim(datadir)//'/tstalk.dat',form='formatted', &
status='unknown')
read(1,*) tstalk, nout
close(1)
else
tstalk=0.0
if (lroot) then
open(1,file=trim(datadir)//'/tstalk.dat',form='formatted')
write(1,'(e14.7,i8)') tstalk, nout
close(1)
endif
endif
!
call keep_compiler_quiet(f)
!
endsubroutine initialize_particles_stalker
!***********************************************************************
subroutine finalize_particles_stalker
!
! Finalize module.
!
! 01-May-2019/PABourdin: coded
!
if (lstalk_xx) deallocate (xp, yp, zp)
if (lstalk_vv) deallocate (vpx, vpy, vpz)
if (lstalk_ap) deallocate (ap)
if (lstalk_npswarm) deallocate (npswarm)
if (lstalk_rhopswarm) deallocate (rhopswarm)
if (lstalk_aps) deallocate (aps)
if (lstalk_uu) deallocate (ux, uy, uz)
if (lstalk_guu) deallocate (duxdx, duxdy, duxdz, duydx, duydy, duydz, duzdx, duzdy, duzdz)
if (lstalk_rho) deallocate (rho)
if (lstalk_grho) deallocate (drhodx, drhody, drhodz)
if (lstalk_bb) deallocate (bx, by, bz)
if (lstalk_potself) deallocate (potself)
if (lstalk_relvel) deallocate (relvel)
if (lstalk_gTT) deallocate (gTTx, gTTy, gTTz)
!
endsubroutine finalize_particles_stalker
!***********************************************************************
subroutine particles_stalker_sub(f,fp,ineargrid)
!
! Find local state of the gas at the position of a limited number of
! particles. This information is saved to a file and can be read in
! with IDL.
!
! 13-nov-07/anders: coded
!
use IO, only: output_stalker_init, output_stalker, output_part_finalize
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mpar_loc,mparray) :: fp
integer, dimension (mpar_loc,3) :: ineargrid
real :: t_sp
!
integer, dimension (npar_stalk) :: k_stalk
integer :: i, k, npar_stalk_loc, ivalue
!
! Only stalk particles every dstalk.
!
if (t<tstalk) then
! Do nothing.
else
t_sp=t
!
! Find out how many stalked particles are at the local processor.
!
npar_stalk_loc=0
if (lstalk_sink_particles) then
do k=1,npar_loc
if (fp(k,iaps)>0.0) then
npar_stalk_loc=npar_stalk_loc+1
if (npar_stalk_loc>npar_stalk) then
print*, 'particles_stalker_sub: too many sink particles '// &
'are stalked'
print*, 'iproc, it, itsub, npar_stalk, npar_stalk_loc=', &
it, itsub, npar_stalk, npar_stalk_loc
call fatal_error_local('particles_stalker_sub','')
else
k_stalk(npar_stalk_loc)=k
endif
endif
enddo
call fatal_error_local_collect()
else
do k=1,npar_loc
if (ipar(k)<=npar_stalk) then
npar_stalk_loc=npar_stalk_loc+1
k_stalk(npar_stalk_loc)=k
endif
enddo
endif
!
! Gather environment information for each stalked particle, starting with the
! position.
!
if (lstalk_xx) then
do i=1,npar_stalk_loc
xp(i)=fp(k_stalk(i),ixp)
yp(i)=fp(k_stalk(i),iyp)
zp(i)=fp(k_stalk(i),izp)
enddo
endif
!
! Velocity vector (only relevant for inertial particles, not for tracers).
!
if (lstalk_vv) then
do i=1,npar_stalk_loc
vpx(i)=fp(k_stalk(i),ivpx)
vpy(i)=fp(k_stalk(i),ivpy)
vpz(i)=fp(k_stalk(i),ivpz)
enddo
endif
!
! Particle radius.
!
if (lstalk_ap) then
do i=1,npar_stalk_loc
ap(i)=fp(k_stalk(i),iap)
enddo
endif
!
! Particle number density.
!
if (lstalk_npswarm) then
do i=1,npar_stalk_loc
npswarm(i)=fp(k_stalk(i),inpswarm)
enddo
endif
!
! Particle mass density.
!
if (lstalk_rhopswarm) then
do i=1,npar_stalk_loc
rhopswarm(i)=fp(k_stalk(i),irhopswarm)
enddo
endif
!
! Sink particle radius.
!
if (lstalk_aps) then
do i=1,npar_stalk_loc
aps(i)=fp(k_stalk(i),iaps)
enddo
endif
!
! Local gas velocity.
!
if (lstalk_uu) then
call stalk_variable(f,fp,k_stalk,npar_stalk_loc,ineargrid,iux,ux)
call stalk_variable(f,fp,k_stalk,npar_stalk_loc,ineargrid,iuy,uy)
call stalk_variable(f,fp,k_stalk,npar_stalk_loc,ineargrid,iuz,uz)
endif
!
! Local gradient of gas velocity.
!
if (lstalk_guu) then
call stalk_gradient(f,fp,k_stalk,npar_stalk_loc,ineargrid,iux,1,duxdx)
call stalk_gradient(f,fp,k_stalk,npar_stalk_loc,ineargrid,iux,2,duxdy)
call stalk_gradient(f,fp,k_stalk,npar_stalk_loc,ineargrid,iux,3,duxdz)
call stalk_gradient(f,fp,k_stalk,npar_stalk_loc,ineargrid,iuy,1,duydx)
call stalk_gradient(f,fp,k_stalk,npar_stalk_loc,ineargrid,iuy,2,duydy)
call stalk_gradient(f,fp,k_stalk,npar_stalk_loc,ineargrid,iuy,3,duydz)
call stalk_gradient(f,fp,k_stalk,npar_stalk_loc,ineargrid,iuz,1,duzdx)
call stalk_gradient(f,fp,k_stalk,npar_stalk_loc,ineargrid,iuz,2,duzdy)
call stalk_gradient(f,fp,k_stalk,npar_stalk_loc,ineargrid,iuz,3,duzdz)
endif
!
! Local gas density (logarithmic or not logarithmic).
!
if (lstalk_rho) then
call stalk_variable(f,fp,k_stalk,npar_stalk_loc,ineargrid,ilnrho,rho)
endif
!
! Local density gradient.
!
if (lstalk_grho) then
call stalk_gradient(f,fp,k_stalk,npar_stalk_loc,ineargrid,ilnrho,1, &
drhodx)
call stalk_gradient(f,fp,k_stalk,npar_stalk_loc,ineargrid,ilnrho,2, &
drhody)
call stalk_gradient(f,fp,k_stalk,npar_stalk_loc,ineargrid,ilnrho,3, &
drhodz)
endif
!
! Local magnetic field.
!
if (lstalk_bb) then
call stalk_magnetic(f,fp,k_stalk,npar_stalk_loc,ineargrid,bx,by,bz)
endif
!
! Local gravitational potential.
!
if (lstalk_potself) then
call stalk_variable(f,fp,k_stalk,npar_stalk_loc,ineargrid,ipotself, &
potself)
endif
!
! Relative velocity of gas and particle
!
if (lstalk_relvel) then
allocate (relvel(npar_stalk))
do i=1,npar_stalk_loc
vpx(i)=fp(k_stalk(i),ivpx)
vpy(i)=fp(k_stalk(i),ivpy)
vpz(i)=fp(k_stalk(i),ivpz)
enddo
call stalk_variable(f,fp,k_stalk,npar_stalk_loc,ineargrid,iux,ux)
call stalk_variable(f,fp,k_stalk,npar_stalk_loc,ineargrid,iuy,uy)
call stalk_variable(f,fp,k_stalk,npar_stalk_loc,ineargrid,iuz,uz)
do i=1,npar_stalk_loc
relvel(i)=sqrt(((vpx(i)-ux(i))**2)+((vpy(i)-uy(i))**2)+((vpz(i)-uz(i))**2))
enddo
endif
!
! Write information to a file
!
call output_stalker_init (nvar_stalk, npar_stalk_loc, nout, ipar(k_stalk(1:npar_stalk_loc)))
if (lstalk_xx) then
call output_stalker ('xp', npar_stalk, npar_stalk_loc, xp)
call output_stalker ('yp', npar_stalk, npar_stalk_loc, yp)
call output_stalker ('zp', npar_stalk, npar_stalk_loc, zp)
endif
if (lstalk_vv) then
call output_stalker ('vpx', npar_stalk, npar_stalk_loc, vpx)
call output_stalker ('vpy', npar_stalk, npar_stalk_loc, vpy)
call output_stalker ('vpz', npar_stalk, npar_stalk_loc, vpz)
endif
if (lstalk_ap) call output_stalker ('ap', npar_stalk, npar_stalk_loc, ap)
if (lstalk_npswarm) call output_stalker ('npswarm', npar_stalk, npar_stalk_loc, npswarm)
if (lstalk_rhopswarm) call output_stalker ('rhopswarm', npar_stalk, npar_stalk_loc, rhopswarm)
if (lstalk_aps) call output_stalker ('aps', npar_stalk, npar_stalk_loc, aps)
if (lstalk_uu) then
call output_stalker ('ux', npar_stalk, npar_stalk_loc, ux)
call output_stalker ('uy', npar_stalk, npar_stalk_loc, uy)
call output_stalker ('uz', npar_stalk, npar_stalk_loc, uz)
endif
if (lstalk_guu) then
call output_stalker ('duxdx', npar_stalk, npar_stalk_loc, duxdx)
call output_stalker ('duxdy', npar_stalk, npar_stalk_loc, duxdy)
call output_stalker ('duxdz', npar_stalk, npar_stalk_loc, duxdz)
call output_stalker ('duydx', npar_stalk, npar_stalk_loc, duydx)
call output_stalker ('duydy', npar_stalk, npar_stalk_loc, duydy)
call output_stalker ('duydz', npar_stalk, npar_stalk_loc, duydz)
call output_stalker ('duzdx', npar_stalk, npar_stalk_loc, duzdx)
call output_stalker ('duzdy', npar_stalk, npar_stalk_loc, duzdy)
call output_stalker ('duzdz', npar_stalk, npar_stalk_loc, duzdz)
endif
if (lstalk_rho) call output_stalker ('rho', npar_stalk, npar_stalk_loc, rho)
if (lstalk_grho) then
call output_stalker ('drhodx', npar_stalk, npar_stalk_loc, drhodx)
call output_stalker ('drhody', npar_stalk, npar_stalk_loc, drhody)
call output_stalker ('drhodz', npar_stalk, npar_stalk_loc, drhodz)
endif
if (lstalk_bb) then
call output_stalker ('bx', npar_stalk, npar_stalk_loc, bx)
call output_stalker ('by', npar_stalk, npar_stalk_loc, by)
call output_stalker ('bz', npar_stalk, npar_stalk_loc, bz)
endif
if (lstalk_potself) call output_stalker ('potself', npar_stalk, npar_stalk_loc, potself)
if (lstalk_relvel) call output_stalker ('relvel', npar_stalk, npar_stalk_loc, relvel)
if (lstalk_gTT) then
call output_stalker ('gTTx', npar_stalk, npar_stalk_loc, gTTx)
call output_stalker ('gTTy', npar_stalk, npar_stalk_loc, gTTy)
call output_stalker ('gTTz', npar_stalk, npar_stalk_loc, gTTz)
endif
call output_part_finalize
!
! Next stalking time is dstalk later.
!
do while (tstalk<=t_sp)
tstalk=tstalk+dstalk
enddo
nout=nout+1
if (lroot) then
open(1,file=trim(datadir)//'/tstalk.dat',form='formatted')
write(1,'(e14.7,i8)') tstalk, nout
close(1)
endif
endif
!
endsubroutine particles_stalker_sub
!***********************************************************************
subroutine stalk_variable(f,fp,k_stalk,npar_stalk_loc,ineargrid,ivar,value)
!
! Find local state of the gas at the position of a limited number of
! particles. This subroutines take care of interpolation all simple
! physical variables, i.e. the ones that are directly present in the f
! array.
!
! 13-nov-07/anders: coded
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mpar_loc,mparray) :: fp
integer, dimension (npar_stalk) :: k_stalk
integer :: npar_stalk_loc
integer, dimension (mpar_loc,3) :: ineargrid
integer :: ivar
real, dimension (npar_stalk) :: value
!
real, dimension (1) :: value_loc
integer :: i, k, ix0, iy0, iz0, iblock
!
do i=1,npar_stalk_loc
k=k_stalk(i)
ix0=ineargrid(k_stalk(i),1)
iy0=ineargrid(k_stalk(i),2)
iz0=ineargrid(k_stalk(i),3)
!
if (lparticles_blocks) then
iblock=inearblock(k)
else
iblock=0
endif
!
! Interpolation is either zeroth, first or second order spline interpolation.
!
if (lparticlemesh_cic) then
call interpolate_linear( &
f,ivar,ivar,fp(k,ixp:izp),value_loc,ineargrid(k,:),iblock,ipar(k))
elseif (lparticlemesh_tsc) then
call interpolate_quadratic_spline( &
f,ivar,ivar,fp(k,ixp:izp),value_loc,ineargrid(k,:),iblock,ipar(k))
else
value_loc=(/f(ix0,iy0,iz0,ivar)/)
endif
value(i)=value_loc(1)
enddo
!
endsubroutine stalk_variable
!***********************************************************************
subroutine stalk_gradient(f,fp,k_stalk,npar_stalk_loc,ineargrid,ivar,ider,value)
!
! Take simple derivative of physical variable in f and store the result
! in scratch slot in f. Then find the local state of the gas using the usual
! subroutine 'stalk_variable'.
!
! 13-nov-07/anders: coded
!
use Boundcond, only: bc_per_x, bc_per_y, bc_per_z
use Deriv, only: der
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mpar_loc,mparray) :: fp
integer, dimension (npar_stalk) :: k_stalk
integer :: npar_stalk_loc
integer, dimension (mpar_loc,3) :: ineargrid
integer :: ivar, ider
real, dimension (npar_stalk) :: value
!
real, dimension (nx) :: der_pencil
!
! Calculate derivative of variable number ivar with respect to
! direction number ider.
!
do n=n1,n2; do m=m1,m2
call der(f(:,:,:,ivar),der_pencil,ider)
f(l1:l2,m,n,iscratch)=der_pencil
enddo; enddo
call bc_per_x(f,'top',iscratch); call bc_per_x(f,'bot',iscratch)
call bc_per_y(f,'top',iscratch); call bc_per_y(f,'bot',iscratch)
call bc_per_z(f,'top',iscratch); call bc_per_z(f,'bot',iscratch)
!
! Now that the derivative is stored in the f array, we can use the usual
! subroutine to find the local state of the gas at the positions of the
! particles.
!
call stalk_variable(f,fp,k_stalk,npar_stalk_loc,ineargrid,iscratch,value)
!
endsubroutine stalk_gradient
!***********************************************************************
subroutine stalk_magnetic(f,fp,k_stalk,npar_stalk_loc,ineargrid,bx,by,bz)
!
! Calculate magnetic field B from vector potential A and interpolate
! to position of stalked particles.
!
! 14-nov-07/anders: coded
!
use Boundcond, only: bc_per_x, bc_per_y, bc_per_z
use Deriv, only: der
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mpar_loc,mparray) :: fp
integer, dimension (npar_stalk) :: k_stalk
integer :: npar_stalk_loc
integer, dimension (mpar_loc,3) :: ineargrid
real, dimension (npar_stalk) :: bx, by, bz
!
real, dimension (nx) :: der1_pencil, der2_pencil
!
do n=n1,n2; do m=m1,m2
call der(f,iaa-1+3,der1_pencil,2)
call der(f,iaa-1+2,der2_pencil,3)
f(l1:l2,m,n,iscratch)=der1_pencil-der2_pencil
enddo; enddo
call bc_per_x(f,'top',iscratch); call bc_per_x(f,'bot',iscratch)
call bc_per_y(f,'top',iscratch); call bc_per_y(f,'bot',iscratch)
call bc_per_z(f,'top',iscratch); call bc_per_z(f,'bot',iscratch)
call stalk_variable(f,fp,k_stalk,npar_stalk_loc,ineargrid,iscratch,bx)
!
do n=n1,n2; do m=m1,m2
call der(f,iaa-1+1,der1_pencil,3)
call der(f,iaa-1+3,der2_pencil,1)
f(l1:l2,m,n,iscratch)=der1_pencil-der2_pencil
enddo; enddo
call bc_per_x(f,'top',iscratch); call bc_per_x(f,'bot',iscratch)
call bc_per_y(f,'top',iscratch); call bc_per_y(f,'bot',iscratch)
call bc_per_z(f,'top',iscratch); call bc_per_z(f,'bot',iscratch)
call stalk_variable(f,fp,k_stalk,npar_stalk_loc,ineargrid,iscratch,by)
!
do n=n1,n2; do m=m1,m2
call der(f,iaa-1+2,der1_pencil,1)
call der(f,iaa-1+1,der2_pencil,2)
f(l1:l2,m,n,iscratch)=der1_pencil-der2_pencil
enddo; enddo
call bc_per_x(f,'top',iscratch); call bc_per_x(f,'bot',iscratch)
call bc_per_y(f,'top',iscratch); call bc_per_y(f,'bot',iscratch)
call bc_per_z(f,'top',iscratch); call bc_per_z(f,'bot',iscratch)
call stalk_variable(f,fp,k_stalk,npar_stalk_loc,ineargrid,iscratch,bz)
!
endsubroutine stalk_magnetic
!***********************************************************************
! subroutine stalk_rep(fp,npar_stalk_loc,rep)
! real, dimension (mpar_loc,mparray) :: fp
! real, dimension (npar_stalk) :: rep,nu
! integer :: npar_stalk_loc
! integer :: k
!
! do k=1,npar_stalk
! rep(k)=2.0*fp(k,iap)*sqrt(sum((interp_uu(k,:)-fp(k,ivpx:ivpz))**2))/nu(k)
! enddo
! endsubroutine stalk_rep
!***********************************************************************
subroutine read_pstalker_init_pars(iostat)
!
use File_io, only: parallel_unit
!
integer, intent(out) :: iostat
!
read(parallel_unit, NML=particles_stalker_init_pars, IOSTAT=iostat)
!
endsubroutine read_pstalker_init_pars
!***********************************************************************
subroutine write_pstalker_init_pars(unit)
!
integer, intent(in) :: unit
!
write(unit, NML=particles_stalker_init_pars)
!
endsubroutine write_pstalker_init_pars
!***********************************************************************
subroutine read_pstalker_run_pars(iostat)
!
use File_io, only: parallel_unit
!
integer, intent(out) :: iostat
!
read(parallel_unit, NML=particles_stalker_run_pars, IOSTAT=iostat)
!
endsubroutine read_pstalker_run_pars
!***********************************************************************
subroutine write_pstalker_run_pars(unit)
!
integer, intent(in) :: unit
!
write(unit, NML=particles_stalker_run_pars)
!
endsubroutine write_pstalker_run_pars
!***********************************************************************
endmodule Particles_stalker