! $Id$
!
! This module takes care of everything related to inertial particles.
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
!
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! MPVAR CONTRIBUTION 6
! MAUX CONTRIBUTION 2
! CPARAM logical, parameter :: lparticles=.true.
! CPARAM character (len=20), parameter :: particles_module="lagrangian"
!
! PENCILS PROVIDED np; rhop; vol; peh
! PENCILS PROVIDED np_rad(ndustrad); npvz(ndustrad); npvz2(ndustrad);
! PENCILS PROVIDED npuz(ndustrad); sherwood
! PENCILS PROVIDED epsp; grhop(3)
!
!***************************************************************
module Particles
!
use Cdata
use Cparam
use General, only: keep_compiler_quiet
use Messages
use Particles_cdata
use Particles_map
use Particles_mpicomm
use Particles_sub
use Particles_radius
!
implicit none
!
include 'particles.h'
include 'particles_common.h'
!
complex, dimension (7) :: coeff=(0.0,0.0)
real, target, dimension (npar_species) :: tausp_species=0.0
real, target, dimension (npar_species) :: tausp1_species=0.0
real, target, dimension (npar_species) :: rpbeta_species=0.0
real, dimension (3) :: temp_grad0=(/0.0,0.0,0.0/)
real, dimension (3) :: pos_sphere=(/0.0,0.0,0.0/)
real, dimension (3) :: pos_ellipsoid=(/0.0,0.0,0.0/)
real :: xp0=0.0, yp0=0.0, zp0=0.0, vpx0=0.0, vpy0=0.0, vpz0=0.0
real :: xp1=0.0, yp1=0.0, zp1=0.0, vpx1=0.0, vpy1=0.0, vpz1=0.0
real :: xp2=0.0, yp2=0.0, zp2=0.0, vpx2=0.0, vpy2=0.0, vpz2=0.0
real :: xp3=0.0, yp3=0.0, zp3=0.0, vpx3=0.0, vpy3=0.0, vpz3=0.0
real :: Lx0=0.0, Ly0=0.0, Lz0=0.0
real :: delta_vp0=1.0, tausp=0.0, tausp1=0.0, rpbeta=0.0
real :: nu_epicycle=0.0, nu_epicycle2=0.0
real :: beta_dPdr_dust=0.0, beta_dPdr_dust_scaled=0.0
real :: tausg_min=0.0, tausg1_max=0.0, epsp_friction_increase=100.0
real :: cdtp=0.2, cdtpgrav=0.1, cdtp_drag=0.2
real :: gravx=0.0, gravz=0.0, gravr=1.0, kx_gg=1.0, kz_gg=1.0
real :: gravsmooth=0.0, gravsmooth2=0.0, Ri0=0.25, eps1=0.5
real :: kx_xxp=0.0, ky_xxp=0.0, kz_xxp=0.0, amplxxp=0.0
real :: kx_vvp=0.0, ky_vvp=0.0, kz_vvp=0.0, amplvvp=0.0
real :: kx_vpx=0.0, kx_vpy=0.0, kx_vpz=0.0
real :: ky_vpx=0.0, ky_vpy=0.0, ky_vpz=0.0
real :: kz_vpx=0.0, kz_vpy=0.0, kz_vpz=0.0
real :: phase_vpx=0.0, phase_vpy=0.0, phase_vpz=0.0
real :: tstart_dragforce_par=0.0
real :: tstart_grav_par=0.0, tstart_grav_x_par=0.0
real :: tstart_grav_z_par=0.0, tstart_grav_r_par=0.0
real :: tstart_liftforce_par=0.0
real :: tstart_brownian_par=0.0
real :: tstart_collisional_cooling=0.0
real :: tstart_sink_par=0.0
real :: tau_coll_min=0.0, tau_coll1_max=0.0
real :: coeff_restitution=0.5, mean_free_path_gas=0.0
real :: rad_sphere=0.0
real :: a_ellipsoid=0.0, b_ellipsoid=0.0, c_ellipsoid=0.0
real :: a_ell2=0.0, b_ell2=0.0, c_ell2=0.0
real :: taucool=0.0, taucool1=0.0, brownian_T0=0.0, thermophoretic_T0=0.0
real :: xsinkpoint=0.0, ysinkpoint=0.0, zsinkpoint=0.0, rsinkpoint=0.0
real :: particles_insert_rate=0.
real :: avg_n_insert, remaining_particles=0.0
real :: max_particle_insert_time=huge1
real :: Deltauy_gas_friction=0.0
real :: cond_ratio=0.0
real :: pscalar_sink_rate=0.0
real :: frac_init_particles=1.0, particles_insert_ramp_time=0.0
real :: tstart_insert_particles=0.0
real :: birthring_r=1.0, birthring_width=0.1
real :: tstart_rpbeta=0.0, birthring_lifetime=huge1
real :: rdiffconst_dragf=0.07,rdiffconst_pass=0.07
integer :: l_hole=0, m_hole=0, n_hole=0
integer :: iffg=0, ifgx=0, ifgy=0, ifgz=0, ibrtime=0
integer :: istep_dragf=3,istep_pass=3
logical :: ldragforce_dust_par=.false., ldragforce_gas_par=.false.
logical :: ldragforce_stiff=.false., ldragforce_radialonly=.false.
logical :: ldragforce_heat=.false., lcollisional_heat=.false.
logical :: lpar_spec=.false., lcompensate_friction_increase=.false.
logical :: lcollisional_cooling_taucool=.false.
logical :: lcollisional_cooling_rms=.false.
logical :: lcollisional_cooling_twobody=.false.
logical :: lcollisional_dragforce_cooling=.false.
logical :: ltau_coll_min_courant=.true.
logical :: ldragforce_equi_global_eps=.false., ldragforce_equi_noback=.false.
logical :: ldraglaw_epstein=.true., ldraglaw_epstein_stokes_linear=.false.
logical :: ldraglaw_simple=.false.
logical :: ldraglaw_steadystate=.false., ldraglaw_variable=.false.
logical :: ldraglaw_purestokes=.false.
logical :: ldraglaw_stokesschiller=.false.
logical :: ldraglaw_epstein_transonic=.false.
logical :: ldraglaw_eps_stk_transonic=.false.
logical :: ldraglaw_variable_density=.false.
logical :: lcoldstart_amplitude_correction=.false.
logical :: luse_tau_ap=.true.
logical :: lbrownian_forces=.false.
logical :: lthermophoretic_forces=.false.
logical :: lenforce_policy=.false., lnostore_uu=.true.
logical :: ldt_grav_par=.true., ldt_adv_par=.true.
logical :: lsinkpoint=.false., lglobalrandom=.false.
logical :: lcoriolis_force_par=.true., lcentrifugal_force_par=.false.
logical :: lshear_accel_par = .true.
logical :: lcalc_uup=.false.
logical :: lparticle_gravity=.true.
logical :: lcylindrical_gravity_par=.false.
logical :: lpscalar_sink=.false.
logical :: lsherwood_const=.false.
logical :: lbubble=.false.
logical :: linsert_as_many_as_possible=.false.
logical :: lwithhold_init_particles=.false.
logical :: lgaussian_birthring=.false.
logical :: lvector_gravity=.false.
logical :: lpeh_radius=.false.
logical :: lbirthring_depletion=.false.
logical :: lstocunn1 = .false.
logical :: lprecalc_cell_volumes=.false.
logical :: ldiffuse_passive = .false.,ldiff_pass=.false.
logical :: ldiffuse_dragf= .false.,ldiff_dragf=.false.
logical :: lsimple_volume=.false.
!
character (len=labellen) :: interp_pol_uu ='none'
character (len=labellen) :: interp_pol_oo ='none'
character (len=labellen) :: interp_pol_TT ='none'
character (len=labellen) :: interp_pol_rho='none'
character (len=labellen) :: interp_pol_pp ='none'
character (len=labellen) :: interp_pol_species='none'
character (len=labellen) :: interp_pol_gradTT='none'
character (len=labellen) :: interp_pol_nu='none'
!
character (len=labellen), dimension (ninit) :: initxxp='nothing'
character (len=labellen), dimension (ninit) :: initvvp='nothing'
character (len=labellen) :: gravx_profile='', gravz_profile=''
character (len=labellen) :: gravr_profile='', interp_default=''
character (len=labellen) :: thermophoretic_eq= 'nothing'
!
integer :: init_repeat=0 !repeat particle initialization for distance statistics
!
! Interactions with special/shell
!
integer :: k_shell=-1 !k associated with minshell (special/shell.f90)
logical :: l_shell=.false. !using special/shell.f90 for gas velocities
!
real, dimension(3) :: uup_shared=0
real :: turnover_shared=0, nu_draglaw=0.
logical :: vel_call=.false., turnover_call=.false.
logical :: lreassign_strat_rhom=.true., lnu_draglaw=.false.
!
logical :: lcdtp_shear = .true.
!
logical :: lcompensate_sedimentation=.false.
real :: compensate_sedimentation=1.
!
real :: A3=0., A2=0.
!
namelist /particles_init_pars/ &
initxxp, initvvp, xp0, yp0, zp0, vpx0, vpy0, vpz0, delta_vp0, &
ldragforce_gas_par, ldragforce_dust_par, bcpx, bcpy, bcpz, tausp, &
beta_dPdr_dust, mpmat, &
nu_epicycle, rp_int, rp_ext, gravx_profile, gravz_profile, &
gravr_profile, gravx, gravz, gravr, gravsmooth, kx_gg, kz_gg, Ri0, &
eps1, lmigration_redo, ldragforce_equi_global_eps, coeff, kx_vvp, &
ky_vvp, kz_vvp, amplvvp, kx_xxp, ky_xxp, kz_xxp, amplxxp, kx_vpx, &
kx_vpy, kx_vpz, ky_vpx, ky_vpy, ky_vpz, kz_vpx, kz_vpy, kz_vpz, &
phase_vpx, phase_vpy, phase_vpz, lcoldstart_amplitude_correction, &
particle_mesh, linterpolate_spline, &
tstart_dragforce_par, tstart_grav_par, lparticle_gravity,&
tstart_grav_x_par, tstart_grav_z_par,tstart_grav_r_par, taucool, &
lcollisional_cooling_taucool, lcollisional_cooling_rms, &
lcollisional_cooling_twobody, tausp_species, tau_coll_min, &
ltau_coll_min_courant, coeff_restitution, tstart_collisional_cooling, &
tausg_min, l_hole, m_hole, n_hole, &
epsp_friction_increase,lcollisional_dragforce_cooling, ldragforce_heat, &
lcollisional_heat, lcompensate_friction_increase, &
lmigration_real_check, ldraglaw_epstein,ldraglaw_simple, &
ldraglaw_epstein_stokes_linear, &
mean_free_path_gas, ldraglaw_epstein_transonic, lcheck_exact_frontier, &
ldraglaw_eps_stk_transonic, dustdensity_powerlaw, rad_sphere, pos_sphere, &
ldragforce_stiff, &
a_ellipsoid, b_ellipsoid, c_ellipsoid, pos_ellipsoid, &
ldraglaw_steadystate, tstart_liftforce_par, &
ldraglaw_purestokes,rpbeta_species, rpbeta, gab_width, &
tstart_brownian_par, tstart_sink_par, &
lbrownian_forces,lthermophoretic_forces,lenforce_policy, &
interp_pol_uu,interp_pol_oo,interp_pol_TT,interp_pol_rho, &
interp_pol_pp,interp_pol_species,brownian_T0, &
thermophoretic_T0, lnostore_uu, ldt_grav_par, ldragforce_radialonly, &
lsinkpoint, xsinkpoint, ysinkpoint, zsinkpoint, rsinkpoint, &
lcoriolis_force_par, lcentrifugal_force_par, ldt_adv_par, Lx0, Ly0, &
Lz0, lglobalrandom, linsert_particles_continuously, &
lrandom_particle_pencils, lnocalc_np, lnocalc_rhop, &
np_const, rhop_const, particle_radius, &
ldragforce_equi_noback, rhopmat, Deltauy_gas_friction, xp1, &
yp1, zp1, vpx1, vpy1, vpz1, xp2, yp2, zp2, vpx2, vpy2, vpz2, &
xp3, yp3, zp3, vpx3, vpy3, vpz3, lsinkparticle_1, rsinkparticle_1, &
lcalc_uup, temp_grad0, thermophoretic_eq, cond_ratio, interp_pol_gradTT, &
lreassign_strat_rhom, lparticlemesh_pqs_assignment, &
lwithhold_init_particles, frac_init_particles, lvector_gravity, &
birthring_r, birthring_width, lgaussian_birthring, ldraglaw_stokesschiller, &
lbirthring_depletion, &
remove_particle_at_time, remove_particle_criteria, remove_particle_criteria_size
!
namelist /particles_run_pars/ &
bcpx, bcpy, bcpz, tausp, dsnap_par_minor, beta_dPdr_dust, &
ldragforce_gas_par, ldragforce_dust_par, &
cdtp, cdtpgrav, lpar_spec, linterp_reality_check, &
nu_epicycle, gravx_profile, gravz_profile, gravr_profile, gravx, gravz, &
gravr, gravsmooth, kx_gg, kz_gg, lmigration_redo, tstart_dragforce_par, &
tstart_grav_par, tstart_grav_x_par, cdtp_drag, &
tstart_grav_z_par, tstart_grav_r_par, lparticle_gravity, &
particle_mesh, taucool, &
lcollisional_cooling_taucool, lcollisional_cooling_rms, &
lcollisional_cooling_twobody, lcollisional_dragforce_cooling, &
tau_coll_min, ltau_coll_min_courant, coeff_restitution, &
tstart_collisional_cooling, tausg_min, epsp_friction_increase, &
ldragforce_heat, lcollisional_heat, lcompensate_friction_increase, &
lmigration_real_check,ldraglaw_variable, luse_tau_ap, &
ldraglaw_epstein, ldraglaw_simple, rdiffconst_pass,&
ldraglaw_epstein_stokes_linear, mean_free_path_gas, &
ldraglaw_epstein_transonic, lcheck_exact_frontier, &
ldraglaw_eps_stk_transonic, ldragforce_stiff, &
ldraglaw_variable_density, ldraglaw_steadystate, tstart_liftforce_par, &
ldraglaw_purestokes, ldiffuse_dragf, ldiffuse_passive, rdiffconst_dragf,&
tstart_brownian_par, tstart_sink_par, ldiff_dragf, ldiff_pass,&
lbrownian_forces, lenforce_policy, &
interp_pol_uu,interp_pol_oo, interp_pol_TT, interp_pol_rho, &
interp_pol_pp,interp_pol_species, &
brownian_T0,thermophoretic_T0, lnostore_uu, ldt_grav_par, &
ldragforce_radialonly, lsinkpoint, xsinkpoint, ysinkpoint, zsinkpoint, &
rsinkpoint, lshear_accel_par, lcoriolis_force_par, &
lcentrifugal_force_par, ldt_adv_par, &
linsert_particles_continuously, particles_insert_rate, &
max_particle_insert_time, lrandom_particle_pencils, lnocalc_np, &
lnocalc_rhop, lstocunn1, istep_dragf, istep_pass,&
np_const, rhop_const, particle_radius, lprecalc_cell_volumes, &
Deltauy_gas_friction, loutput_psize_dist, log_ap_min_dist, &
log_ap_max_dist, nbin_ap_dist, lsinkparticle_1, rsinkparticle_1, &
lthermophoretic_forces, temp_grad0, &
thermophoretic_eq, cond_ratio, interp_pol_gradTT, lcommunicate_rhop, &
lcommunicate_np, lcylindrical_gravity_par, &
l_shell, k_shell, lparticlemesh_pqs_assignment, pscalar_sink_rate, &
lpscalar_sink, lsherwood_const, lnu_draglaw, nu_draglaw,lbubble, &
rpbeta_species, rpbeta, gab_width, initxxp, initvvp, &
particles_insert_ramp_time, tstart_insert_particles, birthring_r, &
birthring_width, lsimple_volume,&
lgaussian_birthring, tstart_rpbeta, linsert_as_many_as_possible, &
lvector_gravity, lcompensate_sedimentation,compensate_sedimentation, &
lpeh_radius, A3, A2, ldraglaw_stokesschiller, lbirthring_depletion, birthring_lifetime, &
remove_particle_at_time, remove_particle_criteria, remove_particle_criteria_size, &
rp_int, rp_ext
!
integer :: idiag_xpm=0, idiag_ypm=0, idiag_zpm=0 ! DIAG_DOC: $x_{part}$
integer :: idiag_xp2m=0, idiag_yp2m=0, idiag_zp2m=0 ! DIAG_DOC: $x^2_{part}$
integer :: idiag_vrelpabsm=0 ! DIAG_DOC: $\rm{Absolute value of mean relative velocity}$
integer :: idiag_rpm=0, idiag_rp2m=0
integer :: idiag_vpxm=0, idiag_vpym=0, idiag_vpzm=0 ! DIAG_DOC: $u_{part}$
integer :: idiag_vpx2m=0, idiag_vpy2m=0, idiag_vpz2m=0 ! DIAG_DOC: $u^2_{part}$
integer :: idiag_ekinp=0 ! DIAG_DOC: $E_{kin,part}$
integer :: idiag_vpxmax=0, idiag_vpymax=0, idiag_vpzmax=0, idiag_vpmax=0 ! DIAG_DOC: $MAX(u_{part})$
integer :: idiag_vpxmin=0, idiag_vpymin=0, idiag_vpzmin=0 ! DIAG_DOC: $MIN(u_{part})$
integer :: idiag_urel=0
integer :: idiag_vpxvpym=0, idiag_vpxvpzm=0, idiag_vpyvpzm=0
integer :: idiag_rhopvpxm=0, idiag_rhopvpym=0, idiag_rhopvpzm=0
integer :: idiag_rhopvpxt=0, idiag_rhopvpyt=0, idiag_rhopvpzt=0
integer :: idiag_rhopvpysm=0
integer :: idiag_lpxm=0, idiag_lpym=0, idiag_lpzm=0
integer :: idiag_lpx2m=0, idiag_lpy2m=0, idiag_lpz2m=0
integer :: idiag_npm=0, idiag_np2m=0, idiag_npmax=0, idiag_npmin=0 ! DIAG_DOC: $\rm{mean particle number density}$
integer :: idiag_dtdragp=0
integer :: idiag_nparmin=0, idiag_nparmax=0, idiag_npargone=0
integer :: idiag_nparsum=0
integer :: idiag_rhopm=0, idiag_rhoprms=0, idiag_rhop2m=0, idiag_rhopmax=0
integer :: idiag_rhopmin=0, idiag_decollp=0, idiag_rhopmphi=0
integer :: idiag_epspmin=0, idiag_epspmax=0, idiag_epspm=0
integer :: idiag_npmx=0, idiag_npmy=0, idiag_npmz=0
integer :: idiag_rhopmx=0, idiag_rhopmy=0, idiag_rhopmz=0
integer :: idiag_rhop2mx=0, idiag_rhop2my=0, idiag_rhop2mz=0
integer :: idiag_epspmx=0, idiag_epspmy=0, idiag_epspmz=0
integer :: idiag_mpt=0, idiag_dedragp=0, idiag_rhopmxy=0, idiag_rhopmr=0
integer :: idiag_sigmap=0
integer :: idiag_dvpx2m=0, idiag_dvpy2m=0, idiag_dvpz2m=0
integer :: idiag_dvpm=0, idiag_dvpmax=0, idiag_epotpm=0
integer :: idiag_rhopmxz=0, idiag_nparpmax=0, idiag_npmxy=0
integer :: idiag_eccpxm=0, idiag_eccpym=0, idiag_eccpzm=0
integer :: idiag_eccpx2m=0, idiag_eccpy2m=0, idiag_eccpz2m=0
integer :: idiag_vprms=0, idiag_vpyfull2m=0, idiag_deshearbcsm=0
integer :: idiag_Shm=0
integer, dimension(ndustrad) :: idiag_npvzmz=0, idiag_npvz2mz=0, idiag_nptz=0
integer, dimension(ndustrad) :: idiag_npuzmz=0
!
contains
!***********************************************************************
subroutine register_particles
!
! Set up indices for access to the fp and dfp arrays
!
! 29-dec-04/anders: coded
!
use FArrayManager, only: farray_register_auxiliary
!
if (lroot) call svn_id( &
"$Id$")
!
! Indices for particle position.
!
ixp=npvar+1
pvarname(npvar+1)='ixp'
iyp=npvar+2
pvarname(npvar+2)='iyp'
izp=npvar+3
pvarname(npvar+3)='izp'
!
! Indices for particle velocity.
!
ivpx=npvar+4
pvarname(npvar+4)='ivpx'
ivpy=npvar+5
pvarname(npvar+5)='ivpy'
ivpz=npvar+6
pvarname(npvar+6)='ivpz'
!
! Increase npvar accordingly.
!
npvar=npvar+6
!
! Set indices for particle assignment.
!
if (.not. lnocalc_np) call farray_register_auxiliary('np',inp,&
communicated=lcommunicate_np)
if (lpeh_radius) then
lnocalc_rhop=.false.
call farray_register_auxiliary('peh',ipeh,communicated=.false.)
endif
if (.not. lnocalc_rhop) call farray_register_auxiliary('rhop',irhop, &
communicated=lparticles_sink.or.lcommunicate_rhop)
if (lcalc_uup .or. ldragforce_stiff) then
call farray_register_auxiliary('uup',iuup,communicated=.true.,vector=3)
iupx=iuup; iupy=iuup+1; iupz=iuup+2
endif
!
! Special variable for stiff drag force equations.
!
if (ldragforce_stiff) then
call farray_register_auxiliary('ffg',iffg,communicated=.true.,vector=3)
ifgx=iffg; ifgy=iffg+1; ifgz=iffg+2
endif
!
! Special variable for diffusion of the passive scalar consumption
!
if (ldiffuse_passive .and. lpscalar_sink) then
call farray_register_auxiliary('dlncc',idlncc,communicated=.true.)
elseif (ldiffuse_passive .and. .not. lpscalar_sink) then
call fatal_error('particles_lagrangian','ldiffuse_passive needs lpscalar_sink')
endif
!
! Special variable for diffusion of particle-gas dragforce
!
if (ldiffuse_dragf .and. ldragforce_gas_par) then
call farray_register_auxiliary('dfg',idfg,communicated=.true.,vector=3)
idfx=idfg; idfy=idfg+1; idfz=idfg+2
elseif (ldiffuse_dragf .and. .not. ldragforce_gas_par) then
call fatal_error('particles_lagrangian','ldiffuse_dragf needs ldragforce_gas_par')
endif
!
! Kill particles that spend enough time in birth ring
if (lbirthring_depletion) then
ibrtime = mpvar+npaux+1
pvarname(ibrtime) = 'ibrtime'
npaux = npaux+1
endif
!
! Check that the fp and dfp arrays are big enough.
!
if (npvar > mpvar) then
if (lroot) write(0,*) 'npvar = ', npvar, ', mpvar = ', mpvar
call fatal_error('register_particles','npvar > mpvar')
endif
!
endsubroutine register_particles
!***********************************************************************
subroutine initialize_particles(f)
!
! Perform any post-parameter-read initialization i.e. calculate derived
! parameters.
!
! 29-dec-04/anders: coded
! 5-mar-15/MR: reference state included in calculation of mean density
!
use EquationOfState, only: rho0, cs0
use SharedVariables, only: put_shared_variable, get_shared_variable
use Density, only: mean_density
!
real, dimension (mx,my,mz,mfarray) :: f
!
real :: rhom
integer :: ierr, jspec
logical, pointer :: lshearadvection_as_shift
real, pointer :: reference_state_mass
!
! This module is incompatible with particle block domain decomposition.
!
if (lparticles_blocks) then
if (lroot) then
print*, 'initialize_particles: must use PARTICLES = PARTICLES_DUST_BLOCKS'
print*, ' with particle block domain decomposition'
endif
call fatal_error('initialize_particles','')
endif
!
! Hand over Coriolis force and shear acceleration to Particles_drag.
!
drag: if (lparticles_drag) then
coriolis: if (lcoriolis_force_par) then
lcoriolis_force_par = .false.
if (lroot) print *, 'initialize_particles: turned off and hand over Coriolis force to Particles_drag. '
endif coriolis
shacc: if (lshear .and. lshear_accel_par) then
lshear_accel_par = .false.
if (lroot) print *, 'initialize_particles: turned off and hand over shear acceleration to Particles_drag. '
endif shacc
endif drag
!
! Check if shear advection is on and decide if it needs to be included in the timestep condition.
!
shear: if (lshear) then
call get_shared_variable('lshearadvection_as_shift', lshearadvection_as_shift)
lcdtp_shear = .not. lshearadvection_as_shift
nullify(lshearadvection_as_shift)
endif shear
!
! Report the particle radius, if set.
!
apar: if (particle_radius /= 0.0) then
if (lparticles_radius) &
call fatal_error('initialize_particles', 'particle_radius /= 0 has no effect when module Particles_radius is on. ')
if (lroot) print *, 'initialize_particles: radius of each constituent particle = ', particle_radius
endif apar
!
! The inverse stopping time is needed for drag force and collisional cooling.
!
if (tausp/=0.0) tausp1=1/tausp
!
! Inverse cooling time.
!
if (taucool/=0.0) taucool1=1/taucool
!
! Inverse material density.
!
if (rhopmat/=0.0) rhopmat1=1/rhopmat
!
! Multiple dust species. Friction time is given in the array tausp_species.
!
if (npar_species>1) then
if (lroot) then
print*, 'initialize_particles: '// &
'Number of particle species = ', npar_species
print*, 'initialize_particles: tausp_species = ', tausp_species
endif
!
! Must have set tausp_species for drag force.
!
if (ldragforce_dust_par .or. ldragforce_gas_par) then
if (any(tausp_species==0.0)) then
if (lroot) print*, &
'initialize_particles: drag force must have tausp_species/=0.0!'
call fatal_error('initialize_particles','')
endif
!
! Inverse friction time is needed for drag force.
!
do jspec=1,npar_species
if (tausp_species(jspec)/=0.0) &
tausp1_species(jspec)=1/tausp_species(jspec)
enddo
endif
else
!
! Single dust species => If tausp_species is set, it is probably an error.
!
if (any(tausp_species/=0.0) .and. tausp/=tausp_species(1)) then
call fatal_error('initialize_particles', 'When there is only '// &
'1 particle species, use tausp instead of tausp_species')
endif
tausp_species(1)=tausp
if (tausp_species(1)/=0.0) tausp1_species(1)=1/tausp_species(1)
endif
!
! Share friction time (but only if Epstein drag regime!).
!
if (ldraglaw_epstein .or. ldraglaw_simple) then
call put_shared_variable( 'tausp_species', tausp_species)
call put_shared_variable('tausp1_species',tausp1_species)
endif
!
! Global gas pressure gradient seen from the perspective of the dust.
!
if (beta_dPdr_dust/=0.0) then
beta_dPdr_dust_scaled=beta_dPdr_dust*Omega/cs0
if (lroot) print*, 'initialize_particles: Global pressure '// &
'gradient with beta_dPdr_dust=', beta_dPdr_dust
endif
!
! Calculate nu_epicycle**2 for gravity.
!
if (gravz_profile=='' .and. nu_epicycle/=0.0) gravz_profile='linear'
nu_epicycle2=nu_epicycle**2
!
! Calculate gravsmooth**2 for gravity.
!
if (gravsmooth/=0.0) gravsmooth2=gravsmooth**2
!
! Share Keplerian gravity.
!
call put_shared_variable('gravr',gravr,ierr)
if (ierr/=0) call fatal_error('initialize_particles', &
'there was a problem when sharing gravr')
!
! Inverse of minimum gas friction time (time-step control).
!
if (tausg_min/=0.0) tausg1_max=1.0/tausg_min
!
! Set minimum collisional time-scale so that time-step is not affected.
!
if (lrun .and. ltau_coll_min_courant) then
if (cs0==impossible) then
tau_coll_min=impossible
else
tau_coll_min=2*dx/cs0
endif
if (lroot) print*, 'initialize particles: set minimum collisional '// &
'time-scale equal to two times the Courant time-step.'
endif
!
! Inverse of minimum collisional time-scale.
!
if (lrun .and. tau_coll_min>0.0) tau_coll1_max=1/tau_coll_min
!
! Gas density is needed for back-reaction friction force.
!
if (ldragforce_gas_par .and. .not. ldensity) then
if (lroot) then
print*, 'initialize_particles: friction force on gas only works '
print*, ' together with gas density module!'
endif
call fatal_error('initialize_particles','')
endif
!
! Need to map particles on the grid for dragforce on gas.
!
if (ldragforce_gas_par) then
!
! When drag force is smoothed, df is also set in the first ghost zone. This
! region needs to be folded back into the df array after pde is finished,
!
if (particle_mesh .eq. 'cic' .or. particle_mesh .eq. 'tsc') &
lfold_df=.true.
endif
!
if (lcollisional_cooling_twobody) then
allocate(kneighbour(mpar_loc))
lshepherd_neighbour=.true.
endif
!
if (ldraglaw_epstein_stokes_linear) ldraglaw_epstein=.false.
if (ldraglaw_epstein_transonic .or.&
ldraglaw_eps_stk_transonic .or.&
ldraglaw_steadystate .or.&
ldraglaw_purestokes .or.&
ldraglaw_stokesschiller .or.&
ldraglaw_simple) then
ldraglaw_epstein=.false.
endif
if (ldraglaw_epstein_transonic .and.&
ldraglaw_eps_stk_transonic) then
print*,'both epstein and epstein-stokes transonic '//&
'drag laws are switched on. You cannot have '//&
'both. Stop and choose only one.'
call fatal_error('initialize_particles','')
endif
!
! Stiff drag force approximation.
!
if (ldragforce_stiff) then
if (ldragforce_dust_par .or. ldragforce_gas_par) then
if (lroot) print*, 'initialize_particles: stiff drag force '// &
'approximation is incompatible with normal drag'
call fatal_error('initialize_particles','')
endif
f(l1:l2,m1:m2,n1:n2,ifgx:ifgz)=0.0
endif
!
! Initialize storage of energy gain released by shearing boundaries.
!
if (idiag_deshearbcsm/=0) energy_gain_shear_bcs=0.0
!
! Fatal error if sink particle radius is zero or negative.
!
if (lsinkparticle_1 .and. rsinkparticle_1<=0.0) then
if (lroot) print*, 'initialize_particles: sink particle radius is '// &
'zero or negative: ', rsinkparticle_1
call fatal_error('initialize_particles','')
endif
!
! Particle self gravity for x,z,r direction
!
if ( tstart_grav_par > 0.0) then
tstart_grav_x_par = tstart_grav_par
tstart_grav_z_par = tstart_grav_par
tstart_grav_r_par = tstart_grav_par
endif
!
! Die if lcompensate_sedimentation is used and the vertical direction is present.
!
if (lcompensate_sedimentation.and.&
( (nzgrid/=1.and.(lcartesian_coords.or.lcylindrical_coords)).or.&
(nygrid/=1.and.lspherical_coords)&
)&
) call fatal_error("initialize_particles",&
"compensate_sedimentation should only be used when the vertical dimension is not present")
!
! Set up interpolation logicals. These logicals can be OR'ed with some logical
! in the other particle modules' initialization subroutines to enable
! interpolation based on some condition local to that module.
! (The particles_spin module will for instance enable interpolation of the
! vorticity oo)
!
uu: if (lnostore_uu) then
if (ldraglaw_steadystate .or. lparticles_spin .or. lsolid_ogrid) &
call fatal_error('initialize_particles', 'lnostore_uu = .false. is required. ')
interp%luu = .false.
else uu
interp%luu = ldragforce_dust_par .or. ldraglaw_steadystate .or. lparticles_spin .or. lsolid_ogrid
endif uu
interp%loo=.false.
interp%lTT=(lbrownian_forces.and.(brownian_T0==0.0))&
.or.(lthermophoretic_forces.and.(thermophoretic_T0==0.0))&
.or.lparticles_temperature
interp%lgradTT=lthermophoretic_forces .and. (temp_grad0(1)==0.0) &
.and. (temp_grad0(2)==0.0) .and. (temp_grad0(3)==0.0)
interp%lrho=lbrownian_forces.or.ldraglaw_steadystate &
.or. lthermophoretic_forces .or. ldraglaw_purestokes &
.or. ldraglaw_stokesschiller
interp%lnu=lchemistry
interp%lpp=lparticles_chemistry
interp%lspecies=lparticles_surfspec
!
! The default interpolation policy is given by the user through the
! input variable "particle_mesh". If this input is not given, the
! default policy is 'ngp' (nearest grid point).
!
if (particle_mesh=='') then
interp_default='ngp'
else
interp_default=particle_mesh
endif
print*,'interp_default=',interp_default,particle_mesh
!
! If individual interpolation policies are given, these will
! overwrite the default one.
!
if (interp_pol_uu .eq. 'none') interp_pol_uu=interp_default
if (interp_pol_oo .eq. 'none') interp_pol_oo=interp_default
if (interp_pol_TT .eq. 'none') interp_pol_TT=interp_default
if (interp_pol_rho .eq. 'none') interp_pol_rho=interp_default
if (interp_pol_pp .eq. 'none') interp_pol_pp=interp_default
if (interp_pol_nu .eq. 'none') interp_pol_nu=interp_default
if (interp_pol_species.eq. 'none') &
interp_pol_species=interp_default
if (interp_pol_gradTT .eq. 'none') &
interp_pol_gradTT=interp_default
!
! Overwrite with new policy variables:
!
select case (interp_pol_uu)
case ('tsc')
interp%pol_uu=tsc
case ('cic')
interp%pol_uu=cic
case ('ngp')
interp%pol_uu=ngp
case default
call fatal_error('initialize_particles','No such value for '// &
'interp_pol_uu: '//trim(interp_pol_uu))
endselect
!
select case (interp_pol_oo)
case ('tsc')
interp%pol_oo=tsc
case ('cic')
interp%pol_oo=cic
case ('ngp')
interp%pol_oo=ngp
case default
call fatal_error('initialize_particles','No such such value for '// &
'interp_pol_oo: '//trim(interp_pol_oo))
endselect
!
select case (interp_pol_TT)
case ('tsc')
interp%pol_TT=tsc
case ('cic')
interp%pol_TT=cic
case ('ngp')
interp%pol_TT=ngp
case default
call fatal_error('initialize_particles','No such such value for '// &
'interp_pol_TT: '//trim(interp_pol_TT))
endselect
!
select case (interp_pol_gradTT)
case ('tsc')
call fatal_error('initialize_particles','Not implemented gradTT'// &
'interp_pol_gradTT: '//trim(interp_pol_gradTT))
case ('cic')
call fatal_error('initialize_particles','Not implemented gradTT'// &
'interp_pol_gradTT: '//trim(interp_pol_gradTT))
case ('ngp')
interp%pol_gradTT=ngp
case default
call fatal_error('initialize_particles','No such such value for '// &
'interp_pol_gradTT: '//trim(interp_pol_gradTT))
endselect
!
select case (interp_pol_rho)
case ('tsc')
interp%pol_rho=tsc
case ('cic')
interp%pol_rho=cic
case ('ngp')
interp%pol_rho=ngp
case default
call fatal_error('initialize_particles','No such such value for '// &
'interp_pol_rho: '//trim(interp_pol_rho))
endselect
!
select case (interp_pol_pp)
case ('tsc')
call fatal_error('initialize_particles','Not implemented pp'// &
'interp_pol_pp: '//trim(interp_pol_pp))
case ('cic')
call fatal_error('initialize_particles','Not implemented pp'// &
'interp_pol_pp: '//trim(interp_pol_pp))
case ('ngp')
interp%pol_pp=ngp
case default
call fatal_error('initialize_particles','No such such value for '// &
'interp_pol_pp: '//trim(interp_pol_pp))
endselect
!
select case (interp_pol_nu)
case ('tsc')
call fatal_error('initialize_particles','Not implemented nu'// &
'interp_pol_nu: '//trim(interp_pol_nu))
case ('cic')
call fatal_error('initialize_particles','Not implemented nu'// &
'interp_pol_nu: '//trim(interp_pol_nu))
case ('ngp')
interp%pol_nu=ngp
case default
call fatal_error('initialize_particles','No such such value for '// &
'interp_pol_nu: '//trim(interp_pol_nu))
endselect
!
if (l_shell) then
if ( k_shell < 0) call fatal_error('initialize_particles','Set k_shell')
call put_shared_variable('uup_shared',uup_shared,ierr)
call put_shared_variable('vel_call',vel_call,ierr)
call put_shared_variable('turnover_call',turnover_call,ierr)
call put_shared_variable('turnover_shared',turnover_shared,ierr)
endif
!
call keep_compiler_quiet(f)
!
endsubroutine initialize_particles
!***********************************************************************
subroutine init_particles(f,fp,ineargrid)
!
! Initial positions and velocities of dust particles.
!
! 29-dec-04/anders: coded
!
use Density, only: beta_glnrho_global
use EquationOfState, only: cs20
use General, only: random_number_wrapper, normal_deviate
use Mpicomm, only: mpireduce_sum, mpibcast_real
use InitialCondition, only: initial_condition_xxp, initial_condition_vvp
use Particles_diagnos_dv, only: repeated_init
!
real, dimension (mx,my,mz,mfarray), intent (out) :: f
real, dimension (mpar_loc,mparray), intent (out) :: fp
integer, dimension (mpar_loc,3), intent (out) :: ineargrid
real, dimension (mpar_loc) :: rr_tmp, az_tmp
!
real, dimension (3) :: uup, Lxyz_par, xyz0_par, xyz1_par
real :: vpx_sum, vpy_sum, vpz_sum
real :: r, p, q, px, py, pz, eps, cs, k2_xxp, rp2
real :: dim1, npar_loc_x, npar_loc_y, npar_loc_z, dx_par, dy_par, dz_par
real :: rad,rad_scl,phi,tht,tmp,OO,xx0,yy0,r2
integer :: l, j, k, ix0, iy0, iz0
logical :: lequidistant=.false.
real :: rpar_int,rpar_ext
!
! Optionally withhold some number of particles, to be inserted in
! insert_particles. The particle indices to be removed are not randomized,
! so any randomization needs to be taken care of in the above initxxp cases.
!
if (lwithhold_init_particles .and. frac_init_particles < 1.0-tini) then
npar_loc = nint(frac_init_particles*real(npar_loc))
endif
!
! 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.
!
do j=1,ninit
!
select case (initxxp(j))
!
case ('nothing')
if (lroot .and. j==1) print*, 'init_particles: nothing'
!
case ('origin')
if (lroot) print*, 'init_particles: All particles at origin'
fp(1:npar_loc,ixp:izp)=0.0
!
case ('zero-z')
if (lroot) print*, 'init_particles: Zero z coordinate'
fp(1:npar_loc,izp)=0.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 ('constant-1')
if (lroot) &
print*, 'init_particles: Particle 1 at x,y,z=', xp1, yp1, zp1
do k=1,npar_loc
if (ipar(k)==1) then
fp(k,ixp)=xp1
fp(k,iyp)=yp1
fp(k,izp)=zp1
endif
enddo
!
case ('constant-2')
if (lroot) &
print*, 'init_particles: Particle 2 at x,y,z=', xp2, yp2, zp2
do k=1,npar_loc
if (ipar(k)==2) then
fp(k,ixp)=xp2
fp(k,iyp)=yp2
fp(k,izp)=zp2
endif
enddo
!
case ('constant-3')
if (lroot) &
print*, 'init_particles: Particle 2 at x,y,z=', xp3, yp3, zp3
do k=1,npar_loc
if (ipar(k)==3) then
fp(k,ixp)=xp3
fp(k,iyp)=yp3
fp(k,izp)=zp3
endif
enddo
!
case ('random-constz')
if (lroot) print*, 'init_particles: Random particle positions'
do k=1,npar_loc
if (nxgrid/=1) then
call random_number_wrapper(r)
fp(k,ixp)=r
endif
if (nygrid/=1) then
call random_number_wrapper(r)
fp(k,iyp)=r
endif
enddo
if (nxgrid/=1) &
fp(1:npar_loc,ixp)=xyz0_par(1)+fp(1:npar_loc,ixp)*Lxyz_par(1)
if (nygrid/=1) &
fp(1:npar_loc,iyp)=xyz0_par(2)+fp(1:npar_loc,iyp)*Lxyz_par(2)
if (nzgrid/=1) &
fp(1:npar_loc,izp)=zp0
!
case ('random')
if (lroot) print*, 'init_particles: Random particle positions'
do k=1,npar_loc
if (nxgrid/=1) then
call random_number_wrapper(r)
fp(k,ixp)=r
endif
if (nygrid/=1) then
call random_number_wrapper(r)
fp(k,iyp)=r
endif
if (nzgrid/=1) then
call random_number_wrapper(r)
fp(k,izp)=r
endif
enddo
if (nxgrid/=1) &
fp(1:npar_loc,ixp)=xyz0_par(1)+fp(1:npar_loc,ixp)*Lxyz_par(1)
if (nygrid/=1) &
fp(1:npar_loc,iyp)=xyz0_par(2)+fp(1:npar_loc,iyp)*Lxyz_par(2)
if (nzgrid/=1) &
fp(1:npar_loc,izp)=xyz0_par(3)+fp(1:npar_loc,izp)*Lxyz_par(3)
case ('random-shift')
if (lroot) print*, 'init_particles: Randomly distribute particle positions, then shift them'
k2_xxp = kx_xxp**2 + ky_xxp**2 + kz_xxp**2
do k=1, npar_loc
if (nxgrid /= 1) then
call random_number_wrapper(r)
fp(k, ixp) = xyz0_par(1) + r * Lxyz_par(1) - &
kx_xxp/k2_xxp * amplxxp * sin(kx_xxp*fp(k, ixp) + ky_xxp*fp(k, iyp) + kz_xxp*fp(k, izp))
endif
if (nygrid /= 1) then
call random_number_wrapper(r)
fp(k, iyp) = xyz0_par(2) + r * Lxyz_par(2) - &
ky_xxp/k2_xxp * amplxxp * sin(kx_xxp*fp(k, ixp) + ky_xxp*fp(k, iyp) + kz_xxp*fp(k, izp))
endif
if (nzgrid /= 1) then
call random_number_wrapper(r)
fp(k, izp) = xyz0_par(3) + r * Lxyz_par(3) - &
kz_xxp/k2_xxp * amplxxp * sin(kx_xxp*fp(k, ixp) + ky_xxp*fp(k, iyp) + kz_xxp*fp(k, izp))
endif
enddo
!
case ('random-circle')
if (lroot) print*, 'init_particles: Random particle positions'
do k=1,npar_loc
call random_number_wrapper(r)
if (zp0>yp0) then
fp(k,ixp)=xp0*cos((zp0-yp0)*r+yp0)
fp(k,iyp)=xp0*sin((zp0-yp0)*r+yp0)
else
fp(k,ixp)=xp0*cos(2*pi*r)
fp(k,iyp)=xp0*sin(2*pi*r)
endif
enddo
!
case ('random-sphere')
if (lroot) print*, 'init_particles: Random particle positions '// &
'in a sphere around (0,0,0) with radius=',rad_sphere
if (rad_sphere==0) then
call fatal_error('init_particles','random-sphere '// &
'radius needs to be larger than zero')
endif
if (-rad_sphere+pos_sphere(1)<xyz0(1) .or. &
+rad_sphere+pos_sphere(1)>xyz1(1) .or. &
-rad_sphere+pos_sphere(2)<xyz0(2) .or. &
+rad_sphere+pos_sphere(2)>xyz1(2) .or. &
-rad_sphere+pos_sphere(3)<xyz0(3) .or. &
+rad_sphere+pos_sphere(3)>xyz1(3)) then
call fatal_error('init_particles','random-sphere '// &
'sphere needs to fit in the box')
endif
if (lcartesian_coords) then
do k=1,npar_loc
rp2=2.*rad_sphere**2
do while (rp2>rad_sphere**2)
call random_number_wrapper(r)
fp(k,ixp)=(r-0.5)*2.*rad_sphere
call random_number_wrapper(r)
fp(k,iyp)=(r-0.5)*2.*rad_sphere
call random_number_wrapper(r)
fp(k,izp)=(r-0.5)*2.*rad_sphere
rp2=fp(k,ixp)**2+fp(k,iyp)**2+fp(k,izp)**2
enddo
fp(k,ixp)=fp(k,ixp)+pos_sphere(1)
fp(k,iyp)=fp(k,iyp)+pos_sphere(2)
fp(k,izp)=fp(k,izp)+pos_sphere(3)
enddo
else
call fatal_error('init_particles','random-sphere '// &
'only implemented for cartesian coordinates')
endif
!
case ('random-ellipsoid')
if (lroot) print*, 'init_particles: Random particle positions '// &
'in an ellipsoid around ', pos_ellipsoid, ' with ' // &
'semi-principal axes a,b,c =',a_ellipsoid,b_ellipsoid,c_ellipsoid
if ((a_ellipsoid==0) .or. (b_ellipsoid==0) .or. (c_ellipsoid==0)) then
call fatal_error('init_particles','random-ellipsoid '// &
'all semi-principal axes need to be larger than zero')
endif
if (-a_ellipsoid+pos_ellipsoid(1)<xyz0(1) .or. &
+a_ellipsoid+pos_ellipsoid(1)>xyz1(1) .or. &
-b_ellipsoid+pos_ellipsoid(2)<xyz0(2) .or. &
+b_ellipsoid+pos_ellipsoid(2)>xyz1(2) .or. &
-c_ellipsoid+pos_ellipsoid(3)<xyz0(3) .or. &
+c_ellipsoid+pos_ellipsoid(3)>xyz1(3)) then
call fatal_error('init_particles','random-ellipsoid '// &
'ellipsoid needs to fit in the box')
endif
if (lcartesian_coords) then
a_ell2=a_ellipsoid**2
b_ell2=b_ellipsoid**2
c_ell2=c_ellipsoid**2
do k=1,npar_loc
rp2=2.
do while (rp2>1.)
call random_number_wrapper(r)
fp(k,ixp)=(r-0.5)*2.*a_ellipsoid
call random_number_wrapper(r)
fp(k,iyp)=(r-0.5)*2.*b_ellipsoid
call random_number_wrapper(r)
fp(k,izp)=(r-0.5)*2.*c_ellipsoid
rp2=fp(k,ixp)**2/a_ell2+fp(k,iyp)**2/b_ell2+fp(k,izp)**2/c_ell2
enddo
fp(k,ixp)=fp(k,ixp)+pos_ellipsoid(1)
fp(k,iyp)=fp(k,iyp)+pos_ellipsoid(2)
fp(k,izp)=fp(k,izp)+pos_ellipsoid(3)
enddo
else
call fatal_error('init_particles','random-ellipsoid '// &
'only implemented for cartesian coordinates')
endif
!
case ('random-line-x')
if (lroot) print*, 'init_particles: Random particle positions'
do k=1,npar_loc
if (nxgrid/=1) then
call random_number_wrapper(r)
fp(k,ixp)=r
endif
enddo
if (nxgrid/=1) &
fp(1:npar_loc,ixp)=xyz0_par(1)+fp(1:npar_loc,ixp)*Lxyz_par(1)
fp(1:npar_loc,iyp)=yp0
fp(1:npar_loc,izp)=zp0
!
case ('random-line-y')
if (lroot) print*, 'init_particles: Random particle positions'
do k=1,npar_loc
if (nygrid/=1) then
call random_number_wrapper(r)
fp(k,iyp)=r
endif
enddo
if (nygrid/=1) &
fp(1:npar_loc,iyp)=xyz0_par(2)+fp(1:npar_loc,iyp)*Lxyz_par(2)
fp(1:npar_loc,ixp)=xp0
fp(1:npar_loc,izp)=zp0
!
case ('random-hole')
if (lroot) print*, 'init_particles: Random particle positions '// &
'with inner hole'
do k=1,npar_loc
rp2=-1.0
do while (rp2<rp_int**2)
rp2=0.0
if (nxgrid/=1) then
call random_number_wrapper(r)
fp(k,ixp)=xyz0(1)+r*Lxyz(1)
rp2=rp2+fp(k,ixp)**2
endif
if (nygrid/=1) then
call random_number_wrapper(r)
fp(k,iyp)=xyz0(2)+r*Lxyz(2)
rp2=rp2+fp(k,iyp)**2
endif
if (nzgrid/=1) then
call random_number_wrapper(r)
fp(k,izp)=xyz0(3)+r*Lxyz(3)
rp2=rp2+fp(k,izp)**2
endif
enddo
enddo
!
case ('random-box')
if (lroot) print*, 'init_particles: Random particle positions '// &
'within a box'
do k=1,npar_loc
if (nxgrid/=1) then
call random_number_wrapper(r)
fp(k,ixp)=r
endif
if (nygrid/=1) then
call random_number_wrapper(r)
fp(k,iyp)=r
endif
if (nzgrid/=1) then
call random_number_wrapper(r)
fp(k,izp)=r
endif
if (lcylindrical_coords) then
xx0=xp0+fp(k,ixp)*Lx0
yy0=yp0+fp(k,iyp)*Ly0
r2=xx0**2+yy0**2
if (nxgrid/=1) fp(k,ixp)=sqrt(r2)
if (nygrid/=1) fp(k,iyp)=atan(yy0/xx0)+pi*(xx0/abs(xx0)-1)*0.5
if (nzgrid/=1) fp(k,izp)=zp0+fp(k,izp)*Lz0
else
if (nxgrid/=1) fp(k,ixp)=xp0+fp(k,ixp)*Lx0
if (nygrid/=1) fp(k,iyp)=yp0+fp(k,iyp)*Ly0
if (nzgrid/=1) fp(k,izp)=zp0+fp(k,izp)*Lz0
endif
enddo
!
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
!
if (lcylindrical_coords.or.lcartesian_coords) then
tmp=2-dustdensity_powerlaw
elseif (lspherical_coords) then
tmp=3-dustdensity_powerlaw
else
call fatal_error("init_particles",&
"The world is flat, and we never got here")
endif
!
if (lcartesian_coords) then
if (nprocx==1) then
rpar_int=rp_int
rpar_ext=rp_ext
else
call fatal_error("init_particles",&
"random-cyl not yet ready for nprocx/=1 in Cartesian. Parallelize in y or z")
endif
else
rpar_int = xyz0_loc(1)
rpar_ext = xyz1_loc(1)
endif
call random_number_wrapper(rad_scl)
rad_scl = rpar_int**tmp + rad_scl*(rpar_ext**tmp-rpar_int**tmp)
rad = rad_scl**(1./tmp)
!
! Random in azimuth
!
if (lcartesian_coords) then
call random_number_wrapper(phi)
phi = 2*pi*phi
if (nxgrid/=1) fp(k,ixp)=rad*cos(phi)
if (nygrid/=1) fp(k,iyp)=rad*sin(phi)
if (nzgrid/=1) then
call random_number_wrapper(r)
fp(k,izp)=xyz0_par(3)+r*Lxyz_par(3)
endif
elseif (lcylindrical_coords) then
if (nxgrid/=1) fp(k,ixp)=rad
if (nygrid/=1) then
call random_number_wrapper(phi)
fp(k,iyp) = xyz0_par(2)+phi*Lxyz_par(2)
endif
if (nzgrid/=1) then
call random_number_wrapper(r)
fp(k,izp)=xyz0_par(3)+r*Lxyz_par(3)
endif
elseif (lspherical_coords) then
if (nxgrid/=1) fp(k,ixp)=rad
if (nygrid/=1) then
call random_number_wrapper(tht)
fp(k,iyp) = xyz0_par(2)+tht*Lxyz_par(2)
endif
if (nzgrid/=1) then
call random_number_wrapper(phi)
fp(k,izp) = xyz0_par(3)+phi*Lxyz_par(3)
endif
endif
!
enddo
!
case ('np-constant')
if (lroot) print*, 'init_particles: Constant number density'
k=1
k_loop: do while (.not. (k>npar_loc))
do l=l1,l2
do m=m1,m2
do n=n1,n2
if (nxgrid/=1) then
call random_number_wrapper(px)
fp(k,ixp)=x(l)+(px-0.5)*dx
endif
if (nygrid/=1) then
call random_number_wrapper(py)
fp(k,iyp)=y(m)+(py-0.5)*dy
endif
if (nzgrid/=1) then
call random_number_wrapper(pz)
fp(k,izp)=z(n)+(pz-0.5)*dz
endif
k=k+1
if (k>npar_loc) exit k_loop
enddo
enddo
enddo
enddo k_loop
!
case ('equidistant')
if (lroot) print*, 'init_particles: Particles placed equidistantly'
dim1=1.0/max(dimensionality,1)
!
! Number of particles per direction. Found by solving the equation system
!
! npar_loc_x/npar_loc_y = Lx_loc/Ly_loc
! npar_loc_x/npar_loc_z = Lx_loc/Lz_loc
! npar_loc_y/npar_loc_z = Ly_loc/Lz_loc
! npar_loc_x*npar_loc_y*npar_loc_z = npar_loc
!
! Found it to be easier to separate in all possible dimensionalities.
! For a missing direction i, set npar_loc_i=1 in the above equations and
! ignore any equation that has Li_loc in it.
!
! Initiate to avoid compiler warnings. Will be overwritten.
npar_loc_x=1;npar_loc_y=1;npar_loc_z=1
!
if (dimensionality==3) then
! 3-D
npar_loc_x=(npar_loc*Lxyz_loc(1)**2/(Lxyz_loc(2)*Lxyz_loc(3)))**dim1
npar_loc_y=(npar_loc*Lxyz_loc(2)**2/(Lxyz_loc(1)*Lxyz_loc(3)))**dim1
npar_loc_z=(npar_loc*Lxyz_loc(3)**2/(Lxyz_loc(1)*Lxyz_loc(2)))**dim1
elseif (dimensionality==2) then
! 2-D
if (nxgrid==1) then
npar_loc_x=1
npar_loc_y=(npar_loc*Lxyz_loc(2)/Lxyz_loc(3))**dim1
npar_loc_z=(npar_loc*Lxyz_loc(3)/Lxyz_loc(2))**dim1
elseif (nygrid==1) then
npar_loc_x=(npar_loc*Lxyz_loc(1)/Lxyz_loc(3))**dim1
npar_loc_y=1
npar_loc_z=(npar_loc*Lxyz_loc(3)/Lxyz_loc(1))**dim1
elseif (nzgrid==1) then
npar_loc_x=(npar_loc*Lxyz_loc(1)/Lxyz_loc(2))**dim1
npar_loc_y=(npar_loc*Lxyz_loc(2)/Lxyz_loc(1))**dim1
npar_loc_z=1
endif
elseif (dimensionality==1) then
! 1-D
if (nxgrid/=1) then
npar_loc_x=npar_loc
npar_loc_y=1
npar_loc_z=1
elseif (nygrid/=1) then
npar_loc_x=1
npar_loc_y=npar_loc
npar_loc_z=1
elseif (nzgrid/=1) then
npar_loc_x=1
npar_loc_y=1
npar_loc_z=npar_loc
endif
endif
! Distance between particles.
dx_par=Lxyz_loc(1)/npar_loc_x
dy_par=Lxyz_loc(2)/npar_loc_y
dz_par=Lxyz_loc(3)/npar_loc_z
! Place first particle.
fp(1,ixp) = x(l1) ; fp(1,iyp) = y(m1) ; fp(1,izp) = z(n1)
if (nxgrid/=1) fp(1,ixp) = xyz0_loc(1)+dx_par/2
if (nygrid/=1) fp(1,iyp) = xyz0_loc(2)+dy_par/2
if (nzgrid/=1) fp(1,izp) = xyz0_loc(3)+dz_par/2
! Place all other particles iteratively.
if (dimensionality==3) then
! 3-D
do k=2,npar_loc
fp(k,ixp)=fp(k-1,ixp)+dx_par
fp(k,iyp)=fp(k-1,iyp)
fp(k,izp)=fp(k-1,izp)
if (fp(k,ixp)>xyz1_loc(1)) then
fp(k,ixp)=fp(1,ixp)
fp(k,iyp)=fp(k,iyp)+dy_par
endif
if (fp(k,iyp)>xyz1_loc(2)) then
fp(k,iyp)=fp(1,iyp)
fp(k,izp)=fp(k,izp)+dz_par
endif
enddo
elseif (dimensionality==2) then
! 2-D
if (nxgrid==1) then
do k=2,npar_loc
fp(k,ixp)=fp(k-1,ixp)
fp(k,iyp)=fp(k-1,iyp)+dy_par
fp(k,izp)=fp(k-1,izp)
if (fp(k,iyp)>xyz1_loc(2)) then
fp(k,iyp)=fp(1,iyp)
fp(k,izp)=fp(k,izp)+dz_par
endif
enddo
elseif (nygrid==1) then
do k=2,npar_loc
fp(k,ixp)=fp(k-1,ixp)+dx_par
fp(k,iyp)=fp(k-1,iyp)
fp(k,izp)=fp(k-1,izp)
if (fp(k,ixp)>xyz1_loc(1)) then
fp(k,ixp)=fp(1,ixp)
fp(k,izp)=fp(k,izp)+dz_par
endif
enddo
elseif (nzgrid==1) then
do k=2,npar_loc
fp(k,ixp)=fp(k-1,ixp)+dx_par
fp(k,iyp)=fp(k-1,iyp)
fp(k,izp)=fp(k-1,izp)
if (fp(k,ixp)>xyz1_loc(1)) then
fp(k,ixp)=fp(1,ixp)
fp(k,iyp)=fp(k,iyp)+dy_par
endif
enddo
endif
elseif (dimensionality==1) then
! 1-D
if (nxgrid/=1) then
do k=2,npar_loc
fp(k,ixp)=fp(k-1,ixp)+dx_par
fp(k,iyp)=fp(k-1,iyp)
fp(k,izp)=fp(k-1,izp)
enddo
elseif (nygrid/=1) then
do k=2,npar_loc
fp(k,ixp)=fp(k-1,ixp)
fp(k,iyp)=fp(k-1,iyp)+dy_par
fp(k,izp)=fp(k-1,izp)
enddo
elseif (nzgrid/=1) then
do k=2,npar_loc
fp(k,ixp)=fp(k-1,ixp)
fp(k,iyp)=fp(k-1,iyp)
fp(k,izp)=fp(k-1,izp)+dz_par
enddo
endif
else
! 0-D
fp(2:npar_loc,ixp)=fp(1,ixp)
fp(2:npar_loc,iyp)=fp(1,iyp)
fp(2:npar_loc,izp)=fp(1,izp)
endif
lequidistant=.true.
!
! Shift particle locations slightly so that a mode appears.
!
case ('shift')
if (lroot) print*, 'init_particles: shift particle positions'
if (.not. lequidistant) then
call fatal_error('init_particles','must place particles equidistantly before shifting!')
endif
k2_xxp=kx_xxp**2+ky_xxp**2+kz_xxp**2
if (k2_xxp==0.0) then
call fatal_error('init_particles','kx_xxp=ky_xxp=kz_xxp=0.0 is not allowed!')
endif
do k=1,npar_loc
fp(k,ixp) = fp(k,ixp) - kx_xxp/k2_xxp*amplxxp* &
sin(kx_xxp*fp(k,ixp)+ky_xxp*fp(k,iyp)+kz_xxp*fp(k,izp))
fp(k,iyp) = fp(k,iyp) - ky_xxp/k2_xxp*amplxxp* &
sin(kx_xxp*fp(k,ixp)+ky_xxp*fp(k,iyp)+kz_xxp*fp(k,izp))
fp(k,izp) = fp(k,izp) - kz_xxp/k2_xxp*amplxxp* &
sin(kx_xxp*fp(k,ixp)+ky_xxp*fp(k,iyp)+kz_xxp*fp(k,izp))
enddo
!
! Shift to egg crate mode 2d, cos(x)cos(z)
!
case ('cosxcosz')
if (lroot) print*, 'init_particles: shift particle positions'
if (.not. lequidistant) then
call fatal_error('init_particles','must place particles equidistantly before shifting!')
endif
k2_xxp=kx_xxp**2+kz_xxp**2
if (k2_xxp==0.0) then
call fatal_error('init_particles','kx_xxp=ky_xxp=kz_xxp=0.0 is not allowed!')
endif
do k=1,npar_loc
fp(k,ixp) = fp(k,ixp) - kx_xxp/k2_xxp*amplxxp* &
sin(kx_xxp*fp(k,ixp))*cos(kz_xxp*fp(k,izp))
fp(k,izp) = fp(k,izp) - kz_xxp/k2_xxp*amplxxp* &
sin(kx_xxp*fp(k,ixp))*cos(kz_xxp*fp(k,izp))
enddo
!
! Shift to egg crate mode 2d, sin(x)sin(z)
!
case ('sinxsinz')
if (lroot) print*, 'init_particles: shift particle positions'
if (.not. lequidistant) then
call fatal_error('init_particles','must place particles equidistantly before shifting!')
endif
k2_xxp=kx_xxp**2+kz_xxp**2
if (k2_xxp==0.0) then
call fatal_error('init_particles','kx_xxp=ky_xxp=kz_xxp=0.0 is not allowed!')
endif
do k=1,npar_loc
fp(k,ixp) = fp(k,ixp) + kx_xxp/k2_xxp*amplxxp* &
cos(kx_xxp*fp(k,ixp))*sin(kz_xxp*fp(k,izp))
fp(k,izp) = fp(k,izp) + kz_xxp/k2_xxp*amplxxp* &
cos(kx_xxp*fp(k,ixp))*sin(kz_xxp*fp(k,izp))
enddo
!
case ('gaussian-z')
if (lroot) print*, 'init_particles: Gaussian particle positions'
do k=1,npar_loc
do while (.true.)
if (nxgrid/=1) then
call random_number_wrapper(r)
fp(k,ixp)=r
endif
if (nygrid/=1) then
call random_number_wrapper(r)
fp(k,iyp)=r
endif
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
if ((fp(k,izp)>=xyz0(3)) .and. (fp(k,izp)<=xyz1(3))) exit
enddo
enddo
if (nxgrid/=1) &
fp(1:npar_loc,ixp)=xyz0_par(1)+fp(1:npar_loc,ixp)*Lxyz_par(1)
if (nygrid/=1) &
fp(1:npar_loc,iyp)=xyz0_par(2)+fp(1:npar_loc,iyp)*Lxyz_par(2)
!
case ('gaussian-x')
if (lroot) print*, 'init_particles: Gaussian particle positions'
do k=1,npar_loc
do while (.true.)
if (nygrid/=1) then
call random_number_wrapper(r)
fp(k,iyp)=r
endif
if (nzgrid/=1) then
call random_number_wrapper(r)
fp(k,izp)=r
endif
call random_number_wrapper(r)
call random_number_wrapper(p)
fp(k,ixp)= xp0*sqrt(-2*alog(r))*cos(2*pi*p)
if ((fp(k,ixp)>=xyz0(1)).and.(fp(k,ixp)<=xyz1(1))) exit
enddo
enddo
if (nygrid/=1) &
fp(1:npar_loc,iyp)=xyz0_par(2)+fp(1:npar_loc,iyp)*Lxyz_par(2)
if (nzgrid/=1) &
fp(1:npar_loc,izp)=xyz0_par(3)+fp(1:npar_loc,izp)*Lxyz_par(3)
!
case ('gaussian-z-pure')
if (lroot) print*, 'init_particles: Gaussian particle positions'
do k=1,npar_loc
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) ! generates a random gaussion number*zp0
endif
enddo
!
case ('gaussian-r')
if (lroot) print*, 'init_particles: Gaussian particle positions'
do k=1,npar_loc
call random_number_wrapper(r)
call random_number_wrapper(p)
call random_number_wrapper(q)
fp(k,ixp)= xp0*sqrt(-2*alog(r))*cos(2*pi*p)*cos(2*pi*q)
fp(k,iyp)= yp0*sqrt(-2*alog(r))*cos(2*pi*p)*sin(2*pi*q)
enddo
!
case ('hole')
!
call map_nearest_grid(fp,ineargrid)
call map_xxp_grid(f,fp,ineargrid)
call sort_particles_imn(fp,ineargrid,ipar)
do k=k1_imn(imn_array(m_hole+m1-1,n_hole+n1-1)), &
k2_imn(imn_array(m_hole+m1-1,n_hole+n1-1))
if (ineargrid(k,1)==l_hole+l1-1) then
print*, k
if (nxgrid/=0) fp(k,ixp)=fp(k,ixp)-dx
endif
enddo
!
case ('birthring')
if (birthring_width>tini) then
if (lgaussian_birthring) then
do k=1,npar_loc
call normal_deviate(rr_tmp(k))
enddo
else
call random_number_wrapper(rr_tmp(1:npar_loc))
endif
rr_tmp(1:npar_loc) = birthring_r+rr_tmp(1:npar_loc)*birthring_width
else
rr_tmp(1:npar_loc) = birthring_r
endif
call random_number_wrapper(az_tmp(1:npar_loc))
az_tmp(1:npar_loc) = -pi+az_tmp(1:npar_loc)*2.0*pi
if (lcartesian_coords) then
fp(1:npar_loc,ixp) = rr_tmp(1:npar_loc)*cos(az_tmp(1:npar_loc))
fp(1:npar_loc,iyp) = rr_tmp(1:npar_loc)*sin(az_tmp(1:npar_loc))
fp(1:npar_loc,izp) = 0.0
else
fp(1:npar_loc,ixp) = rr_tmp(1:npar_loc)
if (lcylindrical_coords) then
fp(1:npar_loc,iyp) = az_tmp(1:npar_loc)
fp(1:npar_loc,izp) = 0.0
elseif (lspherical_coords) then
fp(1:npar_loc,iyp) = pi/2.0
fp(1:npar_loc,izp) = az_tmp(1:npar_loc)
endif
endif
if (lroot .and. nzgrid/=0) print*,"Warning, birthring only implemented for 2D"
!
case default
call fatal_error('init_particles','Unknown value initxxp="'//trim(initxxp(j))//'"')
!
endselect
!
enddo ! do j=1,ninit
!
! Interface for user's own initial condition for position
!
if (linitial_condition) call initial_condition_xxp(f,fp)
!
! 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(nghost+1)
if (nygrid==1) fp(1:npar_loc,iyp)=y(nghost+1)
if (nzgrid==1) fp(1:npar_loc,izp)=z(nghost+1)
!
if (init_repeat/=0) call repeated_init(fp,init_repeat)
!
! Redistribute particles among processors (now that positions are determined).
!
call boundconds_particles(fp,ipar)
!
! Map particle position on the grid.
!
call map_nearest_grid(fp,ineargrid)
call map_xxp_grid(f,fp,ineargrid)
!
! Initial particle velocity.
!
do j=1,ninit
!
select case (initvvp(j))
!
case ('nothing')
if (lroot.and.j==1) print*, 'init_particles: No particle velocity set'
!
case ('zero')
if (lroot) print*, 'init_particles: Zero particle velocity'
fp(1:npar_loc,ivpx:ivpz)=0.0
!
case ('zero-shear')
if (lroot) print*, 'init_particles: Zero particle velocity'
fp(1:npar_loc,ivpy)=-Sshear*fp(1:npar_loc,ixp)
fp(1:npar_loc,ivpx)=0.0
fp(1:npar_loc,ivpz)=0.0
!
case ('constant')
if (lroot) print*, 'init_particles: Constant particle velocity'
if (lroot) &
print*, 'init_particles: vpx0, vpy0, vpz0=', vpx0, vpy0, vpz0
if (lcylindrical_coords) then
fp(1:npar_loc,ivpx)=vpx0*cos(fp(k,iyp))+vpy0*sin(fp(k,iyp))
fp(1:npar_loc,ivpy)=vpy0*cos(fp(k,iyp))-vpx0*sin(fp(k,iyp))
fp(1:npar_loc,ivpz)=vpz0
else
fp(1:npar_loc,ivpx)=vpx0
fp(1:npar_loc,ivpy)=vpy0
fp(1:npar_loc,ivpz)=vpz0
endif
!
case ('constant-1')
if (lroot) &
print*, 'init_particles: Particle 1 velocity vx,vy,vz=', &
vpx1, vpy1, vpz1
do k=1,npar_loc
if (ipar(k)==1) then
fp(k,ivpx)=vpx1
fp(k,ivpy)=vpy1
fp(k,ivpz)=vpz1
endif
enddo
!
case ('constant-2')
if (lroot) &
print*, 'init_particles: Particle 2 velocity vx,vy,vz=', &
vpx2, vpy2, vpz2
do k=1,npar_loc
if (ipar(k)==2) then
fp(k,ivpx)=vpx2
fp(k,ivpy)=vpy2
fp(k,ivpz)=vpz2
endif
enddo
!
case ('constant-3')
if (lroot) &
print*, 'init_particles: Particle 3 velocity vx,vy,vz=', &
vpx3, vpy3, vpz3
do k=1,npar_loc
if (ipar(k)==3) then
fp(k,ivpx)=vpx3
fp(k,ivpy)=vpy3
fp(k,ivpz)=vpz3
endif
enddo
!
case ('sinwave-phase')
if (lroot) print*, 'init_particles: sinwave-phase'
if (lroot) &
print*, 'init_particles: vpx0, vpy0, vpz0=', vpx0, vpy0, vpz0
do k=1,npar_loc
fp(k,ivpx)=fp(k,ivpx)+vpx0*sin(kx_vpx*fp(k,ixp)+ky_vpx*fp(k,iyp)+kz_vpx*fp(k,izp)+phase_vpx)
fp(k,ivpy)=fp(k,ivpy)+vpy0*sin(kx_vpy*fp(k,ixp)+ky_vpy*fp(k,iyp)+kz_vpy*fp(k,izp)+phase_vpy)
fp(k,ivpz)=fp(k,ivpz)+vpz0*sin(kx_vpz*fp(k,ixp)+ky_vpz*fp(k,iyp)+kz_vpz*fp(k,izp)+phase_vpz)
enddo
!
case ('coswave-phase')
if (lroot) print*, 'init_particles: coswave-phase'
if (lroot) &
print*, 'init_particles: vpx0, vpy0, vpz0=', vpx0, vpy0, vpz0
do k=1,npar_loc
fp(k,ivpx)=fp(k,ivpx)+vpx0*cos(kx_vpx*fp(k,ixp)+ky_vpx*fp(k,iyp)+kz_vpx*fp(k,izp)+phase_vpx)
fp(k,ivpy)=fp(k,ivpy)+vpy0*cos(kx_vpy*fp(k,ixp)+ky_vpy*fp(k,iyp)+kz_vpy*fp(k,izp)+phase_vpy)
fp(k,ivpz)=fp(k,ivpz)+vpz0*cos(kx_vpz*fp(k,ixp)+ky_vpz*fp(k,iyp)+kz_vpz*fp(k,izp)+phase_vpz)
enddo
!
case ('random')
if (lroot) print*, 'init_particles: Random particle velocities; '// &
'delta_vp0=', delta_vp0
do k=1,npar_loc
call random_number_wrapper(r)
fp(k,ivpx) = fp(k,ivpx) + delta_vp0*(2*r-1)
call random_number_wrapper(r)
fp(k,ivpy) = fp(k,ivpy) + delta_vp0*(2*r-1)
call random_number_wrapper(r)
fp(k,ivpz) = fp(k,ivpz) + delta_vp0*(2*r-1)
enddo
!
case ('random-x')
if (lroot) print*, 'init_particles: Random particle x-velocity; '// &
'delta_vp0=', delta_vp0
do k=1,npar_loc
call random_number_wrapper(r)
fp(k,ivpx) = fp(k,ivpx) + delta_vp0*(2*r-1)
enddo
!
case ('random-y')
if (lroot) print*, 'init_particles: Random particle y-velocity; '// &
'delta_vp0=', delta_vp0
do k=1,npar_loc
call random_number_wrapper(r)
fp(k,ivpy) = fp(k,ivpy) + delta_vp0*(2*r-1)
enddo
!
case ('random-z')
if (lroot) print*, 'init_particles: Random particle z-velocity; '// &
'delta_vp0=', delta_vp0
do k=1,npar_loc
call random_number_wrapper(r)
fp(k,ivpz) = fp(k,ivpz) + delta_vp0*(2*r-1)
enddo
!
case ('average-to-zero')
call mpireduce_sum(sum(fp(1:npar_loc,ivpx)),vpx_sum)
call mpireduce_sum(sum(fp(1:npar_loc,ivpy)),vpy_sum)
call mpireduce_sum(sum(fp(1:npar_loc,ivpz)),vpz_sum)
call mpibcast_real(vpx_sum)
call mpibcast_real(vpy_sum)
call mpibcast_real(vpz_sum)
fp(1:npar_loc,ivpx)=fp(1:npar_loc,ivpx)-vpx_sum/npar
fp(1:npar_loc,ivpy)=fp(1:npar_loc,ivpy)-vpy_sum/npar
fp(1:npar_loc,ivpz)=fp(1:npar_loc,ivpz)-vpz_sum/npar
!
case ('follow-gas')
if (lroot) &
print*, 'init_particles: Particle velocity equal to gas velocity'
do k=1,npar_loc
call interpolate_linear(f,iux,iuz,fp(k,ixp:izp),uup, &
ineargrid(k,:),0,0)
fp(k,ivpx:ivpz) = uup
enddo
!
case ('jeans-wave-dustpar-x')
! assumes rhs_poisson_const=1 !
do k=1,npar_loc
fp(k,ivpx) = fp(k,ivpx) - amplxxp* &
(sqrt(1+4*1.0*1.0*tausp**2)-1)/ &
(2*kx_xxp*1.0*tausp)*sin(kx_xxp*(fp(k,ixp)))
enddo
!
case ('dragforce_equilibrium','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
cs=sqrt(cs20)
!
if (ldragforce_equi_noback) eps=0.0
!
if (lroot) print*, 'init_particles: average dust-to-gas ratio=', eps
! Set gas velocity field.
if (lhydro) then
do l=l1,l2; do m=m1,m2; do n=n1,n2
! Take either global or local dust-to-gas ratio.
if (.not. ldragforce_equi_global_eps) eps = f(l,m,n,irhop) / get_gas_density(f,l,m,n)
!
f(l,m,n,iux) = f(l,m,n,iux) - &
beta_glnrho_global(1)*eps*Omega*tausp/ &
((1.0+eps)**2+(Omega*tausp)**2)*cs
f(l,m,n,iuy) = f(l,m,n,iuy) + &
beta_glnrho_global(1)*(1+eps+(Omega*tausp)**2)/ &
(2*((1.0+eps)**2+(Omega*tausp)**2))*cs
!
enddo; enddo; enddo
endif
! Set particle velocity field.
do k=1,npar_loc
! Take either global or local dust-to-gas ratio.
if (ldragforce_equi_noback) then
eps=0.0
else
if (.not. ldragforce_equi_global_eps) then
ix0=ineargrid(k,1); iy0=ineargrid(k,2); iz0=ineargrid(k,3)
eps = f(ix0,iy0,iz0,irhop) / get_gas_density(f,ix0,iy0,iz0)
endif
endif
!
fp(k,ivpx) = fp(k,ivpx) + &
beta_glnrho_global(1)*Omega*tausp/ &
((1.0+eps)**2+(Omega*tausp)**2)*cs
fp(k,ivpy) = fp(k,ivpy) + &
beta_glnrho_global(1)*(1+eps)/ &
(2*((1.0+eps)**2+(Omega*tausp)**2))*cs
!
enddo
!
case ('dragforce_equi_nohydro')
!
do k=1,npar_loc
fp(k,ivpx) = fp(k,ivpx) - 2*Deltauy_gas_friction* &
1/(1.0/(Omega*tausp)+Omega*tausp)
fp(k,ivpy) = fp(k,ivpy) - Deltauy_gas_friction* &
1/(1.0+(Omega*tausp)**2)
enddo
!
case ('dragforce_equi_dust')
!
! Equilibrium between drag force and Coriolis force on the dust.
!
if (lroot) then
print*, 'init_particles: drag equilibrium dust'
print*, 'init_particles: beta_dPdr_dust=', beta_dPdr_dust
endif
! Set particle velocity field.
cs=sqrt(cs20)
do k=1,npar_loc
fp(k,ivpx) = fp(k,ivpx) + &
1/(Omega*tausp+1/(Omega*tausp))*beta_dPdr_dust*cs
fp(k,ivpy) = fp(k,ivpy) - &
1/(1.0+1/(Omega*tausp)**2)*beta_dPdr_dust/2*cs
enddo
!
case ('Keplerian','keplerian')
!
! Keplerian velocity based on gravr.
!
if (lroot) then
print*, 'init_particles: Keplerian velocity'
if (lshear) call fatal_error("init_particles",&
"Keplerian initial condition is for global disks, not shearing boxes")
endif
do k=1,npar_loc
if (lcartesian_coords) then
rad=sqrt(fp(k,ixp)**2+fp(k,iyp)**2+fp(k,izp)**2)
OO=sqrt(gravr)*rad**(-1.5)
fp(k,ivpx) = -OO*fp(k,iyp)
fp(k,ivpy) = OO*fp(k,ixp)
fp(k,ivpz) = 0.0
elseif (lcylindrical_coords) then
rad=fp(k,ixp)
OO=sqrt(gravr)*rad**(-1.5)
fp(k,ivpx) = 0.0
fp(k,ivpy) = OO*rad
fp(k,ivpz) = 0.0
elseif (lspherical_coords) then
rad=fp(k,ixp)*sin(fp(k,iyp))
OO=sqrt(gravr)*rad**(-1.5)
fp(k,ivpx) = 0.0
fp(k,ivpy) = 0.0
fp(k,ivpz) = OO*rad
endif
enddo
!
! Explosion.
!
case ('explosion')
do k=1,npar_loc
rad=sqrt(fp(k,ixp)**2+fp(k,iyp)**2+fp(k,izp)**2)
fp(k,ivpx) = delta_vp0*fp(k,ixp)/rp_ext
fp(k,ivpy) = delta_vp0*fp(k,iyp)/rp_ext
fp(k,ivpz) = delta_vp0*fp(k,izp)/rp_ext
enddo
!
!
case default
call fatal_error('init_particles','Unknown value initvvp="'//trim(initvvp(j))//'"')
!
endselect
!
enddo ! do j=1,ninit
!
! Interface for user's own initial condition
!
if (linitial_condition) call initial_condition_vvp(f,fp)
!
! Map particle velocity on the grid.
!
call map_vvp_grid(f,fp,ineargrid)
!
! Sort particles (must happen at the end of the subroutine so that random
! positions and velocities are not displaced relative to when there is no
! sorting).
!
call sort_particles_imn(fp,ineargrid,ipar)
!
! Set the initial auxiliary array for the passive scalar to zero
!
if (ldiffuse_passive) f(:,:,:,idlncc) = 0.0
if (ldiffuse_dragf) f(:,:,:,idfx:idfz) = 0.0
!
endsubroutine init_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), intent (inout) :: fp
integer, dimension (mpar_loc,3), intent (inout) :: 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, normal_deviate
use Particles_diagnos_state, only: insert_particles_diagnos_state
use SharedVariables, only: get_shared_variable
use Mpicomm, only: mpireduce_sum_int
use Particles_number, only: set_particle_number
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mpar_loc,mparray), intent (inout) :: fp
integer, dimension (mpar_loc,3), intent (inout) :: ineargrid
real, dimension (mpar_loc) :: rr_tmp, az_tmp, OO_tmp
!
logical, save :: linsertmore=.true.
real :: xx0, yy0, r2, r, tmp
integer :: j, k, n_insert, npar_loc_old, iii, particles_insert_rate_tmp
real, pointer :: gravr
!
! Insertion of particles is stopped when maximum number of particles is reached,
! unless linsert_as_many_as_possible is set.
! Maximum numer of particles allowed in system is defined by max_particles,
! initialized to npar. Note that this may cause errors at a processor further
! downstream, if particles accumulate and mpar_loc is too small.
! Since root inserts all new particles, make sure
! npar_loc + n_insert < mpar_loc
! so that a processor can not exceed its maximum number of particles.
!
if (t>tstart_insert_particles+particles_insert_ramp_time) then
particles_insert_rate_tmp = particles_insert_rate
else
particles_insert_rate_tmp = int(real(particles_insert_rate)*(t-tstart_insert_particles)/particles_insert_ramp_time)
endif
call mpireduce_sum_int(npar_loc,npar_total)
!
if (lroot) then
avg_n_insert=particles_insert_rate_tmp*dt
n_insert=int(avg_n_insert + remaining_particles)
! Remaining particles saved for subsequent timestep:
remaining_particles=avg_n_insert + remaining_particles - n_insert
!
! Insert particles if maximum count is not reached
!
if ((n_insert+npar_loc <= mpar_loc) .and. (npar_inserted_tot+n_insert <= max_particles) &
.and. (t<max_particle_insert_time) .and. (t>tstart_insert_particles)) then
linsertmore=.true.
else
linsertmore=.false.
endif
!
! Continue inserting even if max_particles is reached, if linsert_as_many_as_possible
! is set.
!
if (linsert_as_many_as_possible) then
n_insert=min((mpar_loc-npar_loc),n_insert)
if ((n_insert+npar_loc <= mpar_loc) &
.and. (t<max_particle_insert_time) .and. (t>tstart_insert_particles)) then
linsertmore=.true.
endif
endif
!
if (linsertmore) then
! Actual (integer) number of particles to be inserted at this timestep:
do iii=npar_loc+1,npar_loc+n_insert
ipar(iii)=npar_inserted_tot+iii-npar_loc
enddo
npar_loc_old=npar_loc
npar_loc=npar_loc + n_insert
!
! Update total number of inserted particles, npar_inserted_tot.
! Not the same as npar_total, which is the number of particles in the system,
! without couting removed particles
!
npar_inserted_tot = n_insert + npar_inserted_tot
!
! Insert particles in chosen position (as in init_particles).
!
do j=1,ninit
select case (initxxp(j))
case ('random-box')
!
do k=npar_loc_old+1,npar_loc
if (nxgrid/=1) then
call random_number_wrapper(r)
fp(k,ixp)=r
endif
if (nygrid/=1) then
call random_number_wrapper(r)
fp(k,iyp)=r
endif
if (nzgrid/=1) then
call random_number_wrapper(r)
fp(k,izp)=r
endif
if (lcylindrical_coords) then
xx0=xp0+fp(k,ixp)*Lx0
yy0=yp0+fp(k,iyp)*Ly0
r2=xx0**2+yy0**2
if (nxgrid/=1) fp(k,ixp)=sqrt(r2)
if (nygrid/=1) fp(k,iyp)=atan(yy0/xx0)+pi*(xx0/abs(xx0)-1)*0.5
if (nzgrid/=1) fp(k,izp)=zp0+fp(k,izp)*Lz0
else
if (nxgrid/=1) fp(k,ixp)=xp0+fp(k,ixp)*Lx0
if (nygrid/=1) fp(k,iyp)=yp0+fp(k,iyp)*Ly0
if (nzgrid/=1) fp(k,izp)=zp0+fp(k,izp)*Lz0
endif
enddo
!
! Maybe random-cylindrical case should be combined with normal initxxp case
!
case ('random-cylindrical','random-cyl')
if (lcylindrical_coords.or.lcartesian_coords) then
tmp=2-dustdensity_powerlaw
elseif (lspherical_coords) then
tmp=3-dustdensity_powerlaw
else
call fatal_error("init_particles",&
"The world is flat, and we never got here")
endif
!
call random_number_wrapper(rr_tmp(npar_loc_old+1:npar_loc))
rr_tmp(npar_loc_old+1:npar_loc) = rp_int**tmp + &
rr_tmp(npar_loc_old+1:npar_loc)*(rp_ext**tmp-rp_int**tmp)
rr_tmp(npar_loc_old+1:npar_loc) = rr_tmp(npar_loc_old+1:npar_loc)**(1./tmp)
if ((lcartesian_coords) .or. (lcylindrical_coords .and. nygrid/=1) .or. (lspherical_coords .and. nzgrid/=1)) then
call random_number_wrapper(az_tmp(npar_loc_old+1:npar_loc))
az_tmp(npar_loc_old+1:npar_loc) = -pi + 2.0*pi*az_tmp(npar_loc_old+1:npar_loc)
endif
if ((lcartesian_coords) .or. (lcylindrical_coords .and. nzgrid/=1) .or. (lspherical_coords .and. nygrid/=1)) then
call random_number_wrapper(fp(npar_loc_old+1:npar_loc,izp))
endif
if (lcartesian_coords) then
fp(npar_loc_old+1:npar_loc,iyp) = az_tmp(npar_loc_old+1:npar_loc)*2.0*pi
if (nxgrid/=1) fp(npar_loc_old+1:npar_loc,ixp) = rr_tmp(npar_loc_old+1:npar_loc) &
*cos(az_tmp(npar_loc_old+1:npar_loc))
if (nygrid/=1) fp(npar_loc_old+1:npar_loc,iyp) = rr_tmp(npar_loc_old+1:npar_loc) &
*sin(az_tmp(npar_loc_old+1:npar_loc))
if (nzgrid/=1) fp(npar_loc_old+1:npar_loc,izp) = xyz0(3)+fp(npar_loc_old+1:npar_loc,izp)*Lxyz(3)
elseif (lcylindrical_coords) then
if (nxgrid/=1) fp(npar_loc_old+1:npar_loc,ixp) = rr_tmp(npar_loc_old+1:npar_loc)
if (nygrid/=1) fp(npar_loc_old+1:npar_loc,iyp) = xyz0(2)+az_tmp(npar_loc_old+1:npar_loc)*Lxyz(2)
if (nzgrid/=1) fp(npar_loc_old+1:npar_loc,izp) = xyz0(3)+fp(npar_loc_old+1:npar_loc,izp)*Lxyz(3)
elseif (lspherical_coords) then
if (nxgrid/=1) fp(npar_loc_old+1:npar_loc,ixp) = rr_tmp(npar_loc_old+1:npar_loc)
if (nygrid/=1) fp(npar_loc_old+1:npar_loc,iyp) = xyz0(2)+az_tmp(npar_loc_old+1:npar_loc)*Lxyz(2)
if (nzgrid/=1) fp(npar_loc_old+1:npar_loc,izp) = xyz0(3)+fp(npar_loc_old+1:npar_loc,izp)*Lxyz(3)
endif
!
case ('birthring')
if (birthring_width>tini) then
if (lgaussian_birthring) then
do k=npar_loc_old+1,npar_loc
call normal_deviate(rr_tmp(k))
enddo
else
call random_number_wrapper(rr_tmp(npar_loc_old+1:npar_loc))
endif
rr_tmp(npar_loc_old+1:npar_loc) = birthring_r+rr_tmp(npar_loc_old+1:npar_loc)*birthring_width
else
rr_tmp(npar_loc_old+1:npar_loc) = birthring_r
endif
call random_number_wrapper(az_tmp(npar_loc_old+1:npar_loc))
az_tmp(npar_loc_old+1:npar_loc) = -pi+az_tmp(npar_loc_old+1:npar_loc)*2.0*pi
if (lcartesian_coords) then
fp(npar_loc_old+1:npar_loc,ixp) = rr_tmp(npar_loc_old+1:npar_loc)*cos(az_tmp(npar_loc_old+1:npar_loc))
fp(npar_loc_old+1:npar_loc,iyp) = rr_tmp(npar_loc_old+1:npar_loc)*sin(az_tmp(npar_loc_old+1:npar_loc))
fp(npar_loc_old+1:npar_loc,izp) = 0.0
else
fp(npar_loc_old+1:npar_loc,ixp) = rr_tmp(npar_loc_old+1:npar_loc)
if (lcylindrical_coords) then
fp(npar_loc_old+1:npar_loc,iyp) = az_tmp(npar_loc_old+1:npar_loc)
fp(npar_loc_old+1:npar_loc,izp) = 0.0
elseif (lspherical_coords) then
fp(npar_loc_old+1:npar_loc,iyp) = pi/2.0
fp(npar_loc_old+1:npar_loc,izp) = az_tmp(npar_loc_old+1:npar_loc)
endif
endif
!
case ('nothing')
if (j==1) print*, 'init_particles: nothing'
!
case default
call fatal_error_local('init_particles','Unknown value initxxp="'//trim(initxxp(j))//'"')
!
endselect
enddo
!
! Initial particle velocity.
!
do j=1,ninit
select case (initvvp(j))
case ('nothing')
if (j==1) print*, 'init_particles: No particle velocity set'
!
case ('Keplerian','keplerian')
call get_shared_variable('gravr',gravr)
if (lcartesian_coords) then
rr_tmp(npar_loc_old+1:npar_loc)=sqrt(fp(npar_loc_old+1:npar_loc,ixp)**2+ &
fp(npar_loc_old+1:npar_loc,iyp)**2+fp(npar_loc_old+1:npar_loc,izp)**2)
OO_tmp(npar_loc_old+1:npar_loc)=sqrt(gravr)*rr_tmp(npar_loc_old+1:npar_loc)**(-1.5)
fp(npar_loc_old+1:npar_loc,ivpx) = -OO_tmp(npar_loc_old+1:npar_loc)*fp(npar_loc_old+1:npar_loc,iyp)
fp(npar_loc_old+1:npar_loc,ivpy) = OO_tmp(npar_loc_old+1:npar_loc)*fp(npar_loc_old+1:npar_loc,ixp)
fp(npar_loc_old+1:npar_loc,ivpz) = 0.0
elseif (lcylindrical_coords) then
rr_tmp(npar_loc_old+1:npar_loc)=fp(npar_loc_old+1:npar_loc,ixp)
OO_tmp(npar_loc_old+1:npar_loc)=sqrt(gravr)*rr_tmp(npar_loc_old+1:npar_loc)**(-1.5)
fp(npar_loc_old+1:npar_loc,ivpx) = 0.0
fp(npar_loc_old+1:npar_loc,ivpy) = OO_tmp(npar_loc_old+1:npar_loc)*rr_tmp(npar_loc_old+1:npar_loc)
fp(npar_loc_old+1:npar_loc,ivpz) = 0.0
elseif (lspherical_coords) then
rr_tmp(npar_loc_old+1:npar_loc)=fp(npar_loc_old+1:npar_loc,ixp)*sin(fp(npar_loc_old+1:npar_loc,iyp))
OO_tmp(npar_loc_old+1:npar_loc)=sqrt(gravr)*rr_tmp(npar_loc_old+1:npar_loc)**(-1.5)
fp(npar_loc_old+1:npar_loc,ivpx) = 0.0
fp(npar_loc_old+1:npar_loc,ivpy) = 0.0
fp(npar_loc_old+1:npar_loc,ivpz) = OO_tmp(npar_loc_old+1:npar_loc)*rr_tmp(npar_loc_old+1:npar_loc)
endif
!
case ('constant')
if (lcylindrical_coords) then
fp(npar_loc_old+1:npar_loc,ivpx)&
=vpx0*cos(fp(npar_loc_old+1:npar_loc,iyp))&
+vpy0*sin(fp(npar_loc_old+1:npar_loc,iyp))
fp(npar_loc_old+1:npar_loc,ivpy)&
=vpy0*cos(fp(npar_loc_old+1:npar_loc,iyp))&
-vpx0*sin(fp(npar_loc_old+1:npar_loc,iyp))
fp(npar_loc_old+1:npar_loc,ivpz)=vpz0
else
fp(npar_loc_old+1:npar_loc,ivpx)=vpx0
fp(npar_loc_old+1:npar_loc,ivpy)=vpy0
fp(npar_loc_old+1:npar_loc,ivpz)=vpz0
endif
!
case default
call fatal_error_local('init_particles','Unknown value initvvp="'//trim(initvvp(j))//'"')
!
endselect
!
enddo ! do j=1,ninit
!
! Initialize particle radius
!
if (lparticles_radius) call set_particle_radius(f,fp,npar_loc_old+1,npar_loc)
if (lparticles_number) call set_particle_number(f,fp,npar_loc_old+1,npar_loc)
if (lbirthring_depletion) fp(npar_loc_old+1:npar_loc,ibrtime) = 0.0
!
! Particles are not allowed to be present in non-existing dimensions.
! This would give huge problems with interpolation later.
!
if (nxgrid==1) fp(npar_loc_old+1:npar_loc,ixp)=x(nghost+1)
if (nygrid==1) fp(npar_loc_old+1:npar_loc,iyp)=y(nghost+1)
if (nzgrid==1) fp(npar_loc_old+1:npar_loc,izp)=z(nghost+1)
!
if (lparticles_diagnos_state) &
call insert_particles_diagnos_state(fp, npar_loc_old)
!
endif
endif ! if (lroot) then
!
! Redistribute particles only when t < max_particle_insert_time
! and t>tstart_insert_particles.
! Could have included some other tests here aswell......
!
if (t<max_particle_insert_time .and. t>tstart_insert_particles) then
!
! Redistribute particles among processors.
!
call boundconds_particles(fp,ipar,linsert=.true.)
!
! Map particle position on the grid.
!
call map_nearest_grid(fp,ineargrid)
call map_xxp_grid(f,fp,ineargrid)
!
! Map particle velocity on the grid.
!
call map_vvp_grid(f,fp,ineargrid)
!
! Sort particles (must happen at the end of the subroutine so that random
! positions and velocities are not displaced relative to when there is no
! sorting).
!
call sort_particles_imn(fp,ineargrid,ipar)
endif
!
if (lbirthring_depletion) then
if (lcartesian_coords) then
rr_tmp(1:npar_loc) = sqrt(fp(1:npar_loc,ixp)**2.0+fp(1:npar_loc,iyp)**2.0)
else
rr_tmp(1:npar_loc) = fp(1:npar_loc,ixp)
endif
do k=1,npar_loc
if (lgaussian_birthring) then
fp(k,ibrtime) = fp(k,ibrtime) + &
dt*exp(-((rr_tmp(k)-birthring_r)**2.0)/(2.0*birthring_width**2.0))
else
if ((rr_tmp(k) .ge. birthring_r-birthring_width) .and. &
(rr_tmp(k) .le. birthring_r+birthring_width)) &
fp(k,ibrtime) = fp(k,ibrtime)+dt
endif
if (fp(k,ibrtime) .ge. birthring_lifetime) &
call remove_particle(fp,ipar,k)
enddo
endif
!
endsubroutine insert_particles
!***********************************************************************
subroutine particles_dragforce_stiff(f,fp,ineargrid)
!
! Force stiff drag force equations towards their equilibrium.
!
! 13-july-17/nils: coded
!
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mpar_loc,mparray) :: fp
integer, dimension (mpar_loc,3) :: ineargrid
!
if (ldragforce_stiff) then
call fatal_error("particles_dragforce_stiff","not implemented")
endif
!NILS: Remove all reference to dragforce stiff - or put it back in again???
!NILS: Email Dhruba about this.
!
endsubroutine particles_dragforce_stiff
!***********************************************************************
subroutine pencil_criteria_particles
!
! All pencils that the Particles module depends on are specified here.
!
! 20-04-06/anders: coded
!
if (ldragforce_gas_par) then
lpenc_requested(i_epsp)=.true.
lpenc_requested(i_np)=.true.
lpenc_requested(i_rho1) = .true.
endif
if (ldraglaw_epstein) then
lpenc_requested(i_cs2)=.true.
lpenc_requested(i_rho)=.true.
endif
if (ldragforce_heat .or. lcollisional_heat) then
lpenc_requested(i_TT1)=.true.
lpenc_requested(i_rho1)=.true.
endif
if (lcollisional_cooling_taucool) then
lpenc_requested(i_np)=.true.
endif
if (lcollisional_cooling_rms) then
lpenc_requested(i_epsp)=.true.
endif
if (lcollisional_cooling_rms .or. lcollisional_dragforce_cooling) then
lpenc_requested(i_np)=.true.
lpenc_requested(i_rho1)=.true.
endif
if (ldraglaw_epstein_transonic .or.&
ldraglaw_eps_stk_transonic) then
lpenc_requested(i_uu)=.true.
lpenc_requested(i_rho)=.true.
lpenc_requested(i_cs2)=.true.
endif
if (ldragforce_stiff) then
lpenc_requested(i_fpres)=.true.
lpenc_requested(i_jxbr)=.true.
lpenc_requested(i_fvisc)=.true.
endif
if (lthermophoretic_forces) then
lpenc_requested(i_gTT)=.true.
endif
!
if (idiag_npm/=0 .or. idiag_np2m/=0 .or. idiag_npmax/=0 .or. &
idiag_npmin/=0 .or. idiag_npmx/=0 .or. idiag_npmy/=0 .or. &
idiag_npmz/=0 .or. idiag_nparpmax/=0) lpenc_diagnos(i_np)=.true.
if (idiag_rhopm/=0 .or. idiag_rhoprms/=0 .or. idiag_rhop2m/=0 .or. &
idiag_rhopmax/=0 .or. idiag_rhopmin/=0 .or. idiag_rhopmphi/=0 .or. &
idiag_rhopmx/=0 .or. idiag_rhopmy/=0 .or. idiag_rhopmz/=0) &
lpenc_diagnos(i_rhop)=.true.
if (idiag_rhop2mx /= 0 .or. idiag_rhop2my /= 0 .or. idiag_rhop2mz /= 0) lpenc_diagnos(i_rhop) = .true.
if (idiag_dedragp/=0 .or. idiag_decollp/=0) then
lpenc_diagnos(i_TT1)=.true.
lpenc_diagnos(i_rho1)=.true.
endif
if (idiag_epspmx/=0 .or. idiag_epspmy/=0 .or. idiag_epspmz/=0 .or. &
idiag_epspmin/=0 .or. idiag_epspmax/=0 .or. idiag_epspm/=0) &
lpenc_diagnos(i_epsp)=.true.
if (idiag_rhopmxy/=0 .or. idiag_rhopmxz/=0 .or. idiag_rhopmphi/=0) &
lpenc_diagnos2d(i_rhop)=.true.
if (idiag_npmxy/=0 ) lpenc_diagnos2d(i_np)=.true.
if (idiag_sigmap /= 0) lpenc_diagnos2d(i_rhop) = .true.
!
if (maxval(idiag_npvzmz) > 0) lpenc_requested(i_npvz)=.true.
if (maxval(idiag_npvz2mz) > 0) lpenc_requested(i_npvz2)=.true.
if (maxval(idiag_npuzmz) > 0) lpenc_requested(i_npuz)=.true.
if (maxval(idiag_nptz) > 0) lpenc_requested(i_np_rad)=.true.
if (idiag_Shm /= 0) lpenc_requested(i_sherwood)=.true.
!
endsubroutine pencil_criteria_particles
!***********************************************************************
subroutine pencil_interdep_particles(lpencil_in)
!
! Interdependency among pencils provided by the Particles module
! is specified here.
!
! 16-feb-06/anders: dummy
!
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_rhop)=.true.
lpencil_in(i_rho1)=.true.
endif
!
endsubroutine pencil_interdep_particles
!***********************************************************************
subroutine calc_pencils_particles(f,p)
!
use Sub, only: grad
!
! Calculate Particles pencils.
! Most basic pencils should come first, as others may depend on them.
!
! 16-feb-06/anders: coded
!
real, dimension (mx,my,mz,mfarray) :: f
type (pencil_case) :: p
!
if (lpencil(i_np)) then
if (inp/=0) then
p%np=f(l1:l2,m,n,inp)
else
p%np=0.0
endif
endif
!
if (ipeh>0) p%peh=f(l1:l2,m,n,ipeh)
!
endsubroutine calc_pencils_particles
!***********************************************************************
subroutine dxxp_dt(f,df,fp,dfp,ineargrid)
!
! Evolution of dust particle position.
!
! 02-jan-05/anders: coded
!
real, dimension (mx,my,mz,mfarray), intent (in) :: f
real, dimension (mx,my,mz,mvar), intent (inout) :: df
real, dimension (mpar_loc,mparray), intent (in) :: fp
real, dimension (mpar_loc,mpvar), intent (inout) :: dfp
integer, dimension (mpar_loc,3), intent (in) :: ineargrid
!
logical :: lheader, lfirstcall=.true.
!
! Print out header information in first time step.
!
lheader=lfirstcall .and. lroot
!
! Identify module and boundary conditions.
!
if (lheader) print*,'dxxp_dt: Calculate dxxp_dt'
if (lheader) then
print*, 'dxxp_dt: Particles boundary condition bcpx=', bcpx
print*, 'dxxp_dt: Particles boundary condition bcpy=', bcpy
print*, 'dxxp_dt: Particles boundary condition bcpz=', bcpz
endif
!
if (lheader) print*, 'dxxp_dt: Set rate of change of particle '// &
'position equal to particle velocity.'
!
! The rate of change of a particle's position is the particle's velocity.
! If pointmasses are used do the evolution in that module, for better
! conservation of the Jacobi constant.
!
if (.not.lpointmasses) then
if (lcartesian_coords) then
!
if (nxgrid/=1) &
dfp(1:npar_loc,ixp) = dfp(1:npar_loc,ixp) + fp(1:npar_loc,ivpx)
if (nygrid/=1) &
dfp(1:npar_loc,iyp) = dfp(1:npar_loc,iyp) + fp(1:npar_loc,ivpy)
if (nzgrid/=1) &
dfp(1:npar_loc,izp) = dfp(1:npar_loc,izp) + fp(1:npar_loc,ivpz)
!
elseif (lcylindrical_coords) then
if (nxgrid/=1) &
dfp(1:npar_loc,ixp) = dfp(1:npar_loc,ixp) + fp(1:npar_loc,ivpx)
if (nygrid/=1) &
dfp(1:npar_loc,iyp) = dfp(1:npar_loc,iyp) + &
fp(1:npar_loc,ivpy)/max(fp(1:npar_loc,ixp),tini)
if (nzgrid/=1) &
dfp(1:npar_loc,izp) = dfp(1:npar_loc,izp) + fp(1:npar_loc,ivpz)
!
elseif (lspherical_coords) then
!
if (nxgrid/=1) &
dfp(1:npar_loc,ixp) = dfp(1:npar_loc,ixp) + fp(1:npar_loc,ivpx)
if (nygrid/=1) &
dfp(1:npar_loc,iyp) = dfp(1:npar_loc,iyp) + &
fp(1:npar_loc,ivpy)/max(fp(1:npar_loc,ixp),tini)
if (nzgrid/=1) &
dfp(1:npar_loc,izp) = dfp(1:npar_loc,izp) + &
fp(1:npar_loc,ivpz)/(max(fp(1:npar_loc,ixp),tini)*sin(fp(1:npar_loc,iyp)))
endif
endif
!
! With shear there is an extra term due to the background shear flow.
!
if (lshear.and.nygrid/=1) dfp(1:npar_loc,iyp) = &
dfp(1:npar_loc,iyp) - qshear*Omega*fp(1:npar_loc,ixp)
!
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,fp,dfp,ineargrid)
!
! Evolution of dust particle velocity.
!
! 29-dec-04/anders: coded
!
use Diagnostics
use EquationOfState, only: cs20
!
real, dimension (mx,my,mz,mfarray), intent (inout) :: f
real, dimension (mx,my,mz,mvar), intent (inout) :: df
real, dimension (mpar_loc,mparray), intent (in) :: fp
real, dimension (mpar_loc,mpvar), intent (inout) :: dfp
integer, dimension (mpar_loc,3), intent (in) :: ineargrid
!
real :: Omega2
integer :: npar_found
logical :: lheader, lfirstcall=.true.
!
! Print out header information in first time step.
!
lheader=lfirstcall .and. lroot
if (lheader) then
print*,'dvvp_dt: Calculate dvvp_dt'
endif
!
! Add Coriolis force from rotating coordinate frame.
!
if (Omega/=0.) then
if (lcoriolis_force_par) then
if (lheader) print*,'dvvp_dt: Add Coriolis force; Omega=', Omega
Omega2=2*Omega
if (.not.lspherical_coords) then
dfp(1:npar_loc,ivpx) = dfp(1:npar_loc,ivpx) + &
Omega2*fp(1:npar_loc,ivpy)
dfp(1:npar_loc,ivpy) = dfp(1:npar_loc,ivpy) - &
Omega2*fp(1:npar_loc,ivpx)
else
call fatal_error('dvvp_dt', &
'Coriolis force on the particles is not yet implemented for spherical coordinates.')
endif
endif
!
! Add centrifugal force.
!
if (lcentrifugal_force_par) then
if (lheader) print*,'dvvp_dt: Add Centrifugal force; Omega=', Omega
if (lcartesian_coords) then
!
dfp(1:npar_loc,ivpx) = dfp(1:npar_loc,ivpx) + &
Omega**2*fp(1:npar_loc,ixp)
!
dfp(1:npar_loc,ivpy) = dfp(1:npar_loc,ivpy) + &
Omega**2*fp(1:npar_loc,iyp)
!
elseif (lcylindrical_coords) then
dfp(1:npar_loc,ivpx) = &
dfp(1:npar_loc,ivpx) + Omega**2*fp(1:npar_loc,ixp)
else
call fatal_error('dvvp_dt', &
'Centrifugal force on the particles is not implemented for spherical coordinates.')
endif
endif
!
! With shear there is an extra term due to the background shear flow.
!
if (lshear .and. lshear_accel_par) &
dfp(1:npar_loc,ivpy) = dfp(1:npar_loc,ivpy) + qshear * Omega * fp(1:npar_loc,ivpx)
endif
!
! Add constant background pressure gradient beta=alpha*H0/r0, where alpha
! comes from a global pressure gradient P = P0*(r/r0)^alpha.
! (the term must be added to the dust equation of motion when measuring
! velocities relative to the shear flow modified by the globalfpressure grad.)
!
if (beta_dPdr_dust/=0.0 .and. t>=tstart_dragforce_par) then
dfp(1:npar_loc,ivpx) = dfp(1:npar_loc,ivpx) + cs20*beta_dPdr_dust_scaled
endif
!
! Gravity on the particles.
! Gravity on particles is implemented only if lparticle_gravity is true which is the default.
!
if (lparticle_gravity) call particle_gravity(f,df,fp,dfp,ineargrid)
!
! The auxiliary has to be set to zero afterwards
!
call diffuse_backreaction(f,df)
!
! Diagnostic output
!
if (ldiagnos) then
if (idiag_nparsum/=0) call sum_name(npar_loc,idiag_nparsum)
if (idiag_nparmin/=0) call max_name(-npar_loc,idiag_nparmin,lneg=.true.)
if (idiag_nparmax/=0) call max_name(+npar_loc,idiag_nparmax)
if (idiag_nparpmax/=0) call max_name(maxval(npar_imn),idiag_nparpmax)
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_rpm/=0) call sum_par_name(sqrt(fp(1:npar_loc,ixp)**2+ &
fp(1:npar_loc,iyp)**2+fp(1:npar_loc,izp)**2),idiag_rpm)
if (idiag_rp2m/=0) call sum_par_name(fp(1:npar_loc,ixp)**2+ &
fp(1:npar_loc,iyp)**2+fp(1:npar_loc,izp)**2,idiag_rp2m)
if (idiag_vpxm/=0) call sum_par_name(fp(1:npar_loc,ivpx),idiag_vpxm)
if (idiag_vpym/=0) call sum_par_name(fp(1:npar_loc,ivpy),idiag_vpym)
if (idiag_vpzm/=0) call sum_par_name(fp(1:npar_loc,ivpz),idiag_vpzm)
if (idiag_vpxvpym/=0) call sum_par_name( &
fp(1:npar_loc,ivpx)*fp(1:npar_loc,ivpy),idiag_vpxvpym)
if (idiag_vpxvpzm/=0) call sum_par_name( &
fp(1:npar_loc,ivpx)*fp(1:npar_loc,ivpz),idiag_vpxvpzm)
if (idiag_vpyvpzm/=0) call sum_par_name( &
fp(1:npar_loc,ivpy)*fp(1:npar_loc,ivpz),idiag_vpyvpzm)
if (idiag_lpxm/=0) call sum_par_name( &
fp(1:npar_loc,iyp)*fp(1:npar_loc,ivpz)- &
fp(1:npar_loc,izp)*fp(1:npar_loc,ivpy),idiag_lpxm)
if (idiag_lpym/=0) call sum_par_name( &
fp(1:npar_loc,izp)*fp(1:npar_loc,ivpx)- &
fp(1:npar_loc,ixp)*fp(1:npar_loc,ivpz),idiag_lpym)
if (idiag_lpzm/=0) call sum_par_name( &
fp(1:npar_loc,ixp)*fp(1:npar_loc,ivpy)- &
fp(1:npar_loc,iyp)*fp(1:npar_loc,ivpx),idiag_lpzm)
if (idiag_lpx2m/=0) call sum_par_name( &
(fp(1:npar_loc,iyp)*fp(1:npar_loc,ivpz)- &
fp(1:npar_loc,izp)*fp(1:npar_loc,ivpy))**2,idiag_lpx2m)
if (idiag_lpy2m/=0) call sum_par_name( &
(fp(1:npar_loc,izp)*fp(1:npar_loc,ivpx)- &
fp(1:npar_loc,ixp)*fp(1:npar_loc,ivpz))**2,idiag_lpy2m)
if (idiag_lpz2m/=0) call sum_par_name( &
(fp(1:npar_loc,ixp)*fp(1:npar_loc,ivpy)- &
fp(1:npar_loc,iyp)*fp(1:npar_loc,ivpx))**2,idiag_lpz2m)
if (idiag_vpx2m/=0) &
call sum_par_name(fp(1:npar_loc,ivpx)**2,idiag_vpx2m)
if (idiag_vpy2m/=0) &
call sum_par_name(fp(1:npar_loc,ivpy)**2,idiag_vpy2m)
if (idiag_vpz2m/=0) &
call sum_par_name(fp(1:npar_loc,ivpz)**2,idiag_vpz2m)
if (idiag_vprms/=0) &
call sum_par_name((fp(1:npar_loc,ivpx)**2 &
+fp(1:npar_loc,ivpy)**2 &
+fp(1:npar_loc,ivpz)**2),idiag_vprms,lsqrt=.true.)
if (idiag_vpyfull2m/=0) &
call sum_par_name((fp(1:npar_loc,ivpy)- &
qshear*Omega*fp(1:npar_loc,ixp))**2,idiag_vpyfull2m)
if (idiag_epotpm/=0) call sum_par_name( &
-gravr/sqrt(sum(fp(1:npar_loc,ixp:izp)**2,dim=2)),idiag_epotpm)
if (idiag_vpmax/=0) call max_par_name( &
sqrt(sum(fp(1:npar_loc,ivpx:ivpz)**2,2)),idiag_vpmax)
if (idiag_vrelpabsm/=0) call calc_relative_velocity(f,fp,ineargrid)
if (idiag_vpxmax/=0) call max_par_name(fp(1:npar_loc,ivpx),idiag_vpxmax)
if (idiag_vpymax/=0) call max_par_name(fp(1:npar_loc,ivpy),idiag_vpymax)
if (idiag_vpzmax/=0) call max_par_name(fp(1:npar_loc,ivpz),idiag_vpzmax)
if (idiag_vpxmin/=0) call max_par_name(-fp(1:npar_loc,ivpx),idiag_vpxmin,lneg=.true.)
if (idiag_vpymin/=0) call max_par_name(-fp(1:npar_loc,ivpy),idiag_vpymin,lneg=.true.)
if (idiag_vpzmin/=0) call max_par_name(-fp(1:npar_loc,ivpz),idiag_vpzmin,lneg=.true.)
if (idiag_eccpxm/=0) call sum_par_name( &
sum(fp(1:npar_loc,ivpx:ivpz)**2,dim=2)*fp(1:npar_loc,ixp)/gravr- &
sum(fp(1:npar_loc,ixp:izp)*fp(1:npar_loc,ivpx:ivpz),dim=2)* &
fp(1:npar_loc,ivpx)/gravr-fp(1:npar_loc,ixp)/ &
sqrt(sum(fp(1:npar_loc,ixp:izp)**2,dim=2)),idiag_eccpxm)
if (idiag_eccpym/=0) call sum_par_name( &
sum(fp(1:npar_loc,ivpx:ivpz)**2,dim=2)*fp(1:npar_loc,iyp)/gravr- &
sum(fp(1:npar_loc,ixp:izp)*fp(1:npar_loc,ivpx:ivpz),dim=2)* &
fp(1:npar_loc,ivpy)/gravr-fp(1:npar_loc,iyp)/ &
sqrt(sum(fp(1:npar_loc,ixp:izp)**2,dim=2)),idiag_eccpym)
if (idiag_eccpzm/=0) call sum_par_name( &
sum(fp(1:npar_loc,ivpx:ivpz)**2,dim=2)*fp(1:npar_loc,izp)/gravr- &
sum(fp(1:npar_loc,ixp:izp)*fp(1:npar_loc,ivpx:ivpz),dim=2)* &
fp(1:npar_loc,ivpz)/gravr-fp(1:npar_loc,izp)/ &
sqrt(sum(fp(1:npar_loc,ixp:izp)**2,dim=2)),idiag_eccpzm)
if (idiag_eccpx2m/=0) call sum_par_name(( &
sum(fp(1:npar_loc,ivpx:ivpz)**2,dim=2)*fp(1:npar_loc,ixp)/gravr- &
sum(fp(1:npar_loc,ixp:izp)*fp(1:npar_loc,ivpx:ivpz),dim=2)* &
fp(1:npar_loc,ivpx)/gravr-fp(1:npar_loc,ixp)/ &
sqrt(sum(fp(1:npar_loc,ixp:izp)**2,dim=2)))**2,idiag_eccpx2m)
if (idiag_eccpy2m/=0) call sum_par_name(( &
sum(fp(1:npar_loc,ivpx:ivpz)**2,dim=2)*fp(1:npar_loc,iyp)/gravr- &
sum(fp(1:npar_loc,ixp:izp)*fp(1:npar_loc,ivpx:ivpz),dim=2)* &
fp(1:npar_loc,ivpy)/gravr-fp(1:npar_loc,iyp)/ &
sqrt(sum(fp(1:npar_loc,ixp:izp)**2,dim=2)))**2,idiag_eccpy2m)
if (idiag_eccpz2m/=0) call sum_par_name(( &
sum(fp(1:npar_loc,ivpx:ivpz)**2,dim=2)*fp(1:npar_loc,izp)/gravr- &
sum(fp(1:npar_loc,ixp:izp)*fp(1:npar_loc,ivpx:ivpz),dim=2)* &
fp(1:npar_loc,ivpz)/gravr-fp(1:npar_loc,izp)/ &
sqrt(sum(fp(1:npar_loc,ixp:izp)**2,dim=2)))**2,idiag_eccpz2m)
if (idiag_npargone/=0) then
call count_particles(ipar,npar_found)
call save_name(float(npar-npar_found),idiag_npargone)
endif
if (idiag_deshearbcsm/=0) then
call sum_name(energy_gain_shear_bcs/npar,idiag_deshearbcsm)
endif
endif
!
if (lfirstcall) lfirstcall=.false.
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(df)
call keep_compiler_quiet(ineargrid)
!
endsubroutine dvvp_dt
!***********************************************************************
subroutine particle_gravity(f,df,fp,dfp,ineargrid)
!
! Contribution of gravity to dvvp_dt
!
! 11-oct-12/dhruba: copied from dvvp_dt
!
use Diagnostics
use SharedVariables, only: get_shared_variable
!
real, dimension (mx,my,mz,mfarray), intent (in) :: f
real, dimension (mx,my,mz,mvar), intent (inout) :: df
real, dimension (mpar_loc,mparray), intent (in) :: fp
real, dimension (mpar_loc,mpvar), intent (inout) :: dfp
integer, dimension (mpar_loc,3), intent (in) :: ineargrid
real, dimension (mpar_loc) :: rpbeta_tmp_arr, OO2_arr, rr_arr, vv_arr
real, dimension(mpar_loc,3) :: gpp_arr
integer, dimension(mpar_loc) :: jspec_arr
!
real, dimension(3) :: ggp
real :: rr=0, vv=0, OO2, rpbeta_tmp=0
integer :: k, jspec
logical :: lheader, lfirstcall=.true.
!
call keep_compiler_quiet(f,df)
!
! Print out header information in first time step.
!
lheader=lfirstcall .and. lroot
if (lheader) then
print*,'dvvp_dt: Calculating gravity'
endif
!
! Gravity in the x-direction.
!
if (t>=tstart_grav_x_par) then
!
select case (gravx_profile)
!
case ('')
if (lheader) print*, 'dvvp_dt: No gravity in x-direction.'
!
case ('zero')
if (lheader) print*, 'dvvp_dt: No gravity in x-direction.'
!
case ('const','plain')
if (lheader) print*, 'dvvp_dt: Constant gravity field in x-direction'
dfp(1:npar_loc,ivpx)=dfp(1:npar_loc,ivpx) + gravx
!
case ('linear')
if (lheader) print*, 'dvvp_dt: Linear gravity field in x-direction.'
dfp(1:npar_loc,ivpx)=dfp(1:npar_loc,ivpx) - &
nu_epicycle2*fp(1:npar_loc,ixp)
!
case ('sinusoidal')
if (lheader) &
print*, 'dvvp_dt: Sinusoidal gravity field in x-direction.'
dfp(1:npar_loc,ivpx)=dfp(1:npar_loc,ivpx) + &
gravx*sin(kx_gg*fp(1:npar_loc,ixp))
!
case default
call fatal_error('dvvp_dt','chosen gravx_profile is not valid!')
!
endselect
!
endif
!
! Gravity in the z-direction.
!
if (t>=tstart_grav_z_par) then
!
select case (gravz_profile)
!
case ('')
if (lheader) print*, 'dvvp_dt: No gravity in z-direction.'
!
case ('zero')
if (lheader) print*, 'dvvp_dt: No gravity in z-direction.'
!
case ('const','plain')
if (lheader) print*, 'dvvp_dt: Constant gravity field in z-direction.'
dfp(1:npar_loc,ivpz)=dfp(1:npar_loc,ivpz) + gravz
!
case ('linear')
if (lheader) print*, 'dvvp_dt: Linear gravity field in z-direction.'
dfp(1:npar_loc,ivpz)=dfp(1:npar_loc,ivpz) - &
nu_epicycle2*fp(1:npar_loc,izp)
!
case ('sinusoidal')
if (lheader) &
print*, 'dvvp_dt: Sinusoidal gravity field in z-direction.'
dfp(1:npar_loc,ivpz)=dfp(1:npar_loc,ivpz) + &
gravz*sin(kz_gg*fp(1:npar_loc,izp))
!
case default
call fatal_error('dvvp_dt','chosen gravz_profile is not valid!')
!
endselect
!
endif
!
! Radial gravity.
!
if (t>=tstart_grav_r_par) then
!
select case (gravr_profile)
!
case ('')
if (lheader) print*, 'dvvp_dt: No radial gravity'
!
case ('zero')
if (lheader) print*, 'dvvp_dt: No radial gravity'
!
case ('newtonian-central','newtonian')
if (lpointmasses) &
call fatal_error('dvvp_dt','You are using massive particles. '//&
'The N-body code should take care of the stellar-like '// &
'gravity on the dust. Switch off the '// &
'gravr_profile=''newtonian'' on particles_init')
if (lheader) &
print*, 'dvvp_dt: Newtonian gravity from a fixed central object'
if (lvector_gravity) then
if (t>=tstart_rpbeta) then
if (lparticles_radius .and. lparticles_radius_rpbeta) then
rpbeta_tmp_arr(1:npar_loc)=fp(1:npar_loc,irpbeta)
elseif (npar_species>1) then
jspec_arr(1:npar_loc)=npar_species*(ipar(1:npar_loc)-1)/npar+1
else
rpbeta_tmp_arr(1:npar_loc)=rpbeta
endif
else
rpbeta_tmp_arr(1:npar_loc)=0.0
endif
if (lcartesian_coords) then
if (lcylindrical_gravity_par) then
rr_arr(1:npar_loc)=sqrt(fp(1:npar_loc,ixp)**2+fp(1:npar_loc,iyp)**2+gravsmooth2)
else
rr_arr(1:npar_loc)=sqrt(fp(1:npar_loc,ixp)**2+fp(1:npar_loc,iyp)**2+fp(1:npar_loc,izp)**2+gravsmooth2)
endif
OO2_arr(1:npar_loc)=rr_arr(1:npar_loc)**(-3.)*gravr*(1.0-rpbeta_tmp_arr(1:npar_loc))
gpp_arr(1:npar_loc,1) = -fp(1:npar_loc,ixp)*OO2_arr(1:npar_loc)
gpp_arr(1:npar_loc,2) = -fp(1:npar_loc,iyp)*OO2_arr(1:npar_loc)
if (lcylindrical_gravity_par) then
gpp_arr(1:npar_loc,3) = 0.
else
gpp_arr(1:npar_loc,3) = -fp(1:npar_loc,izp)*OO2_arr(1:npar_loc)
endif
dfp(1:npar_loc,ivpx:ivpz) = dfp(1:npar_loc,ivpx:ivpz) + gpp_arr(1:npar_loc,:)
elseif (lcylindrical_coords) then
if (lcylindrical_gravity_par) then
rr_arr(1:npar_loc)=sqrt(fp(1:npar_loc,ixp)**2+gravsmooth2)
else
rr_arr(1:npar_loc)=sqrt(fp(1:npar_loc,ixp)**2+fp(1:npar_loc,izp)**2+gravsmooth2)
endif
OO2_arr(1:npar_loc)=rr_arr(1:npar_loc)**(-3.)*gravr*(1.0-rpbeta_tmp_arr(1:npar_loc))
gpp_arr(1:npar_loc,1) = -fp(1:npar_loc,ixp)*OO2_arr(1:npar_loc)
gpp_arr(1:npar_loc,2) = 0.0
if (lcylindrical_gravity_par) then
gpp_arr(1:npar_loc,3) = 0.
else
gpp_arr(1:npar_loc,3) = -fp(1:npar_loc,izp)*OO2_arr(1:npar_loc)
endif
dfp(1:npar_loc,ivpx:ivpz) = dfp(1:npar_loc,ivpx:ivpz) + gpp_arr(1:npar_loc,:)
elseif (lspherical_coords) then
rr_arr(1:npar_loc)=sqrt(fp(1:npar_loc,ixp)**2+gravsmooth2)
OO2_arr(1:npar_loc)=rr_arr(1:npar_loc)**(-3)*gravr*(1.0-rpbeta_tmp_arr(1:npar_loc))
gpp_arr(1:npar_loc,1) = -fp(1:npar_loc,ixp)*OO2_arr(1:npar_loc)
gpp_arr(1:npar_loc,2) = 0.0; gpp_arr(1:npar_loc,3) = 0.0
if (lcylindrical_gravity_par) call fatal_error("dvvp_dt",&
"No cylindrical gravity in spherical coordinates.")
dfp(1:npar_loc,ivpx:ivpz) = dfp(1:npar_loc,ivpx:ivpz) + gpp_arr(1:npar_loc,:)
endif
! Limit time-step if particles close to gravity source.
if (ldt_grav_par.and.(lfirst.and.ldt)) then
if (lcartesian_coords) then
vv_arr(1:npar_loc)=sqrt(fp(1:npar_loc,ivpx)**2+fp(1:npar_loc,ivpy)**2+fp(1:npar_loc,ivpz)**2)
elseif (lcylindrical_coords) then
vv_arr(1:npar_loc)=sqrt(fp(1:npar_loc,ivpx)**2+fp(1:npar_loc,ivpz)**2)
elseif (lspherical_coords) then
vv_arr(1:npar_loc)=abs(fp(1:npar_loc,ivpx))
endif
do k=1,npar_loc
dt1_max(ineargrid(k,1)-nghost) = max(dt1_max(ineargrid(k,1)-nghost),vv_arr(k)/rr_arr(k)/cdtpgrav)
enddo
endif
else
do k=1,npar_loc
if (t>=tstart_rpbeta) then
if (npar_species>1) then
jspec=npar_species*(ipar(k)-1)/npar+1
rpbeta_tmp=rpbeta_species(jspec)
else
rpbeta_tmp=rpbeta
endif
else
rpbeta_tmp=0.0
endif
if (lcartesian_coords) then
if (lcylindrical_gravity_par) then
rr=sqrt(fp(k,ixp)**2+fp(k,iyp)**2+gravsmooth2)
else
rr=sqrt(fp(k,ixp)**2+fp(k,iyp)**2+fp(k,izp)**2+gravsmooth2)
endif
OO2=rr**(-3.)*gravr*(1.0-rpbeta_tmp)
ggp(1) = -fp(k,ixp)*OO2
ggp(2) = -fp(k,iyp)*OO2
if (lcylindrical_gravity_par) then
ggp(3) = 0.
else
ggp(3) = -fp(k,izp)*OO2
endif
dfp(k,ivpx:ivpz) = dfp(k,ivpx:ivpz) + ggp
elseif (lcylindrical_coords) then
if (lcylindrical_gravity_par) then
rr=sqrt(fp(k,ixp)**2+gravsmooth2)
else
rr=sqrt(fp(k,ixp)**2+fp(k,izp)**2+gravsmooth2)
endif
OO2=rr**(-3.)*gravr*(1.0-rpbeta_tmp)
ggp(1) = -fp(k,ixp)*OO2
ggp(2) = 0.0
if (lcylindrical_gravity_par) then
ggp(3) = 0.
else
ggp(3) = -fp(k,izp)*OO2
endif
dfp(k,ivpx:ivpz) = dfp(k,ivpx:ivpz) + ggp
elseif (lspherical_coords) then
rr=sqrt(fp(k,ixp)**2+gravsmooth2)
OO2=rr**(-3)*gravr*(1.0-rpbeta_tmp)
ggp(1) = -fp(k,ixp)*OO2
ggp(2) = 0.0; ggp(3) = 0.0
if (lcylindrical_gravity_par) call fatal_error("dvvp_dt",&
"No cylindrical gravity in spherical coordinates.")
dfp(k,ivpx:ivpz) = dfp(k,ivpx:ivpz) + ggp
endif
! Limit time-step if particles close to gravity source.
if (ldt_grav_par.and.(lfirst.and.ldt)) then
if (lcartesian_coords) then
vv=sqrt(fp(k,ivpx)**2+fp(k,ivpy)**2+fp(k,ivpz)**2)
elseif (lcylindrical_coords) then
vv=sqrt(fp(k,ivpx)**2+fp(k,ivpz)**2)
elseif (lspherical_coords) then
vv=abs(fp(k,ivpx))
endif
dt1_max(ineargrid(k,1)-nghost)= &
max(dt1_max(ineargrid(k,1)-nghost),vv/rr/cdtpgrav)
endif
enddo
endif
!
case default
call fatal_error('dvvp_dt','chosen gravr_profile is not valid!')
!
endselect
!
endif
if (lfirstcall) lfirstcall=.false.
!
endsubroutine particle_gravity
!***********************************************************************
subroutine dxxp_dt_pencil(f,df,fp,dfp,p,ineargrid)
!
! Evolution of particle position (called from main pencil loop).
!
! 25-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
!
integer :: k, ix0, iy0, iz0
real :: dt1_advpx, dt1_advpy, dt1_advpz
!
! Contribution of dust particles to time step.
!
if (lfirst.and.ldt.and.ldt_adv_par) then
if (npar_imn(imn)/=0) then
do k=k1_imn(imn),k2_imn(imn)
ix0=ineargrid(k,1); iy0=ineargrid(k,2); iz0=ineargrid(k,3)
if (nxgrid/=1) then
dt1_advpx=abs(fp(k,ivpx))*dx_1(ix0)
else
dt1_advpx=0.0
endif
if (nygrid/=1) then
if (lshear .and. lcdtp_shear) then
dt1_advpy=(-qshear*Omega*fp(k,ixp)+abs(fp(k,ivpy)))*dy_1(iy0)
else
dt1_advpy=abs(fp(k,ivpy))*dy_1(iy0)
endif
else
dt1_advpy=0.0
endif
if (nzgrid/=1) then
dt1_advpz=abs(fp(k,ivpz))*dz_1(iz0)
else
dt1_advpz=0.0
endif
if (l_shell) then
dt1_advpx=abs(fp(k,ivpx))/k_shell
dt1_advpy=abs(fp(k,ivpy))/k_shell
dt1_advpz=abs(fp(k,ivpz))/k_shell
endif
dt1_max(ix0-nghost)=max(dt1_max(ix0-nghost), &
sqrt(dt1_advpx**2+dt1_advpy**2+dt1_advpz**2)/cdtp)
enddo
endif
endif
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(df)
call keep_compiler_quiet(dfp)
call keep_compiler_quiet(p)
!
endsubroutine dxxp_dt_pencil
!***********************************************************************
subroutine dvvp_dt_pencil(f,df,fp,dfp,p,ineargrid)
!
! Evolution of dust particle velocity (called from main pencil loop).
!
! 25-apr-06/anders: coded
!
use Diagnostics
use Particles_spin, only: calc_liftforce
use Particles_diagnos_dv, only: collisions
use Particles_diagnos_state, only: persistence_check
use SharedVariables, only: get_shared_variable
use Viscosity, only: getnu
!
real, dimension (mx,my,mz,mfarray), intent (inout) :: f
real, dimension (mx,my,mz,mvar), intent (inout) :: df
type (pencil_case) :: p
real, dimension (mpar_loc,mparray), intent (inout) :: fp
real, dimension (mpar_loc,mpvar), intent (inout) :: dfp
integer, dimension (mpar_loc,3), intent (inout) :: ineargrid
!
real, dimension (nx) :: dt1_drag = 0.0, dt1_drag_gas, dt1_drag_dust
real, dimension (nx) :: drag_heat, urel_sum
real, dimension (3) :: dragforce, liftforce, bforce,thermforce, uup, xxp
real, dimension (3) :: adv_der_uup = 0.0
real, dimension(:), allocatable :: rep,stocunn
real :: added_mass_beta = 0.0
real :: rho1_point, tausp1_par, up2, urel
real :: weight, weight_x, weight_y, weight_z
real :: dxp, dyp, dzp, volume_cell
integer :: k, l, ix0, iy0, iz0, ierr, irad
integer :: ixx, iyy, izz, ixx0, iyy0, izz0, ixx1, iyy1, izz1
integer, dimension (3) :: inear
logical :: lsink, lvapour
real, pointer :: pscalar_diff, ap0(:)
real :: Sherwood, mass_trans_coeff, lambda_tilde
real :: dthetadt, mp_vcell
!
real, dimension(k1_imn(imn):k2_imn(imn)) :: nu
character (len=labellen) :: ivis=''
real :: nu_
real :: inversetau
!
! Identify module.
!
if (headtt) then
if (lroot) then
print*,'dvvp_dt_pencil: calculate dvvp_dt'
if (lhydro .and. ldragforce_gas_par) then
print*,'dvvp_dt_pencil: adding feedback from dust to gas'
endif
endif
endif
!
! When the field based handling of passive scalar consumption is enabled, set
! the current pencil of the lncc auxiliary to zero
!
if (ldiffuse_passive) f(:,m,n,idlncc) = 0.0
if (ldiffuse_passive .and. ilncc == 0) &
call fatal_error('particles_lagrangian', &
'ldiffuse_passive needs pscalar_nolog=F')
!
! Initialize the pencils that are calculated within this subroutine
! NILS: Are these initialized somewhere else (with the other pencils)?
!
if (lpenc_requested(i_npvz)) p%npvz=0.
if (lpenc_requested(i_npvz2)) p%npvz2=0.
if (lpenc_requested(i_npuz)) p%npuz=0.
if (lpenc_requested(i_np_rad)) p%np_rad=0.
if (lpenc_requested(i_sherwood)) p%sherwood=0.
!
if (idiag_urel /= 0) urel_sum=0.
!
! Precalculate certain quantities, if necessary.
!
if (npar_imn(imn)/=0) then
!
! Precalculate particle Reynolds numbers.
! NILS: This precalculation, together with the stokes cunningham factor,
! NILS: should be moved to where the interpolation is done.
!
getrep: if (ldraglaw_steadystate .or. lparticles_spin &
.or. ldraglaw_stokesschiller) then
allocate(rep(k1_imn(imn):k2_imn(imn)))
if (.not. allocated(rep)) call fatal_error('dvvp_dt_pencil', &
'unable to allocate sufficient memory for rep', .true.)
call calc_pencil_rep(fp, rep)
endif getrep
!
! Precalculate Stokes-Cunningham factor (only if not ldraglaw_simple or ldraglaw_purestokes)
!
if (.not. (ldraglaw_simple .or. ldraglaw_purestokes &
.or. ldraglaw_stokesschiller)) then
if (ldraglaw_steadystate.or.lbrownian_forces) then
allocate(stocunn(k1_imn(imn):k2_imn(imn)))
if (.not.allocated(stocunn)) then
call fatal_error('dvvp_dt_pencil','unable to allocate sufficient memory for stocunn', .true.)
endif
!
call calc_stokes_cunningham(fp,stocunn)
endif
endif
endif
!
! Drag force on particles and on gas.
!
if (ldragforce_heat .or. (ldiagnos .and. idiag_dedragp/=0)) drag_heat = 0.0
!
if (ldragforce_dust_par .and. t>=tstart_dragforce_par) then
if (headtt) print*,'dvvp_dt: Add drag force; tausp=', tausp
!
if (npar_imn(imn)/=0) then
!
if (lfirst.and.ldt) then
dt1_drag_dust=0.0
if (ldragforce_gas_par) dt1_drag_gas=0.0
endif
!
! Get viscosity used to calculate the pscalar Schmidt number
! NILS: All calculation of the viscosity should be moved to the interpolation
! NILS: part.
!
if (lpscalar .and. lpscalar_sink) then
if (.not. lsherwood_const) then
call getnu(nu_input=nu_,IVIS=ivis)
if (ivis=='nu-const') then
nu=nu_
elseif (ivis=='nu-mixture') then
nu=interp_nu
elseif (ivis=='rho-nu-const') then
nu=nu_/interp_rho(k1_imn(imn):k2_imn(imn))
elseif (ivis=='sqrtrho-nu-const') then
nu=nu_/sqrt(interp_rho(k1_imn(imn):k2_imn(imn)))
elseif (ivis=='nu-therm') then
nu=nu_*sqrt(interp_TT(k1_imn(imn):k2_imn(imn)))
elseif (ivis=='mu-therm') then
nu=nu_*sqrt(interp_TT(k1_imn(imn):k2_imn(imn)))/interp_rho(k1_imn(imn):k2_imn(imn))
else
call fatal_error('dvvp_dt_pencil','No such ivis!')
endif
endif
!
! Get the passive scalar diffusion rate
!
call get_shared_variable('pscalar_diff',pscalar_diff,ierr)
!
endif
!
! Loop over all particles in current pencil.
!
do k=k1_imn(imn),k2_imn(imn)
ix0=ineargrid(k,1)
iy0=ineargrid(k,2)
iz0=ineargrid(k,3)
!
! Calculate required pencils
! NILS: Could this be moved to calc_pencils_particles
if (lpenc_requested(i_npvz) .or. lpenc_requested(i_npvz2) .or. &
lpenc_requested(i_np_rad) .or. lpenc_requested(i_npuz)) then
call get_shared_variable('ap0',ap0,ierr)
do irad=1,npart_radii
if ((fp(k,iap) > ap0(irad)*0.99) .and. (fp(k,iap) < ap0(irad)*1.01)) then
p%npvz(ix0-nghost,irad)=p%npvz(ix0-nghost,irad)+fp(k,ivpz)
p%npvz2(ix0-nghost,irad)=p%npvz2(ix0-nghost,irad)+fp(k,ivpz)**2
p%np_rad(ix0-nghost,irad)=p%np_rad(ix0-nghost,irad)+1.
if (lpenc_requested(i_npuz)) then
!
! DM : I understand that the following is not efficient, it interplolates twice, but seems
! DM: to be the quickest solution now.
!
inear = ineargrid(k,:)
xxp = fp(k,ixp:izp)
if (lparticlemesh_cic) then
call interpolate_linear(f,iux,iuz,xxp,uup,inear,0,ipar(k))
else
if (linterpolate_spline) then
call interpolate_quadratic_spline(f,iux,iuz,xxp,uup,inear,0,ipar(k))
else
call interpolate_quadratic(f,iux,iuz,xxp,uup,inear,0,ipar(k))
endif
endif
p%npuz(ix0-nghost,irad)=p%npuz(ix0-nghost,irad)+uup(3)
endif
endif
enddo
endif
!
! The interpolated gas velocity is either precalculated, and stored in
! interp_uu, or it must be calculated here.
! NILS: In this module uu should always be precalculated - hence, this could be
! NILS: removed.
!
if (.not. interp%luu) then
if (lhydro) then
inear = ineargrid(k,:)
xxp = fp(k,ixp:izp)
if (lparticlemesh_cic) then
call interpolate_linear(f,iux,iuz,xxp,uup,inear,0,ipar(k))
elseif (lparticlemesh_tsc) then
if (linterpolate_spline) then
call interpolate_quadratic_spline(f,iux,iuz,xxp,uup,inear,0,ipar(k))
else
call interpolate_quadratic(f,iux,iuz,xxp,uup,inear,0,ipar(k))
endif
else
uup=f(ix0,iy0,iz0,iux:iuz)
endif
elseif (l_shell) then
call calc_gas_velocity_shell_call(k,uup,fp)
else
uup=0.0
endif
else
uup=interp_uu(k,:)
endif
!
! Track particle state in terms of local gas velocity
!
if (lparticles_diagnos_state .and. lfirst) &
call persistence_check(fp, k, uup)
!
! Get the friction time. For the case of |uup| ~> cs, the Epstein drag law
! is dependent on the relative mach number, hence the need to feed uup as
! an optional argument to get_frictiontime.
!
if (ldraglaw_epstein_transonic .or. ldraglaw_eps_stk_transonic) then
call get_frictiontime(f,fp,p,ineargrid,k,tausp1_par,uup)
elseif (ldraglaw_steadystate) then
call get_frictiontime(f,fp,p,ineargrid,k,tausp1_par,rep=rep(k),stocunn=stocunn(k))
elseif (ldraglaw_stokesschiller) then
call get_frictiontime(f,fp,p,ineargrid,k,tausp1_par,rep=rep(k))
else
call get_frictiontime(f,fp,p,ineargrid,k,tausp1_par)
endif
!
! Calculate and add drag force.
!
dragforce = -tausp1_par*(fp(k,ivpx:ivpz)-uup)
!
! Consider only (spherical) radial component of drag force (for testing).
! NILS: SInce this is for testing only; could it be removed?
!
if (ldragforce_radialonly) then
dragforce=fp(k,ixp:izp)*(dragforce(1)*fp(k,ixp)+ &
dragforce(2)*fp(k,iyp)+dragforce(3)*fp(k,izp))/ &
(fp(k,ixp)**2+fp(k,iyp)**2+fp(k,izp)**2)
endif
!
dfp(k,ivpx:ivpz) = dfp(k,ivpx:ivpz) + dragforce
!
! Account for added mass term beta
! JONAS: The advective derivative of velocity is interpolated for each
! particle
!
if (lbubble) then
if (lhydro) then
inear = ineargrid(k,:)
xxp = fp(k,ixp:izp)
! NILS: SHould do this interpolation as a part of interp!
if (lparticlemesh_cic) then
call interpolate_linear(f,i_adv_derx,i_adv_derz,xxp,adv_der_uup,inear,0,ipar(k))
elseif (lparticlemesh_tsc) then
if (linterpolate_spline) then
call interpolate_quadratic_spline(f,i_adv_derx,i_adv_derz,xxp,adv_der_uup,inear,0,ipar(k))
else
call interpolate_quadratic(f,i_adv_derx,i_adv_derz,xxp,adv_der_uup,inear,0,ipar(k))
endif
else
adv_der_uup=f(ix0,iy0,iz0,i_adv_derx:i_adv_derz)
endif
else
adv_der_uup=0.0
endif
!
! Calculate the beta for the current particle
!
call calc_added_mass_beta(fp,k,added_mass_beta)
!
! Add the contribution of the added mass/virtual mass term to the velocity evolution
!
dfp(k,ivpx:ivpz) = dfp(k,ivpx:ivpz) + added_mass_beta * adv_der_uup
endif
!
! Back-reaction friction force from particles on gas. Three methods are
! implemented for assigning a particle to the mesh (see Hockney & Eastwood):
!
!
! NILS: Must clean up the backreations A LOT. SHould use
! NILS: "find_interpolation_indeces". See e.g. particles_mass for example.
!
!
! 0. NGP (Nearest Grid Point)
! The entire effect of the particle goes to the nearest grid point.
! 1. CIC (Cloud In Cell)
! The particle has a region of influence with the size of a grid cell.
! This is equivalent to a first order (spline) interpolation scheme.
! 2. TSC (Triangular Shaped Cloud)
! The particle is spread over a length of two grid cells, but with
! a density that falls linearly outwards.
! This is equivalent to a second order spline interpolation scheme.
!
if (ldragforce_gas_par .or. (lpscalar_sink .and. lpscalar)) then
!
! Check if the particles consume passive scalar, and calculate the
! consumption rate
!
if (lpscalar_sink .and. lpscalar) then
!
! JONAS: michaelides 2006,p122
! Nu/Sh = 0.922+Pe**0.33+0.1*Pe**0.33*Re**0.33
! Present: rep(k), needed: Pe(k)
! From Multiphase flows with Droplets and particles, p.62:
! Sh = 2+0.69*Re_rel**0.5 * Sc**0.33
! The long number is 0.7**(1/3)
! Sc = nu/pscalar_diff implement with getnu
!
if (lsherwood_const) then
Sherwood = 2.
else
Sherwood = 2.0 + 0.69*sqrt(rep(k))*(nu(k)/pscalar_diff)**(1./3.)
endif
!
if (lpenc_requested(i_sherwood)) then
p%sherwood(ix0-nghost)=p%sherwood(ix0-nghost)+Sherwood
endif
if (idiag_urel /= 0) then
urel=sqrt(sum((interp_uu(k,:) - fp(k,ivpx:ivpz))**2))
urel_sum(ix0-nghost)=urel_sum(ix0-nghost)+urel
endif
!
mass_trans_coeff=Sherwood*pscalar_diff/ &
(2*fp(k,iap))
lambda_tilde=pscalar_sink_rate*mass_trans_coeff/ &
(pscalar_sink_rate+mass_trans_coeff)
dthetadt=lambda_tilde*4.*pi*fp(k,iap)**2
endif
!
! Cloud In Cell (CIC) scheme.
!
if (lparticlemesh_cic) then
ixx0=ix0; iyy0=iy0; izz0=iz0
ixx1=ix0; iyy1=iy0; izz1=iz0
!
! Particle influences param_iothe 8 surrounding grid points. The reference point is
! the grid point at the lower left corner.
!
if ( (x(ix0)>fp(k,ixp)) .and. nxgrid/=1) ixx0=ixx0-1
if ( (y(iy0)>fp(k,iyp)) .and. nygrid/=1) iyy0=iyy0-1
if ( (z(iz0)>fp(k,izp)) .and. nzgrid/=1) izz0=izz0-1
if (nxgrid/=1) ixx1=ixx0+1
if (nygrid/=1) iyy1=iyy0+1
if (nzgrid/=1) izz1=izz0+1
!
do izz=izz0,izz1; do iyy=iyy0,iyy1; do ixx=ixx0,ixx1
weight=1.0
if (nxgrid/=1) weight=weight*( 1.0-abs(fp(k,ixp)-x(ixx))*dx_1(ixx) )
if (nygrid/=1) weight=weight*( 1.0-abs(fp(k,iyp)-y(iyy))*dy_1(iyy) )
if (nzgrid/=1) weight=weight*( 1.0-abs(fp(k,izp)-z(izz))*dz_1(izz) )
! Save the calculation of rho1 when inside pencil.
if ( (iyy/=m).or.(izz/=n).or.(ixx<l1).or.(ixx>l2) ) then
rho1_point = 1.0 / get_gas_density(f,ixx,iyy,izz)
else
rho1_point = p%rho1(ixx-nghost)
endif
!
! Add friction force to grid point.
!NILS: The grid volume should be put into a pencil when required
! if ((lpscalar_sink .and. lpscalar) .or. &
! (ldragforce_gas_par .and. ldraglaw_steadystate)) &
!DM : the above if condition is always true.
!DM : the following call is being made twice.
! call find_grid_volume(ixx,iyy,izz,volume_cell)
! alexrichert: above call to find_grid_volume is superfluous, not sure why
! conditions are different from call below. Perhaps lparticles_radius or
! iap>0 would be a better condition than eps_dtog/ldraglaw_steadystate?
!
if (lhydro .and. ldragforce_gas_par) then
call find_grid_volume(ixx,iyy,izz,volume_cell)
mp_vcell=4.*pi*fp(k,iap)**3*rhopmat/(3.*volume_cell)
if (.not. ldiffuse_dragf) then
df(ixx,iyy,izz,iux:iuz)=df(ixx,iyy,izz,iux:iuz) - &
mp_vcell*rho1_point*dragforce*weight
endif
endif
!
if (lpscalar_sink .and. lpscalar) then
if (ilncc == 0) then
call fatal_error('dvvp_dt_pencil','lpscalar_sink not allowed for pscalar_nolog!')
else
if (.not. ldiffuse_passive) then
call find_grid_volume(ixx,iyy,izz,volume_cell)
df(ixx,iyy,izz,ilncc) = df(ixx,iyy,izz,ilncc) - weight*dthetadt/volume_cell
endif
endif
endif
enddo; enddo; enddo
!
! Triangular Shaped Cloud (TSC) scheme.
!
elseif (lparticlemesh_tsc) then
if (.not. lparticlemesh_pqs_assignment) then
!
! Particle influences the 27 surrounding grid points, but has a density that
! decreases with the distance from the particle centre.
!
if (nxgrid/=1) then
ixx0=ix0-1; ixx1=ix0+1
else
ixx0=ix0 ; ixx1=ix0
endif
if (nygrid/=1) then
iyy0=iy0-1; iyy1=iy0+1
else
iyy0=iy0 ; iyy1=iy0
endif
if (nzgrid/=1) then
izz0=iz0-1; izz1=iz0+1
else
izz0=iz0 ; izz1=iz0
endif
!
! The nearest grid point is influenced differently than the left and right
! neighbours are. A particle that is situated exactly on a grid point gives
! 3/4 contribution to that grid point and 1/8 to each of the neighbours.
!
do izz=izz0,izz1; do iyy=iyy0,iyy1; do ixx=ixx0,ixx1
if ( ((ixx-ix0)==-1) .or. ((ixx-ix0)==+1) ) then
weight_x=1.125-1.5* abs(fp(k,ixp)-x(ixx))*dx_1(ixx) + &
0.5*(abs(fp(k,ixp)-x(ixx))*dx_1(ixx))**2
else
if (nxgrid/=1) &
weight_x=0.75-((fp(k,ixp)-x(ixx))*dx_1(ixx))**2
endif
if ( ((iyy-iy0)==-1) .or. ((iyy-iy0)==+1) ) then
weight_y=1.125-1.5* abs(fp(k,iyp)-y(iyy))*dy_1(iyy) + &
0.5*(abs(fp(k,iyp)-y(iyy))*dy_1(iyy))**2
else
if (nygrid/=1) &
weight_y=0.75-((fp(k,iyp)-y(iyy))*dy_1(iyy))**2
endif
if ( ((izz-iz0)==-1) .or. ((izz-iz0)==+1) ) then
weight_z=1.125-1.5* abs(fp(k,izp)-z(izz))*dz_1(izz) + &
0.5*(abs(fp(k,izp)-z(izz))*dz_1(izz))**2
else
if (nzgrid/=1) &
weight_z=0.75-((fp(k,izp)-z(izz))*dz_1(izz))**2
endif
!
weight=1.0
!
if (nxgrid/=1) weight=weight*weight_x
if (nygrid/=1) weight=weight*weight_y
if (nzgrid/=1) weight=weight*weight_z
! Save the calculation of rho1 when inside pencil.
if ( (iyy/=m).or.(izz/=n).or.(ixx<l1).or.(ixx>l2) ) then
rho1_point = 1.0 / get_gas_density(f,ixx,iyy,izz)
else
rho1_point = p%rho1(ixx-nghost)
endif
!NILS: The grid volume should be put into a pencil when required
if ((lpscalar_sink .and. lpscalar) .or. &
(ldragforce_gas_par .and. ldraglaw_steadystate)) &
call find_grid_volume(ixx,iyy,izz,volume_cell)
! Add friction force to grid point.
if (lhydro .and. ldragforce_gas_par) then
! Calculate the particle mass divided by the cell volume
call find_grid_volume(ixx,iyy,izz,volume_cell)
mp_vcell=4.*pi*fp(k,iap)**3*rhopmat/(3.*volume_cell)
if (.not.lcompensate_sedimentation) then
if (.not. ldiffuse_dragf) then
df(ixx,iyy,izz,iux:iuz)=df(ixx,iyy,izz,iux:iuz) - &
mp_vcell*rho1_point*dragforce*weight
endif
else
df(ixx,iyy,izz,iux:iuz)=df(ixx,iyy,izz,iux:iuz) - &
mp_vcell*rho1_point*dragforce*weight*compensate_sedimentation
endif
endif
if (lpscalar_sink .and. lpscalar) then
if (ilncc == 0) then
call fatal_error('dvvp_dt_pencil',&
'lpscalar_sink not allowed for pscalar_nolog!')
else
if (.not. ldiffuse_passive) then
call find_grid_volume(ixx,iyy,izz,volume_cell)
df(ixx,iyy,izz,ilncc) = df(ixx,iyy,izz,ilncc) - &
weight*dthetadt/volume_cell
endif
endif
endif
enddo; enddo; enddo
else
!
! QPS assignment, one order higher than TSC. Experimental.
!
if (nxgrid/=1) then
ixx0=ix0-2; ixx1=ix0+2
else
ixx0=ix0 ; ixx1=ix0
endif
if (nygrid/=1) then
iyy0=iy0-2; iyy1=iy0+2
else
iyy0=iy0 ; iyy1=iy0
endif
if (nzgrid/=1) then
izz0=iz0-2; izz1=iz0+2
else
izz0=iz0 ; izz1=iz0
endif
do izz=izz0,izz1; do iyy=iyy0,iyy1; do ixx=ixx0,ixx1
dxp=(x(ixx)-fp(k,ixp))*dx_1(ixx)
dyp=(y(iyy)-fp(k,iyp))*dy_1(iyy)
dzp=(z(izz)-fp(k,izp))*dz_1(izz)
!
if (abs(dxp)<=1.0) then
weight_x=(2/3.0)+dxp**2*(abs(dxp)/2-1.0)
elseif (abs(dxp)<=2.0) then
weight_x=(1/6.0)*(2.0-abs(dxp))**3
else
weight_x=0.0
endif
!
if (abs(dyp)<=1.0) then
weight_y=(2/3.0)+dyp**2*(abs(dyp)/2-1.0)
elseif (abs(dyp)<=2.0) then
weight_y=(1/6.0)*(2.0-abs(dyp))**3
else
weight_y=0.0
endif
!
if (abs(dzp)<=1.0) then
weight_z=(2/3.0)+dzp**2*(abs(dzp)/2-1.0)
elseif (abs(dzp)<=2.0) then
weight_z=(1/6.0)*(2.0-abs(dzp))**3
else
weight_z=0.0
endif
!
weight=1.0
if (nxgrid/=1) weight=weight*weight_x
if (nygrid/=1) weight=weight*weight_y
if (nzgrid/=1) weight=weight*weight_z
!
if ( (iyy/=m).or.(izz/=n).or.(ixx<l1).or.(ixx>l2) ) then
rho1_point = 1.0 / get_gas_density(f,ixx,iyy,izz)
else
rho1_point = p%rho1(ixx-nghost)
endif
!
!NILS: The grid volume should be put into a pencil when required
if (ldragforce_gas_par .and. ldraglaw_steadystate) then
call find_grid_volume(ixx,iyy,izz,volume_cell)
endif
if (lhydro .and. ldragforce_gas_par) then
! Calculate the particle mass divided by the cell volume
call find_grid_volume(ixx,iyy,izz,volume_cell)
mp_vcell=4.*pi*fp(k,iap)**3*rhopmat/(3.*volume_cell)
if (.not.lcompensate_sedimentation) then
df(ixx,iyy,izz,iux:iuz)=df(ixx,iyy,izz,iux:iuz) - &
mp_vcell*rho1_point*dragforce*weight
else
df(ixx,iyy,izz,iux:iuz)=df(ixx,iyy,izz,iux:iuz) - &
mp_vcell*rho1_point*dragforce*weight*compensate_sedimentation
endif
endif
if (lpscalar_sink .and. lpscalar) then
if (ilncc == 0) then
call fatal_error('dvvp_dt_pencil',&
'lpscalar_sink not allowed for pscalar_nolog!')
else
if (.not. ldiffuse_passive) then
call find_grid_volume(ixx,iyy,izz,volume_cell)
df(ixx,iyy,izz,ilncc) = df(ixx,iyy,izz,ilncc) - &
weight*dthetadt/volume_cell
endif
endif
endif
enddo; enddo; enddo
endif
!
! GAussian Box (GAB) scheme.
!
elseif (lparticlemesh_gab) then
!
!
!
call find_grid_volume(ix0,iy0,iz0,volume_cell)
!
! In case the diffusion of particle-fluid interaction is activated, the values are only taken from
! the nearest grid point
!
if (ldiffuse_dragf .and. ldiffuse_passive) then
if ( (iy0/=m).or.(iz0/=n).or.(ix0<l1).or.(ix0>l2) ) then
rho1_point = 1.0 / get_gas_density(f,ix0,iy0,iz0)
else
rho1_point = p%rho1(ix0-nghost)
endif
if (lhydro .and. ldragforce_gas_par) then
if (ldraglaw_steadystate) then
mp_vcell=4.*pi*fp(k,iap)**3*rhopmat/(3.*volume_cell)
endif
else
call fatal_error('particles_lagrangian','diffusion of dragforce needs &
& ldragforce_gas_par and ldraglaw_steadystate')
endif
else
!
ixx0=ix0-3
ixx1=ix0+3
iyy0=iy0-3
iyy1=iy0+3
izz0=iz0-3
izz1=iz0+3
!
do izz=izz0,izz1; do iyy=iyy0,iyy1; do ixx=ixx0,ixx1
weight = 1.
if (nxgrid/=1) weight=weight*gab_weights(abs(ixx-ix0)+1)
if (nygrid/=1) weight=weight*gab_weights(abs(iyy-iy0)+1)
if (nzgrid/=1) weight=weight*gab_weights(abs(izz-iz0)+1)
!
if ( (iyy/=m).or.(izz/=n).or.(ixx<l1).or.(ixx>l2) ) then
rho1_point = 1.0 / get_gas_density(f,ixx,iyy,izz)
else
rho1_point = p%rho1(ixx-nghost)
endif
! Add friction force to grid point.
! alexrichert: above call to find_grid_volume is superfluous, not sure why
! conditions are different from call below. Perhaps lparticles_radius or
! iap>0 would be a better condition than eps_dtog/ldraglaw_steadystate?
if (lhydro .and. ldragforce_gas_par) then
mp_vcell=4.*pi*fp(k,iap)**3*rhopmat/(3.*volume_cell)
if (.not. ldiffuse_dragf) then
df(ixx,iyy,izz,iux:iuz)=df(ixx,iyy,izz,iux:iuz) - &
mp_vcell*rho1_point*dragforce*weight
endif
endif
if (lpscalar_sink .and. lpscalar) then
if (ilncc == 0) then
call fatal_error('dvvp_dt_pencil','lpscalar_sink not allowed for pscalar_nolog!')
else
if (.not. ldiffuse_passive) then
df(ixx,iyy,izz,ilncc) = df(ixx,iyy,izz,ilncc) - &
weight*dthetadt/volume_cell
endif
endif
endif
enddo
enddo
enddo
endif
else
!
! Nearest Grid Point (NGP) scheme.
!
l=ineargrid(k,1)
!NILS: The grid volume should be put into a pencil when required
if ((lpscalar_sink .and. lpscalar) .or. &
(ldragforce_gas_par .and. ldraglaw_steadystate)) &
call find_grid_volume(ix0,iy0,iz0,volume_cell)
if (lhydro .and. ldragforce_gas_par) then
! Calculate the particle mass divided by the cell volume
call find_grid_volume(ix0,iy0,iz0,volume_cell)
mp_vcell=4.*pi*fp(k,iap)**3*rhopmat/(3.*volume_cell)
if (.not.lcompensate_sedimentation) then
df(l,m,n,iux:iuz) = df(l,m,n,iux:iuz) - &
mp_vcell*p%rho1(l-nghost)*dragforce
else
df(l,m,n,iux:iuz) = df(l,m,n,iux:iuz) - &
mp_vcell*p%rho1(l-nghost)*dragforce*compensate_sedimentation
endif
endif
if (lpscalar_sink .and. lpscalar) then
if (ilncc == 0) then
call fatal_error('dvvp_dt_pencil',&
'lpscalar_sink not allowed for pscalar_nolog!')
else
if (.not. ldiffuse_passive) then
call find_grid_volume(l,m,n,volume_cell)
df(l,m,n,ilncc) = df(l,m,n,ilncc) - dthetadt/volume_cell
endif
endif
endif
endif
!
! Adding the passive scalar contribution to the auxiliary grid
!
if (ldiffuse_passive .and. lpscalar_sink .and. &
.not. lpencil_check_at_work) then
call find_grid_volume(ix0,iy0,iz0,volume_cell)
f(ix0,iy0,iz0,idlncc) = f(ix0,iy0,iz0,idlncc) - dthetadt/volume_cell
endif
!
! Adding the particle dragforce contribution to the auxiliary grid
!
if (ldiffuse_dragf .and. ldragforce_gas_par .and. &
.not. lpencil_check_at_work) then
call find_grid_volume(ix0,iy0,iz0,volume_cell)
f(ix0,iy0,iz0,idfx:idfz) = f(ix0,iy0,iz0,idfx:idfz) &
-mp_vcell*rho1_point*dragforce
endif
!
endif
!
! Calculate particle mass density in grid cell
!
call find_grid_volume(ix0,iy0,iz0,volume_cell)
mp_vcell=4.*pi*fp(k,iap)**3*rhopmat/(3.*volume_cell)
!
! Heating of gas due to drag force.
!
if (ldragforce_heat .or. (ldiagnos .and. idiag_dedragp/=0)) then
if (ldragforce_gas_par) then
up2=sum((fp(k,ivpx:ivpz)-uup)**2)
else
up2=sum(fp(k,ivpx:ivpz)*(fp(k,ivpx:ivpz)-uup))
endif
!
! WL: Check if this is right. drag_heat is of dimension nx, not a scalar
!
drag_heat(ix0-nghost)=drag_heat(ix0-nghost) + &
mp_vcell*tausp1_par*up2
endif
!
! The minimum friction time of particles in a grid cell sets the local friction
! time-step when there is only drag force on the dust,
! dt1_drag = max(1/tausp)
!
! With drag force on the gas as well, the maximum time-step is set as
! dt1_drag = Sum_k[eps_k/tau_k]
!
getdt1: if (lfirst .and. ldt) then
dt1_drag_dust(ix0-nghost) = max(dt1_drag_dust(ix0-nghost), tausp1_par)
if (ldragforce_gas_par) &
dt1_drag_gas(ix0-nghost) = dt1_drag_gas(ix0-nghost) + mp_vcell * p%rho1(ix0-nghost) * tausp1_par
endif getdt1
enddo
!
! Add drag force heating in pencils.
!
if (lentropy .and. ldragforce_heat) &
df(l1:l2,m,n,iss) = df(l1:l2,m,n,iss) + p%rho1*p%TT1*drag_heat
!
! Contribution of friction force to time-step. Dust and gas inverse friction
! time-steps are added up to give a valid expression even when the two are
! of similar magnitude.
!
if (lfirst.and.ldt) then
if (ldragforce_gas_par) then
dt1_drag=dt1_drag_dust+dt1_drag_gas
else
dt1_drag=dt1_drag_dust
endif
dt1_drag=dt1_drag/cdtp_drag
dt1_max=max(dt1_max,dt1_drag)
endif
else
!
! No particles in this pencil.
!
if (lfirst.and.ldt) dt1_drag=0.0
endif
endif
!
! Add friction force from gas that moves systematically slower. Can be used
! to mimic e.g. a sub-Keplerian gas flow without using the Hydro module.
!
if (Deltauy_gas_friction/=0.0 .and. t>=tstart_dragforce_par) then
if (npar_imn(imn)/=0) then
do k=k1_imn(imn),k2_imn(imn)
call get_frictiontime(f,fp,p,ineargrid,k,tausp1_par)
dfp(k,ivpy) = dfp(k,ivpy) - Deltauy_gas_friction*tausp1_par
enddo
endif
endif
!
! Collisional cooling is in a separate subroutine.
!
if ( (lcollisional_cooling_taucool .or. lcollisional_cooling_rms .or. &
lcollisional_cooling_twobody .or. lcollisional_dragforce_cooling) &
.and. t>=tstart_collisional_cooling) &
call collisional_cooling(f,df,fp,dfp,p,ineargrid)
!
! Compensate for increased friction time by appying extra friction force to
! particles.
!
if (lcompensate_friction_increase) &
call compensate_friction_increase(f,fp,dfp,p,ineargrid)
!
! Add lift forces.
!
if (lparticles_spin .and. t>=tstart_liftforce_par) then
if (npar_imn(imn)/=0) then
do k=k1_imn(imn),k2_imn(imn)
call calc_liftforce(fp(k,:), k, rep(k), liftforce)
dfp(k,ivpx:ivpz)=dfp(k,ivpx:ivpz)+liftforce
enddo
endif
endif
!
! Add Brownian forces.
!
if (lbrownian_forces .and. t>=tstart_brownian_par) then
if (npar_imn(imn)/=0) then
do k=k1_imn(imn),k2_imn(imn)
call calc_brownian_force(fp,k,ineargrid(k,:),stocunn(k),bforce)
dfp(k,ivpx:ivpz)=dfp(k,ivpx:ivpz)+bforce
enddo
endif
endif
!
! Add Thermophoretic forces.
!
if (lthermophoretic_forces) then
if (npar_imn(imn)/=0) then
do k=k1_imn(imn),k2_imn(imn)
call calc_thermophoretic_force(fp,k,ineargrid(k,:),thermforce)
dfp(k,ivpx:ivpz)=dfp(k,ivpx:ivpz)+thermforce
enddo
endif
endif
!
! For stiff drag force equations we need to store the forces that are
! unique to the gas.
!
if (ldragforce_stiff .and. .not. lpencil_check_at_work) then
f(l1:l2,m,n,ifgx:ifgz)=p%fpres+p%jxbr+p%fvisc
endif
!
!
! Diagnostic output.
!
if (ldiagnos) then
if (idiag_npm/=0) call sum_mn_name(p%np,idiag_npm)
if (idiag_np2m/=0) call sum_mn_name(p%np**2,idiag_np2m)
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_rhopm/=0) call sum_mn_name(p%rhop,idiag_rhopm)
if (idiag_rhop2m/=0 ) call sum_mn_name(p%rhop**2,idiag_rhop2m)
if (idiag_rhoprms/=0) call sum_mn_name(p%rhop**2,idiag_rhoprms,lsqrt=.true.)
if (idiag_rhopmax/=0) call max_mn_name(p%rhop,idiag_rhopmax)
if (idiag_rhopmin/=0) call max_mn_name(-p%rhop,idiag_rhopmin,lneg=.true.)
if (idiag_epspmax/=0) call max_mn_name(p%epsp,idiag_epspmax)
if (idiag_epspmin/=0) call max_mn_name(-p%epsp,idiag_epspmin,lneg=.true.)
if (idiag_epspm/=0) call sum_mn_name(p%epsp,idiag_epspm)
if (idiag_dedragp/=0) call sum_mn_name(drag_heat,idiag_dedragp)
if (idiag_Shm/=0) call sum_mn_name(p%sherwood/npar*nwgrid,idiag_Shm)
if (idiag_urel/=0) call sum_mn_name(urel_sum/npar*nwgrid,idiag_urel)
if (idiag_dvpx2m/=0 .or. idiag_dvpx2m/=0 .or. idiag_dvpx2m/=0 .or. &
idiag_dvpm /=0 .or. idiag_dvpmax/=0) &
call calculate_rms_speed(fp,ineargrid,p)
if (idiag_dtdragp/=0.and.(lfirst.and.ldt)) &
call max_mn_name(dt1_drag,idiag_dtdragp,l_dt=.true.)
endif
!
! 1d-averages. Happens at every it1d timesteps, NOT at every it1
!
if (l1davgfirst) then
if (idiag_npmx/=0) call yzsum_mn_name_x(p%np,idiag_npmx)
if (idiag_npmy/=0) call xzsum_mn_name_y(p%np,idiag_npmy)
if (idiag_npmz/=0) call xysum_mn_name_z(p%np,idiag_npmz)
if (idiag_rhopmx/=0) call yzsum_mn_name_x(p%rhop,idiag_rhopmx)
if (idiag_rhopmy/=0) call xzsum_mn_name_y(p%rhop,idiag_rhopmy)
if (idiag_rhopmz/=0) call xysum_mn_name_z(p%rhop,idiag_rhopmz)
if (idiag_rhop2mx/=0) call yzsum_mn_name_x(p%rhop**2,idiag_rhop2mx)
if (idiag_rhop2my/=0) call xzsum_mn_name_y(p%rhop**2,idiag_rhop2my)
if (idiag_rhop2mz/=0) call xysum_mn_name_z(p%rhop**2,idiag_rhop2mz)
if (idiag_epspmx/=0) call yzsum_mn_name_x(p%epsp,idiag_epspmx)
if (idiag_epspmy/=0) call xzsum_mn_name_y(p%epsp,idiag_epspmy)
if (idiag_epspmz/=0) call xysum_mn_name_z(p%epsp,idiag_epspmz)
if (idiag_rhopmr/=0) call phizsum_mn_name_r(p%rhop,idiag_rhopmr)
!
do k=1,ndustrad
if (idiag_npvzmz(k)/=0) call xysum_mn_name_z(p%npvz(:,k),idiag_npvzmz(k))
if (idiag_npvz2mz(k)/=0) call xysum_mn_name_z(p%npvz2(:,k),idiag_npvz2mz(k))
if (idiag_npuzmz(k)/=0) call xysum_mn_name_z(p%npuz(:,k),idiag_npuzmz(k))
if (idiag_nptz(k)/=0) call xysum_mn_name_z(p%np_rad(:,k),idiag_nptz(k))
enddo
endif
!
if (l2davgfirst) then
if (idiag_npmxy/=0) call zsum_mn_name_xy(p%np,idiag_npmxy)
if (idiag_rhopmphi/=0) call phisum_mn_name_rz(p%rhop,idiag_rhopmphi)
if (idiag_rhopmxy/=0) call zsum_mn_name_xy(p%rhop,idiag_rhopmxy)
if (idiag_rhopmxz/=0) call ysum_mn_name_xz(p%rhop,idiag_rhopmxz)
if (idiag_sigmap /= 0) call zsum_mn_name_xy(p%rhop, idiag_sigmap, lint=.true.)
endif
!
! particle-particle separation and relative velocity diagnostics
!
if (lparticles_diagnos_dv .and. lfirstpoint .and. lfirst) then
if (t > t_nextcol) call collisions(fp)
endif
!
! Clean up (free allocated memory).
!
if (allocated(rep)) deallocate(rep)
if (allocated(stocunn)) deallocate(stocunn)
!
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)
!
! Subroutine for taking particles out of the simulation due to their proximity
! to a sink particle or sink point.
!
! 25-sep-08/anders: coded
!
use Mpicomm
use Solid_Cells
!
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
!
real, dimension(3) :: momp_swarm_removed, momp_swarm_removed_send
real :: rp, rp_box, rhop_swarm_removed, rhop_swarm_removed_send
real :: xsinkpar, ysinkpar, zsinkpar
integer :: k, ksink, iproc_sink, iproc_sink_send
integer :: ix, ix1, ix2, iy, iy1, iy2, iz, iz1, iz2
integer, parameter :: itag1=100, itag2=101
real :: apar
!
! Sinkparticle activated at time tstart_sink_par
if (t <= tstart_sink_par) return
!
if (lsinkpoint) then
k=1
do while (k<=npar_loc)
rp=sqrt((fp(k,ixp)-xsinkpoint)**2+(fp(k,iyp)-ysinkpoint)**2+ &
(fp(k,izp)-zsinkpoint)**2)
if (rp<rsinkpoint) then
call remove_particle(fp,ipar,k,dfp,ineargrid)
else
k=k+1
endif
enddo
endif
!
if (lsinkparticle_1) then
!
! Search for position of particle with index 1.
!
xsinkpoint=0.0
ysinkpoint=0.0
zsinkpoint=0.0
iproc_sink=-1
do k=1,npar_loc
if (ipar(k)==1) then
ksink=k
iproc_sink=iproc
xsinkpar=fp(k,ixp)
ysinkpar=fp(k,iyp)
zsinkpar=fp(k,izp)
endif
enddo
!
! Communicate position of sink particle to all processors.
!
iproc_sink_send=iproc_sink
call mpireduce_max_int(iproc_sink_send,iproc_sink)
call mpibcast_int(iproc_sink)
call mpibcast_real(xsinkpar,proc=iproc_sink)
call mpibcast_real(ysinkpar,proc=iproc_sink)
call mpibcast_real(zsinkpar,proc=iproc_sink)
!
! Remove particles that are too close to sink particle.
!
if (lparticles_density) then
rhop_swarm_removed=0.0
momp_swarm_removed=0.0
endif
!
k=1
do while (k<=npar_loc)
rp=-1.0
!
! Take into account periodic boundary conditions.
!
if (bcpx=='p' .and. nxgrid/=1) then
ix1=-1; ix2=+1
else
ix1=0; ix2=0
endif
if (bcpy=='p' .and. nygrid/=1) then
iy1=-1; iy2=+1
else
iy1=0; iy2=0
endif
if (bcpz=='p' .and. nzgrid/=1) then
iz1=-1; iz2=+1
else
iz1=0; iz2=0
endif
do iz=iz1,iz2; do iy=iy1,iy2; do ix=ix1,ix2
if (lshear) then
rp_box=sqrt((fp(k,ixp)+Lx*ix-xsinkpar)**2+ &
(fp(k,iyp)+Ly*iy-deltay*ix-ysinkpar)**2+ &
(fp(k,izp)+Lz*iz-zsinkpar)**2)
else
rp_box=sqrt((fp(k,ixp)+Lx*ix-xsinkpar)**2+ &
(fp(k,iyp)+Ly*iy-ysinkpar)**2+ &
(fp(k,izp)+Lz*iz-zsinkpar)**2)
endif
if (rp<0.0) then
rp=rp_box
else
rp=min(rp,rp_box)
endif
enddo; enddo; enddo
if (ipar(k)/=1 .and. rp<rsinkparticle_1) then
if (lparticles_density) then
rhop_swarm_removed = rhop_swarm_removed + fp(k,irhopswarm)
momp_swarm_removed = momp_swarm_removed + &
fp(k,ivpx:ivpz)*fp(k,irhopswarm)
! if (lshear) then
! momp_swarm_removed(2) = momp_swarm_removed(2) + &
! Sshear*fp(k,ixp)*fp(k,irhopswarm)
! endif
endif
call remove_particle(fp,ipar,k,dfp,ineargrid)
else
k=k+1
endif
enddo
!
! Add mass and momentum to sink particle.
!
if (lparticles_density) then
rhop_swarm_removed_send=rhop_swarm_removed
momp_swarm_removed_send=momp_swarm_removed
call mpireduce_sum(rhop_swarm_removed_send,rhop_swarm_removed)
call mpireduce_sum(momp_swarm_removed_send,momp_swarm_removed,3)
if (lroot) then
call mpisend_real(rhop_swarm_removed,iproc_sink,itag1)
call mpisend_real(momp_swarm_removed,3,iproc_sink,itag2)
endif
if (iproc==iproc_sink) then
call mpirecv_real(rhop_swarm_removed,0,itag1)
call mpirecv_real(momp_swarm_removed,3,0,itag2)
!
! Need to find sink particle again since particle removal may have shifted
! its index.
!
do k=1,npar_loc
if (ipar(k)==1) then
ksink=k
iproc_sink=iproc
xsinkpar=fp(k,ixp)
ysinkpar=fp(k,iyp)
zsinkpar=fp(k,izp)
endif
enddo
! if (lshear) then
! fp(ksink,ivpx) = (fp(ksink,ivpx)* &
! fp(ksink,irhopswarm) + momp_swarm_removed(1))/ &
! (fp(ksink,irhopswarm) + rhop_swarm_removed)
! fp(ksink,ivpy) = ((fp(ksink,ivpy)+Sshear*fp(ksink,ixp))* &
! fp(ksink,irhopswarm) + momp_swarm_removed(2))/ &
! (fp(ksink,irhopswarm) + rhop_swarm_removed) - &
! Sshear*fp(ksink,ixp)
! fp(ksink,ivpz) = (fp(ksink,ivpz)* &
! fp(ksink,irhopswarm) + momp_swarm_removed(3))/ &
! (fp(ksink,irhopswarm) + rhop_swarm_removed)
! else
fp(ksink,ivpx:ivpz) = (fp(ksink,ivpx:ivpz)* &
fp(ksink,irhopswarm) + momp_swarm_removed)/ &
(fp(ksink,irhopswarm) + rhop_swarm_removed)
! endif
fp(ksink,irhopswarm) = fp(ksink,irhopswarm) + rhop_swarm_removed
endif
endif
endif
!
! Remove particles if they are within a solid geometry
!
if (lsolid_cells) then
k=1
do while (k<=npar_loc)
if (lparticles_radius) then
apar = fp(k,iap)
else
apar = particle_radius
endif
if (in_solid_cell(fp(k,ixp:izp),apar)) then
call remove_particle(fp,ipar,k,dfp,ineargrid)
else
k=k+1
endif
enddo
endif
!
call keep_compiler_quiet(f)
!
endsubroutine remove_particles_sink_simple
!***********************************************************************
subroutine create_particles_sink_simple(f,fp,dfp,ineargrid)
!
! Subroutine for creating new sink particles or sink points.
!
! Just a dummy routine for now.
!
! 25-sep-08/anders: coded
!
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 get_frictiontime(f,fp,p,ineargrid,k,tausp1_par,uup,&
nochange_opt,rep,stocunn)
!
! Calculate the friction time.
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mpar_loc,mparray) :: fp
type (pencil_case) :: p
integer, dimension (mpar_loc,3) :: ineargrid
integer :: k
real :: tausp1_par
real, optional, dimension(3), intent(in) :: uup
logical, optional :: nochange_opt
real, optional, intent(in) :: rep
real, optional :: stocunn
!
real :: tausg1_point,OO, tmp
integer :: ix0, iy0, iz0, inx0, jspec
logical :: nochange=.true.
!
if (present(nochange_opt)) then
if (nochange_opt) then
nochange=.true.
else
nochange=.false.
endif
else
nochange=.false.
endif
!
ix0=ineargrid(k,1); iy0=ineargrid(k,2); iz0=ineargrid(k,3)
inx0=ix0-nghost
!
! Epstein drag law.
!
if (ldraglaw_epstein) then
if (iap/=0) then
if (fp(k,iap)/=0.0) then
tausp1_par = (sqrt(p%cs2(inx0))*p%rho(inx0))/(fp(k,iap)*rhopmat)
!
! DM : 10 Nov 2016
! For the usual Epstein drag we need also the thermal velocity and the
! density at this point. None of them fluctuate by large amounts, so it is better
! to get them from the same pencil without any interpolation. This may be even
! a better prescription for flows with shocks if the particle sits very close to
! a shock.
!
endif
else
! Check if we are using multiple or single particle species.
if (npar_species>1) then
jspec=npar_species*(ipar(k)-1)/npar+1
tmp=tausp1_species(jspec)
else
tmp=tausp1
endif
!
! Scale friction time with local density.
!
if (ldraglaw_variable_density) then
tausp1_par = tmp * get_gas_density(f,ix0,iy0,iz0)
!
! Discriminate between constant tau and special case for
! 1/tau=omega when omega is not constant (as, for instance,
! global Keplerian disks, for which omega=rad**(-3/2)
!
elseif (ldraglaw_variable) then
if (lcartesian_coords) then
OO=(fp(k,ixp)**2 + fp(k,iyp)**2)**(-0.75)
elseif (lcylindrical_coords) then
OO=fp(k,ixp)**(-1.5)
elseif (lspherical_coords) then
OO=(fp(k,ixp)*sin(fp(k,iyp)))**(-1.5)
else
call fatal_error("get_frictiontime", "no valid coord system")
OO=0.
endif
tausp1_par=tmp*OO
else
!
! Constant friction time.
!
tausp1_par=tmp
endif
endif
else if (ldraglaw_epstein_stokes_linear) then
!
! When the particle radius is larger than 9/4 times the mean free path
! of the gas molecules one must use the Stokes drag law rather than the
! Epstein law.
!
! We need here to know the mean free path of the gas molecules:
! lambda = mu_mol/(rhog*sigma_mol)
!
! The quantities are:
! mu_mol = mean molecular weight [=3.9e-24 g for H_2]
! rhog = gas density
! sigma_mol = cross section of gas molecules [=2e-15 cm^2 for H_2]
!
! Actually need to know the mean free path in units of the gas scale
! height H [if H=1]. Inserting the mid-plane expression
! rhog=Sigmag/[sqrt(2*pi)*H]
! gives
! lambda/H = sqrt(2*pi)*mu_mol/(Sigmag*sigma_mol)
! ~= 4.5e-9/Sigmag
! when Sigmag is given in g/cm^2.
!
if (iap==0) then
call fatal_error('get_frictiontime','need particle radius as dynamical variable for Stokes law.')
endif
if (fp(k,iap)<2.25*mean_free_path_gas) then
tausp1_par = 1/(fp(k,iap)*rhopmat)
else
tausp1_par = 1/(fp(k,iap)*rhopmat)*2.25*mean_free_path_gas/fp(k,iap)
endif
!
else if (ldraglaw_epstein_transonic) then
!
! Draw laws for intermediate mach number. This is for pure Epstein drag...
!
call calc_draglaw_parameters(fp,k,uup,p,inx0,tausp1_par)
!
else if (ldraglaw_eps_stk_transonic) then
!
! ...and this is for a linear combination of Esptein and Stokes drag at
! intermediate mach number. Pure Stokes drag is not implemented. (implemented now, see below)
!
call calc_draglaw_parameters(fp,k,uup,p,inx0,tausp1_par,lstokes=.true.)
!
! draglaw_purestokes implemented below, it is simple stokes drag with no
! dependence on partile's reynolds number.
!
elseif (ldraglaw_purestokes) then
call calc_draglaw_purestokes(fp,k,tausp1_par)
!
! stokes drag with schiller_nauman, but without epstein
!
elseif (ldraglaw_stokesschiller) then
if (.not.present(rep)) then
call fatal_error('get_frictiontime','need particle reynolds '// &
'number, rep, to calculate the stokes_schiller drag '// &
'relaxation time!')
else
call calc_draglaw_steadystate(fp,k,rep,1.0,tausp1_par)
endif
elseif (ldraglaw_steadystate) then
if (.not.present(rep)) then
call fatal_error('get_frictiontime','need particle reynolds '// &
'number, rep, to calculate the steady state drag '// &
'relaxation time!')
elseif (.not.present(stocunn)) then
call fatal_error('get_frictiontime','need particle stokes '// &
'cunningham factor, stocunn, to calculate the steady '// &
' state drag relaxation time!')
else
call calc_draglaw_steadystate(fp,k,rep,stocunn,tausp1_par)
endif
!
! Simple drag law, drag force = 1/\taup (v-u) where \taup is an input parameter.
!
else if (ldraglaw_simple) then
! write(*,*)'DM','simple drag'
! Check if we are using multiple or single particle species.
if (npar_species>1) then
jspec=npar_species*(ipar(k)-1)/npar+1
tmp=tausp1_species(jspec)
else
tmp=tausp1
endif
tausp1_par=tmp
endif
!
! Change friction time artificially.
!
if (.not. nochange) then
!
! Increase friction time to avoid very small time-steps where the
! dust-to-gas ratio is high.
!
if (tausg_min/=0.0) then
tausg1_point=tausp1_par*p%epsp(ix0-nghost)
if (tausg1_point>tausg1_max) &
tausp1_par=tausg1_max/p%epsp(ix0-nghost)
endif
!
! Increase friction time linearly with dust density where the dust-to-gas
! ratio is higher than a chosen value. Supposed to mimick decreased cooling
! when the gas follows the dust.
!
if (epsp_friction_increase/=0.0) then
if (p%epsp(ix0-nghost)>epsp_friction_increase) &
tausp1_par=tausp1_par/(p%epsp(ix0-nghost)/epsp_friction_increase)
endif
!
endif
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(ineargrid)
!
endsubroutine get_frictiontime
!***********************************************************************
subroutine calc_draglaw_parameters(fp,k,uup,p,inx0,tausp1_par,lstokes)
!
use EquationOfState, only: rho0,cs0
!
real, dimension (mpar_loc,mparray) :: fp
real, dimension(3) :: uup,duu
type (pencil_case) :: p
real :: tausp1_par,tmp,tmp1
integer :: k, inx0, jspec
real :: kd,fd,mach,mach2,fac,OO
real :: knudsen,reynolds,lambda
real :: inv_particle_radius,kn_crit
logical, optional :: lstokes
logical, save :: lfirstcall
!
! Epstein drag away from the limit of subsonic particle motion. The drag
! force is given by (Schaaf 1963)
!
! Feps=-pi*a**2 * rhog * |Delta(u)| * Delta(u) & (1)
! *[(1+1/m**2+1/(4*m**4))*erf(m) + (1/m+1/(2*m**3)*exp(-m**2)/sqrt(pi)]
!
! where Delta(u) is the relative dust-to-gas velocity (vector)
! and m=|Delta(u)|/cs (scalar) is the relative mach number of the flow
!
! As erf is too cumbersome a function to implement numerically, an interpolation
! between the limits of
!
! subsonic: Feps=-sqrt(128*pi)/3*a**2*rhog*cs*Delta(u) (2)
! supersonic: Feps=-pi*a**2*rhog*|Delta(u)|*Delta(u) (3)
!
! is used, leading to an expression that can be used for arbitrary velocities
! as derived by Kwok (1975).
!
! transonic: Feps=-sqrt(128*pi)/3*a**2*rhog*cs*fd*Delta(u) (4)
!
! where fd=sqrt(1 + 9*pi/128*m**2) (5)
!
! The force Feps is divided by the mass of the particle mp=4/3*pi*a**3*rhopmat
! to yield the acceleration feps=Feps/mp
!
! feps = -sqrt(8/pi)*rhog*cs*fd*Delta(u)/[a*rhopmat] (6)
!
! Epstein drag ceases to work when the particle diameter becomes comparable
! to the mean free path (lambda) of the gas molecules. In this case, the force
! is given by Stokes friction in the viscous case (low dust Reynolds number)
!
! Fsto=-6*pi*a*mu_kin*Delta(u) (7)
!
! where mu_kin is the kinematic viscosity of the gas
!
! mu_kin=1/3*rhog*vth*lambda (8)
!
! and vth=sqrt(8/pi)*cs is the mean thermal velocity of the gas. For high dust
! Reynolds numbers the viscosity if uninmportant and the drag force of the tur-
! bulent flow past the particle is given by Newtonian friction
!
! Fnew=-1.3*pi*a**2*rhog*|Delta(u)|*Delta(u)
!
! The two cases are once again connected by an interpolating factor
!
! F'sto=-6*pi*a*kd*mu_kin*Delta(u)
!
! where kd is a factor that contains the Reynolds number of the flow over the
! particle (defined in the code, some lines below).
!
! The following interpolation then works for flows of arbitrary Knudsen, Mach and Reynolds
! numbers
!
! Fdrag = [Kn'/(Kn'+1)]**2 * Feps + [1/(Kn'+1)]**2 * F'sto
!
! Where Kn'=3*Kn is the critical Knudsen number where the viscous (Stokes) drag and the subsonic
! Epstein drag are equal.
!
! (The discussion above was taken from Paardekooper 2006, Woite & Helling 2003 and Kwok 1975)
!
! In the 2D case, the density rhog is to be replaced by
!
! rhog=Sigmag/[sqrt(2*pi)H]
! =Sigmag*Omega/[sqrt(2*pi)*cs]
!
! which removes the dependence of (6) on cs. We are left with
!
! feps = -2/pi*sigmag*Omega*fd*Delta(u)/[a*rhopmat]
!
! the constant terms are tausp1. The same follows for Stokes drag
!
! Friction time for different species
!
if (npar_species==1) then
tmp=tausp
tmp1=tausp1
else
jspec=npar_species*(ipar(k)-1)/npar+1
tmp=tausp_species(jspec)
tmp1=tausp1_species(jspec)
endif
!
! Relative velocity
!
duu=fp(k,ivpx:ivpz)-uup
!
if (nzgrid==1) then
! then omega is needed
if (ldraglaw_variable) then
!these omegas assume GM=1
if (lcartesian_coords) then
OO=(fp(k,ixp)**2 + fp(k,iyp)**2)**(-0.75)
elseif (lcylindrical_coords) then
OO=fp(k,ixp)**(-1.5)
elseif (lspherical_coords) then
OO=(fp(k,ixp)*sin(fp(k,iyp)))**(-1.5)
else
call fatal_error("calc_draglaw_parameters", &
"no valid coord system")
OO=0.
endif
else
OO=nu_epicycle
endif
endif
!
! Possibility to include the transition from Epstein to Stokes drag
!
if (present(lstokes)) then
!
if (lfirstcall) &
print*, 'get_frictiontime: Epstein-Stokes transonic drag law'
!
! The mach number and the correction fd to flows of arbitrary mach number
!
mach=sqrt((duu(1)**2+duu(2)**2+duu(3)**2)/p%cs2(inx0))
fd=sqrt(1+(9.0*pi/128)*mach**2)
!
! For Stokes drag, the mean free path is needed
!
! lambda = 1/rhog*(mu/sigma_coll)_H2
!
! were mu_H2 is the mean molecular weight of the hydrogen molecule (3.9e-24 g),
! and sigma_coll its cross section (2e-15 cm^2).
! Assume that (mu/sigma_coll) is the input parameter mean_free_path_gas
!
if (mean_free_path_gas==0) then
call fatal_error("calc_draglaw_parameters","You want to use Stokes drag"// &
"but you forgot to set 'mean_free_path_gas' in the *.in files.")
endif
!
if (nzgrid==1) then
!the sqrt(2pi) factor is inside the mean_free_path_gas constant
lambda=mean_free_path_gas * sqrt(p%cs2(inx0))*rho0/(p%rho(inx0)*OO*cs0)
else
lambda=mean_free_path_gas * rho0/p%rho(inx0)
endif
!
! The Knudsen number is the ratio of the mean free path to the particle
! radius, 2s. To keep consistency with the formulation evolving for radius,
! tausp1 is C/(s*rhopmat) where C is 2/pi for 2d runs and sqrt(8/pi) for 3D
! runs (because of the sqrt(2*pi) factor coming from the substitution
! Sigma=rho/(sqrt(2*pi)*H). 's' is the particle radius
if (iap/=0) then
inv_particle_radius=1/fp(k,iap)
else
if (luse_tau_ap) then
! use tausp as the radius (in meter) to make life easier
inv_particle_radius=tmp1
else
if (nzgrid==1) then
inv_particle_radius=0.5*pi*tmp1 !rhopmat=1, particle_radius in meters
else
inv_particle_radius=sqrt(pi/8)*tmp1 !rhopmat=1, particle_radius in meters
endif
endif
endif
!
knudsen=0.5*lambda*inv_particle_radius
!
! The Stokes drag depends non-linearly on
!
! Re = 2*s*rho_g*|delta(v)|/mu_kin
!
reynolds=3*sqrt(pi/8)*mach/knudsen
!
! the Reynolds number of the flow over the particle. It can parameterized by
!
if (reynolds<=500) then
kd=1.0+0.15*reynolds**0.687
elseif ((reynolds>500).and.(reynolds<=1500)) then
kd=3.96e-6*reynolds**2.4
elseif (reynolds>1500) then
kd=0.11*reynolds
else
call fatal_error_local("calc_draglaw_parameters", "'reynolds' seems to be NaN!")
kd=0.
endif
!
! And we finally have the Stokes correction to intermediate Knudsen numbers
! kn_crit is the critical knudsen number where viscous (low reynolds)
! Stokes and subsonic Epstein friction are equal (Woitke & Helling, 2003)
!
kn_crit=3*knudsen
fac=kn_crit/(kn_crit+1)**2 * (kn_crit*fd + kd)
!
else
!
! Only use Epstein drag
!
if (lfirstcall) &
print*,'get_frictiontime: Epstein transonic drag law'
!
mach2=(duu(1)**2+duu(2)**2+duu(3)**2)/p%cs2(inx0)
fd=sqrt(1+(9.0*pi/128)*mach2)
fac=fd
!
endif
!
! Calculate tausp1_par for 2d and 3d cases with and without particle_radius
! as a dynamical variable
!
if (iap/=0) then
if (fp(k,iap)/=0.0) then
if (nzgrid==1) then
tausp1_par= 2*pi_1*OO* &
p%rho(inx0)*fac/(fp(k,iap)*rhopmat)
else
tausp1_par=sqrt(8*pi_1*p%cs2(inx0))*p%rho(inx0)* &
fac/(fp(k,iap)*rhopmat)
endif
endif
else
!normalize to make tausp1 not dependent on cs0 or rho0
!bad because it comes at the expense of evil divisions
if (nzgrid==1) then
if (luse_tau_ap) then
tausp1_par=tmp1*2*pi_1*OO*p%rho(inx0)*fac/(rho0*rhopmat)
else
tausp1_par=tmp1*OO*p%rho(inx0)*fac/ rho0
endif
else
if (luse_tau_ap) then
tausp1_par=tmp1*sqrt(8*pi_1*p%cs2(inx0))*p%rho(inx0)*fac/(rho0*cs0)
else
tausp1_par=tmp1*sqrt(p%cs2(inx0))*p%rho(inx0)*fac/(rho0*cs0)
endif
endif
endif
!
if (lfirstcall) lfirstcall=.false.
!
endsubroutine calc_draglaw_parameters
!***********************************************************************
subroutine collisional_cooling(f,df,fp,dfp,p,ineargrid)
!
! Reduce relative speed between particles due to inelastic collisions.
!
! 23-sep-06/anders: coded
!
use Diagnostics
!
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
!
real, dimension (nx,npar_species,npar_species) :: tau_coll1_species
real, dimension (nx,3,npar_species) :: vvpm_species
real, dimension (nx,npar_species) :: np_species, vpm_species
real, dimension (nx,npar_species) :: tau_coll1_tot
real, dimension (nx,3) :: vvpm
real, dimension (nx) :: vpm, tau_coll1, tausp1m, vcoll, coll_heat
real, dimension (nx) :: rhop_swarm_mn
real, dimension (3) :: deltavp_vec, vbar_jk
real :: deltavp, tau_cool1_par, dt1_cool
real :: tausp1_par, tausp1_parj, tausp1_park, tausp_parj, tausp_park
real :: tausp_parj3, tausp_park3, rhop_swarm_par
integer :: j, k, l, ix0
integer :: ispecies, jspecies
!
if (ldiagnos .or. lentropy .and. lcollisional_heat) coll_heat=0.0
!
! Add collisional cooling of the rms speed.
!
if (lcollisional_cooling_taucool) then
if (npar_imn(imn)/=0) then
vvpm=0.0
do k=k1_imn(imn),k2_imn(imn)
ix0=ineargrid(k,1)
vvpm(ix0-nghost,:) = vvpm(ix0-nghost,:) + fp(k,ivpx:ivpz)
enddo
do l=1,nx
if (p%np(l)>1.0) vvpm(l,:)=vvpm(l,:)/p%np(l)
enddo
!
do k=k1_imn(imn),k2_imn(imn)
ix0=ineargrid(k,1)
dfp(k,ivpx:ivpz) = dfp(k,ivpx:ivpz) - &
taucool1*(fp(k,ivpx:ivpz)-vvpm(ix0-nghost,:))
enddo
!
if (lfirst.and.ldt) dt1_max=max(dt1_max,taucool1/cdtp)
endif
endif
!
! Add collisional cooling of the rms speed, with cooling time-scale based
! on friction time and local rms speed of particles.
!
if (lcollisional_cooling_rms) then
if (npar_imn(imn)/=0) then
! When multiple friction times are present, the average is used for the
! number density in each superparticle.
if (npar_species>1) then
tausp1m=0.0
else
call get_frictiontime(f,fp,p,ineargrid,1,tausp1_par, &
nochange_opt=.true.)
endif
! Need vpm=<|vvp-<vvp>|> to calculate the collisional time-scale.
vvpm=0.0; vpm=0.0
do k=k1_imn(imn),k2_imn(imn)
ix0=ineargrid(k,1)
vvpm(ix0-nghost,:) = vvpm(ix0-nghost,:) + fp(k,ivpx:ivpz)
if (npar_species>1) then
call get_frictiontime(f,fp,p,ineargrid,k,tausp1_par, &
nochange_opt=.true.)
tausp1m(ix0-nghost) = tausp1m(ix0-nghost) + tausp1_par
endif
enddo
do l=1,nx
if (p%np(l)>1.0) then
vvpm(l,:)=vvpm(l,:)/p%np(l)
if (npar_species>1) tausp1m=tausp1m/p%np
endif
enddo
! vpm
do k=k1_imn(imn),k2_imn(imn)
ix0=ineargrid(k,1)
vpm(ix0-nghost) = vpm(ix0-nghost) + &
sqrt( (fp(k,ivpx)-vvpm(ix0-nghost,1))**2 + &
(fp(k,ivpy)-vvpm(ix0-nghost,2))**2 + &
(fp(k,ivpz)-vvpm(ix0-nghost,3))**2 )
enddo
do l=1,nx
if (p%np(l)>1.0) then
vpm(l)=vpm(l)/p%np(l)
endif
enddo
! The collisional time-scale is 1/tau_coll=nd*vrms*sigma_coll.
! Inserting Epstein friction time gives 1/tau_coll=3*rhod/rho*vprms/tauf.
if (npar_species>1) then
tau_coll1=(1.0-coeff_restitution)*p%epsp*vpm*tausp1m
else
tau_coll1=(1.0-coeff_restitution)*p%epsp*vpm*tausp1_par
endif
! Limit inverse time-step of collisional cooling if requested.
if (tau_coll_min>0.0) then
where (tau_coll1>tau_coll1_max) tau_coll1=tau_coll1_max
endif
!
do k=k1_imn(imn),k2_imn(imn)
ix0=ineargrid(k,1)
dfp(k,ivpx:ivpz) = dfp(k,ivpx:ivpz) - &
tau_coll1(ix0-nghost)*(fp(k,ivpx:ivpz)-vvpm(ix0-nghost,:))
if (lcollisional_heat .or. ldiagnos) then
call get_rhopswarm(mp_swarm,fp,k,ineargrid(k,:),rhop_swarm_par)
coll_heat(ix0-nghost) = coll_heat(ix0-nghost) + &
rhop_swarm_par*tau_coll1(ix0-nghost)*&
sum(fp(k,ivpx:ivpz)*(fp(k,ivpx:ivpz)-vvpm(ix0-nghost,:)))
endif
enddo
!
if (lfirst.and.ldt) dt1_max=max(dt1_max,tau_coll1/cdtp)
endif
endif
!
! More advanced collisional cooling model. Collisions are considered for
! every possible two-body process in a grid cell.
!
if (lcollisional_cooling_twobody) then
do l=1,nx
! Collisions between particle k and all other particles in the grid cell.
k=kshepherd(l)
if (k>0) then
! Limit inverse time-step of collisional cooling if requested.
do while (k/=0)
dt1_cool=0.0
call get_frictiontime(f,fp,p,ineargrid,k,tausp1_park, &
nochange_opt=.true.)
tausp_park=1/tausp1_park
tausp_park3=tausp_park**3
j=k
do while (kneighbour(j)/=0)
! Collide with the neighbours of k and their neighbours.
j=kneighbour(j)
call get_frictiontime(f,fp,p,ineargrid,j,tausp1_parj, &
nochange_opt=.true.)
tausp_parj=1/tausp1_parj
tausp_parj3=tausp_parj**3
! Collision velocity.
deltavp_vec=fp(k,ivpx:ivpz)-fp(j,ivpx:ivpz)
deltavp=sqrt( deltavp_vec(1)**2 + deltavp_vec(2)**2 + &
deltavp_vec(3)**2 )
vbar_jk= &
(tausp_parj3*fp(k,ivpx:ivpz)+tausp_park3*fp(j,ivpx:ivpz))/ &
(tausp_parj3+tausp_park3)
! Cooling time-scale.
call get_rhopswarm(mp_swarm,fp,k,ineargrid(k,:),rhop_swarm_par)
tau_cool1_par= &
(1.0-coeff_restitution)* &
rhop_swarm_par*deltavp*(tausp_parj+tausp_park)**2/ &
(tausp_parj3+tausp_park3)
dt1_cool=dt1_cool+tau_cool1_par
! if (tau_coll_min>0.0) then
! if (tau_cool1_par>tau_coll1_max) tau_cool1_par=tau_coll1_max
! endif
dfp(j,ivpx:ivpz) = dfp(j,ivpx:ivpz) - &
tau_cool1_par*(fp(j,ivpx:ivpz)-vbar_jk)
dfp(k,ivpx:ivpz) = dfp(k,ivpx:ivpz) - &
tau_cool1_par*(fp(k,ivpx:ivpz)-vbar_jk)
enddo
if (lfirst.and.ldt) dt1_max=max(dt1_max(l),dt1_cool/cdtp)
! Go through all possible k.
k=kneighbour(k)
enddo
endif
enddo
!
endif
!
! Treat collisions as a drag force that damps the rms speed at the same
! time-scale.
!
if (lcollisional_dragforce_cooling) then
if (npar_imn(imn)/=0) then
vvpm_species=0.0; vpm_species=0.0; np_species=0.0
! Calculate mean velocity and number of particles for each species.
do k=k1_imn(imn),k2_imn(imn)
ix0=ineargrid(k,1)
ispecies=npar_species*(ipar(k)-1)/npar+1
vvpm_species(ix0-nghost,:,ispecies) = &
vvpm_species(ix0-nghost,:,ispecies) + fp(k,ivpx:ivpz)
np_species(ix0-nghost,ispecies) = &
np_species(ix0-nghost,ispecies) + 1.0
enddo
do l=1,nx
do ispecies=1,npar_species
if (np_species(l,ispecies)>1.0) then
vvpm_species(l,:,ispecies)=vvpm_species(l,:,ispecies)/np_species(l,ispecies)
endif
enddo
enddo
! Calculate rms speed for each species.
do k=k1_imn(imn),k2_imn(imn)
ix0=ineargrid(k,1)
ispecies=npar_species*(ipar(k)-1)/npar+1
vpm_species(ix0-nghost,ispecies) = &
vpm_species(ix0-nghost,ispecies) + sqrt( &
(fp(k,ivpx)-vvpm_species(ix0-nghost,1,ispecies))**2 + &
(fp(k,ivpy)-vvpm_species(ix0-nghost,2,ispecies))**2 + &
(fp(k,ivpz)-vvpm_species(ix0-nghost,3,ispecies))**2 )
enddo
do l=1,nx
do ispecies=1,npar_species
if (np_species(l,ispecies)>1.0) then
vpm_species(l,ispecies)=vpm_species(l,ispecies)/np_species(l,ispecies)
endif
enddo
enddo
!
! Collisional drag force time-scale between particles i and j with R_i < R_j.
!
! tau_ji = tau_j^3/(tau_i+tau_j)^2/(deltav_ij/cs*rhoi/rhog)
! tau_ij = tau_ji*rho_j/rho_i
!
do ispecies=1,npar_species; do jspecies=ispecies,npar_species
vcoll= &
sqrt(vpm_species(:,ispecies)**2+vpm_species(:,ispecies)**2 + &
(vvpm_species(:,1,ispecies)-vvpm_species(:,1,jspecies))**2 + &
(vvpm_species(:,2,ispecies)-vvpm_species(:,2,jspecies))**2 + &
(vvpm_species(:,3,ispecies)-vvpm_species(:,3,jspecies))**2)
call get_rhopswarm(mp_swarm,fp,k,m,n,rhop_swarm_mn)
tau_coll1_species(:,jspecies,ispecies) = &
vcoll*np_species(:,ispecies)*rhop_swarm_mn*p%rho1 / ( &
tausp_species(jspecies)**3/ &
(tausp_species(ispecies)+tausp_species(jspecies))**2 )
where (np_species(:,ispecies)/=0.0) &
tau_coll1_species(:,ispecies,jspecies)= &
tau_coll1_species(:,jspecies,ispecies)*np_species(:,jspecies)/np_species(:,ispecies)
enddo; enddo
!
tau_coll1_tot=0.0
do ispecies=1,npar_species; do jspecies=1,npar_species
tau_coll1_tot(:,ispecies)=tau_coll1_tot(:,ispecies)+tau_coll1_species(:,ispecies,jspecies)
enddo; enddo
! Limit inverse time-step of collisional cooling if requested.
if (tau_coll_min>0.0) then
do ispecies=1,npar_species; do l=1,nx
if (tau_coll1_tot(l,ispecies) > tau_coll1_max) then
tau_coll1_species(l,ispecies,:)=tau_coll1_species(l,ispecies,:)* &
tau_coll1_max/tau_coll1_tot(l,ispecies)
endif
enddo; enddo
tau_coll1_tot=0.0
do ispecies=1,npar_species; do jspecies=1,npar_species
tau_coll1_tot(:,ispecies)=tau_coll1_tot(:,ispecies)+tau_coll1_species(:,ispecies,jspecies)
enddo; enddo
endif
if (lfirst.and.ldt) then
do ispecies=1,npar_species
dt1_max=max(dt1_max,tau_coll1_tot(:,ispecies)/cdtp)
enddo
endif
! Add to equation of motion.
do k=k1_imn(imn),k2_imn(imn)
ix0=ineargrid(k,1)
ispecies=npar_species*(ipar(k)-1)/npar+1
do jspecies=1,npar_species
dfp(k,ivpx:ivpz) = dfp(k,ivpx:ivpz) - &
tau_coll1_species(ix0-nghost,ispecies,jspecies)* &
(fp(k,ivpx:ivpz)-vvpm_species(ix0-nghost,:,jspecies))
if (lcollisional_heat .or. ldiagnos) then
call get_rhopswarm(mp_swarm,fp,k,ineargrid(k,:),rhop_swarm_par)
coll_heat(ix0-nghost) = coll_heat(ix0-nghost) + &
rhop_swarm_par* &
tau_coll1_species(ix0-nghost,ispecies,jspecies)* &
sum(fp(k,ivpx:ivpz)*(fp(k,ivpx:ivpz) - &
vvpm_species(ix0-nghost,:,jspecies)))
endif
enddo
enddo
endif
endif
!
! Heating of the gas due to dissipative collisions.
!
if (lentropy .and. lcollisional_heat) &
df(l1:l2,m,n,iss) = df(l1:l2,m,n,iss) + p%rho1*p%TT1*coll_heat
!
! Diagnostics.
if (ldiagnos) then
if (idiag_decollp/=0) &
call sum_mn_name(coll_heat,idiag_decollp)
endif
!
endsubroutine collisional_cooling
!***********************************************************************
subroutine compensate_friction_increase(f,fp,dfp,p,ineargrid)
!
! Compensate for increased friction time in regions of high solids-to-gas
! ratio by applying missing friction force to particles only.
!
! 26-feb-07/anders: coded
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mpar_loc,mparray) :: fp
real, dimension (mpar_loc,mpvar) :: dfp
type (pencil_case) :: p
integer, dimension (mpar_loc,3) :: ineargrid
!
real, dimension (nx,3) :: vvpm
real :: tausp1_par, tausp1_par_mod, tausp1_par_org
integer :: k, l, ix0
!
if (npar_imn(imn)/=0) then
! Calculate mean particle velocity.
vvpm=0.0
do k=k1_imn(imn),k2_imn(imn)
ix0=ineargrid(k,1)
vvpm(ix0-nghost,:) = vvpm(ix0-nghost,:) + fp(k,ivpx:ivpz)
enddo
do l=1,nx
if (p%np(l)>1.0) vvpm(l,:)=vvpm(l,:)/p%np(l)
enddo
!
! Compare actual and modified friction time and apply the difference as
! friction force relative to mean particle velocity.
!
do k=k1_imn(imn),k2_imn(imn)
call get_frictiontime(f,fp,p,ineargrid,k,tausp1_par_mod, &
nochange_opt=.false.)
call get_frictiontime(f,fp,p,ineargrid,k,tausp1_par_org, &
nochange_opt=.true.)
tausp1_par=tausp1_par_org-tausp1_par_mod
ix0=ineargrid(k,1)
dfp(k,ivpx:ivpz) = dfp(k,ivpx:ivpz) - &
tausp1_par*(fp(k,ivpx:ivpz)-vvpm(ix0-nghost,:))
enddo
endif
!
endsubroutine compensate_friction_increase
!***********************************************************************
subroutine calc_gas_velocity_shell_call(k1, uup, fp)
!
use Special, only: special_calc_particles
!
real, dimension (mpar_loc,mparray) :: fp
real, dimension(3) :: uup
integer :: k1
!
vel_call=.true.
uup_shared=fp(k1,ixp:izp)
! call special_calc_particles(fp,ineargrid)
uup=uup_shared
!
endsubroutine calc_gas_velocity_shell_call
!***********************************************************************
subroutine calculate_rms_speed(fp,ineargrid,p)
!
use Diagnostics
!
! Calculate the rms speed dvpm=sqrt(<(vvp-<vvp>)^2>) of the
! particle for diagnostic purposes
!
! 08-04-08/wlad: coded
!
real, dimension (mpar_loc,mparray) :: fp
integer, dimension (mpar_loc,3) :: ineargrid
real,dimension(nx,3) :: vvpm,dvp2m
integer :: inx0,k,l
type (pencil_case) :: p
!
! Initialize the variables
!
vvpm=0.0; dvp2m=0.0
!
! Calculate the average velocity at each cell
! if there are particles in the pencil only
!
if (npar_imn(imn)/=0) then
!
do k=k1_imn(imn),k2_imn(imn)
inx0=ineargrid(k,1)-nghost
vvpm(inx0,:) = vvpm(inx0,:) + fp(k,ivpx:ivpz)
enddo
do l=1,nx
if (p%np(l)>1.0) vvpm(l,:)=vvpm(l,:)/p%np(l)
enddo
!
! Get the residual in quadrature, dvp2m. Need vvpm calculated above.
!
do k=k1_imn(imn),k2_imn(imn)
inx0=ineargrid(k,1)-nghost
dvp2m(inx0,1)=dvp2m(inx0,1)+(fp(k,ivpx)-vvpm(inx0,1))**2
dvp2m(inx0,2)=dvp2m(inx0,2)+(fp(k,ivpy)-vvpm(inx0,2))**2
dvp2m(inx0,3)=dvp2m(inx0,3)+(fp(k,ivpz)-vvpm(inx0,3))**2
enddo
do l=1,nx
if (p%np(l)>1.0) dvp2m(l,:)=dvp2m(l,:)/p%np(l)
enddo
!
endif
!
! Output the diagnostics
!
if (idiag_dvpx2m/=0) call sum_mn_name(dvp2m(:,1),idiag_dvpx2m)
if (idiag_dvpy2m/=0) call sum_mn_name(dvp2m(:,2),idiag_dvpy2m)
if (idiag_dvpz2m/=0) call sum_mn_name(dvp2m(:,3),idiag_dvpz2m)
if (idiag_dvpm/=0) call sum_mn_name(dvp2m(:,1)+dvp2m(:,2)+dvp2m(:,3),&
idiag_dvpm,lsqrt=.true.)
if (idiag_dvpmax/=0) call max_mn_name(dvp2m(:,1)+dvp2m(:,2)+dvp2m(:,3),&
idiag_dvpmax,lsqrt=.true.)
!
endsubroutine calculate_rms_speed
!***********************************************************************
subroutine calc_pencil_rep(fp,rep)
!
! Calculate particle Reynolds numbers
!
! 16-jul-08/kapelrud: coded
!
use Viscosity, only: getnu
!
real, dimension (mpar_loc,mparray), intent(in) :: fp
real,dimension(k1_imn(imn):k2_imn(imn)), intent(inout) :: rep
!
real,dimension(k1_imn(imn):k2_imn(imn)) :: nu
character (len=labellen) :: ivis=''
real :: nu_
integer :: k
!
! Find the kinematic viscosity.
! Check whether we want to override the usual viscosity for the drag law.
!
if (lnu_draglaw) then
nu=nu_draglaw
else
call getnu(nu_input=nu_,IVIS=ivis)
if (ivis=='nu-const') then
nu=nu_
elseif (ivis=='nu-mixture') then
nu=interp_nu
elseif (ivis=='rho-nu-const') then
nu=nu_/interp_rho(k1_imn(imn):k2_imn(imn))
elseif (ivis=='sqrtrho-nu-const') then
nu=nu_/sqrt(interp_rho(k1_imn(imn):k2_imn(imn)))
elseif (ivis=='nu-therm') then
nu=nu_*sqrt(interp_TT(k1_imn(imn):k2_imn(imn)))
elseif (ivis=='mu-therm') then
nu=nu_*sqrt(interp_TT(k1_imn(imn):k2_imn(imn)))&
/interp_rho(k1_imn(imn):k2_imn(imn))
else
call fatal_error('calc_pencil_rep','No such ivis!')
endif
endif
!
if (maxval(nu) == 0.0) call fatal_error('calc_pencil_rep', 'nu (kinematic visc.) must be non-zero!')
!
do k=k1_imn(imn),k2_imn(imn)
rep(k) = 2.0 * sqrt(sum((interp_uu(k,:) - fp(k,ivpx:ivpz))**2)) / nu(k)
enddo
!
if (lparticles_radius) then
rep = rep * fp(k1_imn(imn):k2_imn(imn),iap)
elseif (particle_radius > 0.0) then
rep = rep * particle_radius
else
call fatal_error('calc_pencil_rep', 'unable to calculate the particle Reynolds number without a particle radius. ')
endif
!
endsubroutine calc_pencil_rep
!***********************************************************************
subroutine calc_stokes_cunningham(fp,stocunn)
!
! Calculate the Stokes-Cunningham factor
!
! 12-aug-08/kapelrud: coded
!
use Particles_radius
!
real, dimension (mpar_loc,mparray) :: fp
real,dimension(k1_imn(imn):k2_imn(imn)) :: stocunn
!
real :: dia
integer :: k
!
do k=k1_imn(imn),k2_imn(imn)
!
! Particle diameter
!
dia=2.0*fp(k,iap)
!
if (lstocunn1) then
stocunn(k) = 1.
else
stocunn(k)=1.+2.*mean_free_path_gas/dia* &
(1.257+0.4*exp(-0.55*dia/mean_free_path_gas))
endif
!
enddo
!
endsubroutine calc_stokes_cunningham
!**********************************************************************
subroutine calc_added_mass_beta(fp,k,added_mass_beta)
!
real, dimension (mpar_loc,mparray), intent(in) :: fp
integer, intent(in) :: k
real, intent(out) :: added_mass_beta
!
! beta for added mass according to beta=3rho_fluid/(2rho_part+rho_fluid)
! problem: we would have to calculate beta every time for every particle
added_mass_beta=3*interp_rho(k)/(2*rhopmat+interp_rho(k))
!
end subroutine calc_added_mass_beta
!***********************************************************************
subroutine calc_draglaw_purestokes(fp,k,tausp1_par)
!
! Calculate relaxation time for particles under Pure Stokes drag
!
! 6-Aug-15/nils+dhruba: coded
!
use Viscosity, only: getnu
use Particles_radius
!
real, dimension (mpar_loc,mparray), intent(in) :: fp
integer, intent(in) :: k
real, intent(out) :: tausp1_par
!
character (len=labellen) :: ivis=''
real :: dia,nu,nu_
!
! Find the kinematic viscosity
!
call getnu(nu_input=nu_,ivis=ivis)
if (ivis=='nu-const') then
nu=nu_
elseif (ivis=='nu-mixture') then
nu=interp_nu(k)
elseif (ivis=='rho-nu-const') then
nu=nu_/interp_rho(k)
elseif (ivis=='sqrtrho-nu-const') then
nu=nu_/sqrt(interp_rho(k))
elseif (ivis=='nu-therm') then
nu=nu_*sqrt(interp_TT(k))
elseif (ivis=='mu-therm') then
nu=nu_*sqrt(interp_TT(k))&
/interp_rho(k)
else
call fatal_error('calc_draglaw_purestokes','No such ivis!')
endif
!
! Particle diameter
!
if (.not.lparticles_radius) then
call fatal_error('calc_draglaw_purestokes', &
'need particles_radius module to calculate the relaxation time!')
endif
!
dia=2.0*fp(k,iap)
!
! Relaxation time:
!
tausp1_par=18.0*nu/((rhopmat/interp_rho(k))*dia**2)
!
endsubroutine calc_draglaw_purestokes
!***********************************************************************
subroutine calc_draglaw_steadystate(fp,k,rep,stocunn,tausp1_par)
!
! Calculate relaxation time for particles under steady state drag.
!
! 15-jul-08/kapelrud: coded
!
use Viscosity, only: getnu
use Particles_radius
!
real, dimension (mpar_loc,mparray), intent(in) :: fp
integer, intent(in) :: k
real, intent(in) :: rep, stocunn
real, intent(out) :: tausp1_par
!
character (len=labellen) :: ivis=''
real :: cdrag,dia,nu,nu_
!
! Find the kinematic viscosity.
! Check whether we want to override the usual viscosity for the drag law.
!
if (lnu_draglaw) then
nu=nu_draglaw
else
!
! Use usual viscosity for the drag law.
!
call getnu(nu_input=nu_,ivis=ivis)
if (ivis=='nu-const') then
nu=nu_
elseif (ivis=='nu-mixture') then
nu=interp_nu(k)
elseif (ivis=='rho-nu-const') then
nu=nu_/interp_rho(k)
elseif (ivis=='sqrtrho-nu-const') then
nu=nu_/sqrt(interp_rho(k))
elseif (ivis=='nu-therm') then
nu=nu_*sqrt(interp_TT(k))
elseif (ivis=='mu-therm') then
nu=nu_*sqrt(interp_TT(k))&
/interp_rho(k)
else
call fatal_error('calc_draglaw_steadystate','No such ivis!')
endif
endif
!
! Particle diameter
!
if (.not.lparticles_radius) then
call fatal_error('calc_draglaw_steadystate', &
'need particles_radius module to calculate the relaxation time!')
endif
!
dia=2.0*fp(k,iap)
!
! Calculate drag coefficent pre-factor:
!
if (rep<1) then
cdrag=1.0
elseif (rep>1000) then
cdrag=0.44*rep/24.0
else
cdrag=(1.+0.15*rep**0.687)
endif
!
! Relaxation time:
!
tausp1_par=18.0*cdrag*nu/((rhopmat/interp_rho(k))*stocunn*dia**2)
!
endsubroutine calc_draglaw_steadystate
!***********************************************************************
subroutine calc_brownian_force(fp,k,ineark,stocunn,force)
!
! Calculate the Brownian force contribution due to the random thermal motions
! of the gas molecules.
!
! 28-jul-08/kapelrud: coded
!
use General, only: normal_deviate
use Viscosity, only: getnu
!
real, dimension (mpar_loc,mparray), intent(in) :: fp
integer, intent(in) :: k
integer, dimension(3) :: ineark
real :: stocunn
real, dimension(3), intent(out) :: force
!
character (len=labellen) :: ivis=''
real :: Szero,dia,TT,rhop_swarm_par,nu,nu_
!
! Find kinematic viscosity
!
call getnu(nu_input=nu_,IVIS=ivis)
if (ivis=='nu-const') then
nu=nu_
elseif (ivis=='nu-mixture') then
nu=interp_nu(k)
elseif (ivis=='rho-nu-const') then
nu=nu_/interp_rho(k)
elseif (ivis=='sqrtrho-nu-const') then
nu=nu_/sqrt(interp_rho(k))
elseif (ivis=='nu-therm') then
nu=nu_*sqrt(interp_TT(k))
elseif (ivis=='mu-therm') then
nu=nu_*sqrt(interp_TT(k))&
/interp_rho(k)
else
call fatal_error('calc_brownian_force','No such ivis!')
endif
!
! Particle diameter:
!
dia=2.0*fp(k,iap)
!
! Get zero mean, unit variance Gaussian random numbers:
!
call normal_deviate(force(1))
call normal_deviate(force(2))
call normal_deviate(force(3))
!
if (interp%lTT) then
TT=interp_TT(k)
else
TT=brownian_T0
endif
!
call get_rhopswarm(mp_swarm,fp,k,ineark,rhop_swarm_par)
!
Szero=216*nu*k_B*TT*pi_1/ &
(dia**5*stocunn*rhop_swarm_par**2/interp_rho(k))
!
! https://en.wikipedia.org/wiki/Brownian_motion
! .3 * urms * lmfp = D=kB*T/(6pi*eta*a)
! .3 * urms^2 * tau = D=kB*T/(6pi*eta*a)
!
!print*,'AXEL: nu,k_B,dia,stocunn,rhop_swarm_par,interp_rho(k)=', nu,k_B,dia,stocunn,rhop_swarm_par,interp_rho(k)
!
! du/dt = sqrt(dt)
! .3 * urms^2 * tau = D=kB*T/(6pi*eta*a)
!
if (dt==0.0) then
force=0.0
else
force=force*sqrt(Szero/dt)
endif
!
endsubroutine calc_brownian_force
!***********************************************************************
subroutine calc_thermophoretic_force(fp,k,ineark,force)
!
! Calculate the Thermophoretic force due to a temperature gradient in the gas.
!
! Henrik Lutro, testing
!
use Viscosity, only: getnu
!
real, dimension (mpar_loc,mparray), intent(in) :: fp
integer, intent(in) :: k
integer, dimension(3) :: ineark
real, dimension(3), intent(out) :: force
!
real, dimension(3) :: temp_grad
real TT,mu,nu_,Kn,phi_therm,mass_p
real Ktc,Ce,Cm,Cint
character (len=labellen) :: ivis=''
!
call keep_compiler_quiet(ineark)
!
Ktc=1.10
Cm=1.13
Ce=2.17
if (interp%lTT) then
TT=interp_TT(k)
else
TT=thermophoretic_T0
endif
!Find dynamic viscosity
call getnu(nu_input=nu_,ivis=ivis)
if (ivis=='nu-const') then
mu=nu_*interp_rho(k)
elseif (ivis=='nu-mixture') then
mu=interp_nu(k)*interp_rho(k)
elseif (ivis=='rho-nu-const') then
mu=nu_
elseif (ivis=='sqrtrho-nu-const') then
mu=nu_*sqrt(interp_rho(k))
elseif (ivis=='nu-therm') then
mu=nu_*interp_rho(k)*sqrt(TT)
elseif (ivis=='mu-therm') then
mu=nu_*sqrt(TT)
else
call fatal_error('calc_thermophoretic_force','No such ivis!')
endif
Cint=0.5
if (interp%lgradTT) then
temp_grad=interp_gradTT(k,:)
else
temp_grad=temp_grad0
endif
!
Kn=mean_free_path_gas/fp(k,iap)
!
mass_p=(4.0*pi/3.0)*rhopmat*fp(k,iap)**3
if (thermophoretic_eq=='near_continuum') then
phi_therm=-9*pi/cond_ratio
elseif (thermophoretic_eq=='transition') then
phi_therm=-12.0*pi*(Ktc*(1.0+cond_ratio*Ce*Kn)+3.0*Cm*Kn*(1.0-cond_ratio+cond_ratio*Ce*Kn))&
/((1.0+3.0*Kn*exp(-Cint/Kn))*(1.0+3.0*Cm*Kn)*(2.0+cond_ratio+2.0*cond_ratio*Ce*Kn))
elseif (thermophoretic_eq=='free_molecule') then
phi_therm=0.0
else
call fatal_error('calc_pencil_rep','No thermoporetic range chosen')
endif
!
force=(fp(k,iap)*temp_grad*mu**2*phi_therm)/(TT*interp_rho(k)*mass_p)
!
endsubroutine calc_thermophoretic_force
!***********************************************************************
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)
!
! If we have bubbles, the advective derivative has to be saved in
! an auxiliary variable
! COMMENT: This would be better to do in a step between registering and
! initializing. Such a hook does not exist at the moment.
!
if (lbubble) ladv_der_as_aux=.true.
!
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 dust particle variables.
!
! 01-jan-06/anders: coded
!
use Power_spectrum, only: power_1d
!
real, dimension (mx,my,mz,mfarray) :: f
!
if (lpar_spec) call power_1d(f,'p',0,irhop)
!
endsubroutine powersnap_particles
!***********************************************************************
subroutine rprint_particles(lreset,lwrite)
!
! Read and register print parameters relevant for particles.
!
! 29-dec-04/anders: coded
!
use Diagnostics
use General, only: itoa
!
logical :: lreset
logical, optional :: lwrite
!
integer :: iname,inamez,inamey,inamex,inamexy,inamexz,inamer,inamerz
integer :: k
logical :: lwr
character (len=intlen) :: srad
!
! Write information to index.pro.
!
lwr = .false.
if (present(lwrite)) lwr=lwrite
!
if (lwr) then
write(3,*) 'ixp=', ixp
write(3,*) 'iyp=', iyp
write(3,*) 'izp=', izp
write(3,*) 'ivpx=', ivpx
write(3,*) 'ivpy=', ivpy
write(3,*) 'ivpz=', ivpz
write(3,*) 'inp=', inp
write(3,*) 'irhop=', irhop
write(3,*) 'iupx=', iupx
write(3,*) 'iupy=', iupy
write(3,*) 'iupz=', iupz
write(3,*) 'ifgx=', ifgx
write(3,*) 'ifgy=', ifgy
write(3,*) 'ifgz=', ifgz
endif
!
! Reset everything in case of reset.
!
if (lreset) then
idiag_xpm=0; idiag_ypm=0; idiag_zpm=0
idiag_vrelpabsm=0
idiag_xp2m=0; idiag_yp2m=0; idiag_zp2m=0; idiag_rpm=0; idiag_rp2m=0
idiag_vpxm=0; idiag_vpym=0; idiag_vpzm=0; idiag_urel=0
idiag_vpxvpym=0; idiag_vpxvpzm=0; idiag_vpyvpzm=0
idiag_vpx2m=0; idiag_vpy2m=0; idiag_vpz2m=0; idiag_ekinp=0
idiag_vpxmax=0; idiag_vpymax=0; idiag_vpzmax=0
idiag_vpxmin=0; idiag_vpymin=0; idiag_vpzmin=0
idiag_rhopvpxm=0; idiag_rhopvpym=0; idiag_rhopvpzm=0; idiag_rhopvpysm=0
idiag_rhopvpxt=0; idiag_rhopvpyt=0; idiag_rhopvpzt=0
idiag_lpxm=0; idiag_lpym=0; idiag_lpzm=0
idiag_lpx2m=0; idiag_lpy2m=0; idiag_lpz2m=0
idiag_npm=0; idiag_np2m=0; idiag_npmax=0; idiag_npmin=0
idiag_dtdragp=0; idiag_dedragp=0
idiag_rhopm=0; idiag_rhoprms=0; idiag_rhop2m=0; idiag_rhopmax=0
idiag_rhopmin=0; idiag_decollp=0; idiag_rhopmphi=0
idiag_epspmin=0; idiag_epspmax=0; idiag_epspm=0;
idiag_nparmin=0; idiag_nparmax=0; idiag_nparsum=0
idiag_nmigmax=0; idiag_nmigmmax=0; idiag_mpt=0
idiag_npmx=0; idiag_npmy=0; idiag_npmz=0; idiag_epotpm=0
idiag_rhopmx=0; idiag_rhopmy=0; idiag_rhopmz=0
idiag_rhop2mx=0; idiag_rhop2my=0; idiag_rhop2mz=0
idiag_epspmx=0; idiag_epspmy=0; idiag_epspmz=0
idiag_rhopmxy=0; idiag_rhopmxz=0; idiag_rhopmr=0
idiag_sigmap = 0
idiag_dvpx2m=0; idiag_dvpy2m=0; idiag_dvpz2m=0
idiag_dvpmax=0; idiag_dvpm=0; idiag_nparpmax=0
idiag_eccpxm=0; idiag_eccpym=0; idiag_eccpzm=0
idiag_eccpx2m=0; idiag_eccpy2m=0; idiag_eccpz2m=0
idiag_npargone=0; idiag_vpyfull2m=0; idiag_deshearbcsm=0
idiag_npmxy=0; idiag_vprms=0
idiag_npvzmz=0; idiag_npvz2mz=0; idiag_nptz=0; idiag_Shm=0
idiag_npuzmz=0
endif
!
! Run through all possible names that may be listed in print.in.
!
if (lroot .and. ip<14) print*,'rprint_particles: run through parse list'
do iname=1,nname
call parse_name(iname,cname(iname),cform(iname),'nparsum',idiag_nparsum)
call parse_name(iname,cname(iname),cform(iname),'nparmin',idiag_nparmin)
call parse_name(iname,cname(iname),cform(iname),'nparmax',idiag_nparmax)
call parse_name(iname,cname(iname),cform(iname),'nparpmax',idiag_nparpmax)
call parse_name(iname,cname(iname),cform(iname),'xpm',idiag_xpm)
call parse_name(iname,cname(iname),cform(iname),'vrelpabsm',idiag_vrelpabsm)
call parse_name(iname,cname(iname),cform(iname),'ypm',idiag_ypm)
call parse_name(iname,cname(iname),cform(iname),'zpm',idiag_zpm)
call parse_name(iname,cname(iname),cform(iname),'xp2m',idiag_xp2m)
call parse_name(iname,cname(iname),cform(iname),'yp2m',idiag_yp2m)
call parse_name(iname,cname(iname),cform(iname),'zp2m',idiag_zp2m)
call parse_name(iname,cname(iname),cform(iname),'rpm',idiag_rpm)
call parse_name(iname,cname(iname),cform(iname),'rp2m',idiag_rp2m)
call parse_name(iname,cname(iname),cform(iname),'urel',idiag_urel)
call parse_name(iname,cname(iname),cform(iname),'vpxm',idiag_vpxm)
call parse_name(iname,cname(iname),cform(iname),'vpym',idiag_vpym)
call parse_name(iname,cname(iname),cform(iname),'vpzm',idiag_vpzm)
call parse_name(iname,cname(iname),cform(iname),'vpxvpym',idiag_vpxvpym)
call parse_name(iname,cname(iname),cform(iname),'vpxvpzm',idiag_vpxvpzm)
call parse_name(iname,cname(iname),cform(iname),'vpyvpzm',idiag_vpyvpzm)
call parse_name(iname,cname(iname),cform(iname),'vpx2m',idiag_vpx2m)
call parse_name(iname,cname(iname),cform(iname),'vpy2m',idiag_vpy2m)
call parse_name(iname,cname(iname),cform(iname),'vpz2m',idiag_vpz2m)
call parse_name(iname,cname(iname),cform(iname),'ekinp',idiag_ekinp)
call parse_name(iname,cname(iname),cform(iname),'vpxmax',idiag_vpxmax)
call parse_name(iname,cname(iname),cform(iname),'vpymax',idiag_vpymax)
call parse_name(iname,cname(iname),cform(iname),'vpzmax',idiag_vpzmax)
call parse_name(iname,cname(iname),cform(iname),'vpxmin',idiag_vpxmin)
call parse_name(iname,cname(iname),cform(iname),'vpymin',idiag_vpymin)
call parse_name(iname,cname(iname),cform(iname),'vpzmin',idiag_vpzmin)
call parse_name(iname,cname(iname),cform(iname),'vpmax',idiag_vpmax)
call parse_name(iname,cname(iname),cform(iname),'rhopvpxm',idiag_rhopvpxm)
call parse_name(iname,cname(iname),cform(iname),'rhopvpym',idiag_rhopvpym)
call parse_name(iname,cname(iname),cform(iname),'rhopvpzm',idiag_rhopvpzm)
call parse_name(iname,cname(iname),cform(iname),'rhopvpysm',idiag_rhopvpysm)
call parse_name(iname,cname(iname),cform(iname),'rhopvpxt',idiag_rhopvpxt)
call parse_name(iname,cname(iname),cform(iname),'rhopvpyt',idiag_rhopvpyt)
call parse_name(iname,cname(iname),cform(iname),'rhopvpzt',idiag_rhopvpzt)
call parse_name(iname,cname(iname),cform(iname),'rhopvpysm',idiag_rhopvpysm)
call parse_name(iname,cname(iname),cform(iname),'lpxm',idiag_lpxm)
call parse_name(iname,cname(iname),cform(iname),'lpym',idiag_lpym)
call parse_name(iname,cname(iname),cform(iname),'lpzm',idiag_lpzm)
call parse_name(iname,cname(iname),cform(iname),'lpx2m',idiag_lpx2m)
call parse_name(iname,cname(iname),cform(iname),'lpy2m',idiag_lpy2m)
call parse_name(iname,cname(iname),cform(iname),'lpz2m',idiag_lpz2m)
call parse_name(iname,cname(iname),cform(iname),'eccpxm',idiag_eccpxm)
call parse_name(iname,cname(iname),cform(iname),'eccpym',idiag_eccpym)
call parse_name(iname,cname(iname),cform(iname),'eccpzm',idiag_eccpzm)
call parse_name(iname,cname(iname),cform(iname),'eccpx2m',idiag_eccpx2m)
call parse_name(iname,cname(iname),cform(iname),'eccpy2m',idiag_eccpy2m)
call parse_name(iname,cname(iname),cform(iname),'eccpz2m',idiag_eccpz2m)
call parse_name(iname,cname(iname),cform(iname),'dtdragp',idiag_dtdragp)
call parse_name(iname,cname(iname),cform(iname),'npm',idiag_npm)
call parse_name(iname,cname(iname),cform(iname),'np2m',idiag_np2m)
call parse_name(iname,cname(iname),cform(iname),'npmax',idiag_npmax)
call parse_name(iname,cname(iname),cform(iname),'npmin',idiag_npmin)
call parse_name(iname,cname(iname),cform(iname),'rhopm',idiag_rhopm)
call parse_name(iname,cname(iname),cform(iname),'rhoprms',idiag_rhoprms)
call parse_name(iname,cname(iname),cform(iname),'rhop2m',idiag_rhop2m)
call parse_name(iname,cname(iname),cform(iname),'rhopmin',idiag_rhopmin)
call parse_name(iname,cname(iname),cform(iname),'rhopmax',idiag_rhopmax)
call parse_name(iname,cname(iname),cform(iname),'epspm',idiag_epspm)
call parse_name(iname,cname(iname),cform(iname),'epspmin',idiag_epspmin)
call parse_name(iname,cname(iname),cform(iname),'epspmax',idiag_epspmax)
call parse_name(iname,cname(iname),cform(iname),'rhopmphi',idiag_rhopmphi)
call parse_name(iname,cname(iname),cform(iname),'nmigmax',idiag_nmigmax)
call parse_name(iname,cname(iname),cform(iname),'nmigmmax',idiag_nmigmmax)
call parse_name(iname,cname(iname),cform(iname),'mpt',idiag_mpt)
call parse_name(iname,cname(iname),cform(iname),'dvpx2m',idiag_dvpx2m)
call parse_name(iname,cname(iname),cform(iname),'dvpy2m',idiag_dvpy2m)
call parse_name(iname,cname(iname),cform(iname),'dvpz2m',idiag_dvpz2m)
call parse_name(iname,cname(iname),cform(iname),'dvpm',idiag_dvpm)
call parse_name(iname,cname(iname),cform(iname),'dvpmax',idiag_dvpmax)
call parse_name(iname,cname(iname),cform(iname),'dedragp',idiag_dedragp)
call parse_name(iname,cname(iname),cform(iname),'decollp',idiag_decollp)
call parse_name(iname,cname(iname),cform(iname),'epotpm',idiag_epotpm)
call parse_name(iname,cname(iname),cform(iname),'npargone',idiag_npargone)
call parse_name(iname,cname(iname),cform(iname),'vpyfull2m',idiag_vpyfull2m)
call parse_name(iname,cname(iname),cform(iname),'vprms',idiag_vprms)
call parse_name(iname,cname(iname),cform(iname),'Shm',idiag_Shm)
call parse_name(iname,cname(iname),cform(iname),'deshearbcsm',idiag_deshearbcsm)
enddo
!
! Check for those quantities for which we want x-averages.
!
do inamex=1,nnamex
call parse_name(inamex,cnamex(inamex),cformx(inamex),'npmx',idiag_npmx)
call parse_name(inamex,cnamex(inamex),cformx(inamex),'rhopmx',idiag_rhopmx)
call parse_name(inamex,cnamex(inamex),cformx(inamex),'rhop2mx',idiag_rhop2mx)
call parse_name(inamex,cnamex(inamex),cformx(inamex),'epspmx',idiag_epspmx)
enddo
!
! Check for those quantities for which we want y-averages.
!
do inamey=1,nnamey
call parse_name(inamey,cnamey(inamey),cformy(inamey),'npmy',idiag_npmy)
call parse_name(inamey,cnamey(inamey),cformy(inamey),'rhopmy',idiag_rhopmy)
call parse_name(inamey,cnamey(inamey),cformy(inamey),'rhop2my',idiag_rhop2my)
call parse_name(inamey,cnamey(inamey),cformy(inamey),'epspmy',idiag_epspmy)
enddo
!
! Check for those quantities for which we want z-averages.
!
do inamez=1,nnamez
call parse_name(inamez,cnamez(inamez),cformz(inamez),'npmz',idiag_npmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'rhopmz',idiag_rhopmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'rhop2mz',idiag_rhop2mz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'epspmz',idiag_epspmz)
do k=1,ndustrad
srad=itoa(k)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'npvzmz'//trim(srad),idiag_npvzmz(k))
call parse_name(inamez,cnamez(inamez),cformz(inamez),'npvz2mz'//trim(srad),idiag_npvz2mz(k))
call parse_name(inamez,cnamez(inamez),cformz(inamez),'npuzmz'//trim(srad),idiag_npuzmz(k))
call parse_name(inamez,cnamez(inamez),cformz(inamez),'nptz'//trim(srad),idiag_nptz(k))
enddo
enddo
!
! Check for those quantities for which we want xy-averages.
!
do inamexy=1,nnamexy
call parse_name(inamexy,cnamexy(inamexy),cformxy(inamexy),'npmxy',idiag_npmxy)
call parse_name(inamexy,cnamexy(inamexy),cformxy(inamexy),'rhopmxy',idiag_rhopmxy)
call parse_name(inamexy,cnamexy(inamexy),cformxy(inamexy),'sigmap',idiag_sigmap)
enddo
!
! Check for those quantities for which we want xz-averages.
!
do inamexz=1,nnamexz
call parse_name(inamexz,cnamexz(inamexz),cformxz(inamexz),'rhopmxz',idiag_rhopmxz)
enddo
!
! Check for those quantities for which we want phiz-averages.
!
do inamer=1,nnamer
call parse_name(inamer,cnamer(inamer),cformr(inamer),'rhopmr',idiag_rhopmr)
enddo
!
! Check for those quantities for which we want phi-averages.
!
do inamerz=1,nnamerz
call parse_name(inamerz,cnamerz(inamerz),cformrz(inamerz),'rhopmphi',idiag_rhopmphi)
enddo
!
endsubroutine rprint_particles
!***********************************************************************
subroutine particles_final_clean_up
!
! cleanup (dummy)
!
print*,'particles_tracer: Nothing to clean up'
endsubroutine particles_final_clean_up
!***********************************************************************
subroutine periodic_boundcond_on_aux(f)
!
!
! Impose periodic boundary condition on gradu as auxiliary variable
!
use Boundcond, only: set_periodic_boundcond_on_aux
!
real, dimension(mx,my,mz,mfarray), intent(in) :: f
!
if (lparticles_grad) then
if (igradu .ne. 0) then
call set_periodic_boundcond_on_aux(f,igradu11)
call set_periodic_boundcond_on_aux(f,igradu12)
call set_periodic_boundcond_on_aux(f,igradu13)
call set_periodic_boundcond_on_aux(f,igradu21)
call set_periodic_boundcond_on_aux(f,igradu22)
call set_periodic_boundcond_on_aux(f,igradu23)
call set_periodic_boundcond_on_aux(f,igradu31)
call set_periodic_boundcond_on_aux(f,igradu32)
call set_periodic_boundcond_on_aux(f,igradu33)
else
call fatal_error('periodic_boundcond_on_aux','particles_grad demands igradu ne 0')
endif
endif
!
endsubroutine periodic_boundcond_on_aux
!***********************************************************************
subroutine calc_relative_velocity(f,fp,ineargrid)
!
use Diagnostics
!
real, dimension(mx,my,mz,mfarray), intent(in) :: f
real, dimension(mpar_loc,mparray), intent(in) :: fp
real, dimension(3) :: uup,rel_vel_sing
real, dimension(:), allocatable :: rel_vel
integer :: k
integer, dimension(mpar_loc,3), intent(in) :: ineargrid
!
! Calculate particle relative velocity
!
allocate(rel_vel(npar_loc))
!
rel_vel = 0.0
!
do k = 1,npar_loc
call interpolate_linear(f,iux,iuz,fp(k,ixp:izp),uup,ineargrid(k,:),0,0)
rel_vel_sing = (fp(k,ivpx:ivpz)-uup)**2
rel_vel(k) = sqrt(sum(rel_vel_sing))
enddo
!
call sum_par_name(rel_vel(1:npar_loc),idiag_vrelpabsm)
if (allocated(rel_vel)) deallocate(rel_vel)
!
endsubroutine calc_relative_velocity
!***********************************************************************
subroutine diffuse_backreaction(f,df)
!
use Boundcond
use Mpicomm
!
integer :: i
real, dimension(mx,my,mz,mfarray), intent(inout) :: f
real, dimension(mx,my,mz,mvar), intent(inout) :: df
!
! passive scalar consumption
!
if (.not. lpencil_check_at_work) then
if (ldiffuse_passive .and. lpscalar_sink) then
do i = 1,istep_pass
call boundconds_x(f,idlncc,idlncc)
call initiate_isendrcv_bdry(f,idlncc,idlncc)
call finalize_isendrcv_bdry(f,idlncc,idlncc)
call boundconds_y(f,idlncc,idlncc)
call boundconds_z(f,idlncc,idlncc)
call diffuse_interaction(f(:,:,:,idlncc),ldiff_pass,.False.,rdiffconst_pass)
enddo
df(l1:l2,m1:m2,n1:n2,ilncc) = df(l1:l2,m1:m2,n1:n2,ilncc) + &
f(l1:l2,m1:m2,n1:n2,idlncc)
f(:,:,:,idlncc) = 0.0
endif
!
! particle fluid dragforce
!
if (ldiffuse_dragf .and. ldragforce_gas_par) then
do i = 1,istep_dragf
call boundconds_x(f,idfx,idfz)
call initiate_isendrcv_bdry(f,idfx,idfz)
call finalize_isendrcv_bdry(f,idfx,idfz)
call boundconds_y(f,idfx,idfz)
call boundconds_z(f,idfx,idfz)
call diffuse_interaction(f(:,:,:,idfx),ldiff_dragf,.False.,rdiffconst_dragf)
call diffuse_interaction(f(:,:,:,idfy),ldiff_dragf,.False.,rdiffconst_dragf)
call diffuse_interaction(f(:,:,:,idfz),ldiff_dragf,.False.,rdiffconst_dragf)
enddo
df(l1:l2,m1:m2,n1:n2,iux:iuz) = df(l1:l2,m1:m2,n1:n2,iux:iuz) + &
f(l1:l2,m1:m2,n1:n2,idfx:idfz)
f(:,:,:,idfx:idfz) = 0.0
endif
endif
!
endsubroutine diffuse_backreaction
!***********************************************************************
endmodule Particles