! $Id$
!
! This modules takes care of instantaneous coagulation, shattering,
! erosion, and bouncing of superparticles.
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: lparticles_coagulation = .true.
!
! MVAR CONTRIBUTION 0
! MAUX CONTRIBUTION 0
!
!***************************************************************
module Particles_coagulation
!
use Cdata
use Cparam
use Messages
use Particles_cdata
use Particles_map
use Particles_mpicomm
use Particles_sub
use General, only: keep_compiler_quiet
!
implicit none
!
include 'particles_coagulation.h'
!
real :: cdtpcoag=0.2, cdtpcoag1=5.0
real :: kernel_cst=1.0, kernel_lin=1.0, kernel_pro=1.0
real :: four_pi_rhopmat_over_three2=0.0
real :: three_over_four_pi_rhopmat=0.0
real :: GNewton=6.67428e-11, deltav_grav_floor=0.0
real :: critical_mass_ratio_sticking=1.0
real :: minimum_particle_mass=0.0, minimum_particle_radius=0.0
real, pointer :: rhs_poisson_const
real :: tstart_droplet_coagulation=impossible
logical :: ldroplet_coagulation_runtime=.false.
logical :: lcoag_simultaneous=.false., lnoselfcollision=.true.
logical :: lshear_in_vp=.true.
logical :: lkernel_test=.false., lconstant_kernel_test=.false.
logical :: llinear_kernel_test=.false., lproduct_kernel_test=.false.
logical :: lgravitational_cross_section=.false.
logical :: lzsomdullemond=.false. ! use Zsom and Dullemond method
logical :: lconstant_deltav=.false. ! use constant relative velocity
logical :: lmaxwell_deltav=.false. ! use maxwellian relative velocity
logical :: ldroplet_coagulation=.false., normal_coagulation=.false.
logical :: lcollision_output=.false., lcollision_output_swapped=.false.,luser_random_number_wrapper=.false.
logical :: lrelabelling=.false.
logical :: kernel_output=.false., radius_output=.false.
logical :: sphericalKernel=.false.
logical :: lremove_particle=.false., lremove_particle2=.false.
character (len=labellen) :: droplet_coagulation_model='standard'
!
real, dimension(:,:), allocatable :: r_ik_mat, cum_func_sec_ik
real, dimension(:), allocatable :: r_i_tot, cum_func_first_i
real :: r_total = 0.0
real :: delta_r = 0.0
real :: deltav = 1.0 ! relative velocity
real :: maxwell_param = 1.0 ! alpha parameter for maxwell distribution
real :: rdifference = 1.0
!
integer :: idiag_ncoagpm=0, idiag_ncoagpartpm=0, idiag_dt1_coag_par=0
!
real :: deltad = 1., a0 = 1.
real :: r1, r2, r3, r4, r5, r6, r7, r8, r_diff
integer :: idiag_k100_100, idiag_k100_80, idiag_k100_60, idiag_k100_50, &
idiag_k100_40, idiag_k100_30, idiag_k100_20, idiag_k100_10, &
idiag_k80_80, idiag_k80_60, idiag_k80_50, idiag_k80_40, &
idiag_k80_30, idiag_k80_20, idiag_k80_10, idiag_k60_60, &
idiag_k60_50, idiag_k60_40, idiag_k60_30, idiag_k60_20, &
idiag_k60_10, idiag_k50_50, idiag_k50_40, idiag_k50_30, &
idiag_k50_20, idiag_k50_10, idiag_k40_40, idiag_k40_30, &
idiag_k40_20, idiag_k40_10, idiag_k30_30, idiag_k30_20, &
idiag_k30_10, idiag_k20_20, idiag_k20_10, idiag_k10_10
integer :: rbin = 0
!
namelist /particles_coag_run_pars/ &
cdtpcoag, lcoag_simultaneous, lshear_in_vp, lconstant_kernel_test, &
kernel_cst, llinear_kernel_test, kernel_lin, lproduct_kernel_test, &
kernel_pro, lnoselfcollision, lgravitational_cross_section, &
GNewton, deltav_grav_floor, critical_mass_ratio_sticking, &
minimum_particle_mass, minimum_particle_radius, lzsomdullemond, &
lconstant_deltav, lmaxwell_deltav, deltav, maxwell_param, &
ldroplet_coagulation, droplet_coagulation_model, lcollision_output, &
luser_random_number_wrapper, lrelabelling, rdifference, &
kernel_output, deltad, a0, rbin, &
radius_output, r1, r2, r3, r4, r5, r6, r7, r8, r_diff, &
sphericalKernel, normal_coagulation, tstart_droplet_coagulation, &
lremove_particle, lremove_particle2, &
lcollision_output_swapped
!
contains
!***********************************************************************
subroutine initialize_particles_coag(f)
!
! Perform any post-parameter-read initialization i.e. calculate derived
! parameters.
!
! 24-nov-10/anders: coded
! 15-nov-12/KWJ: modified
!
use SharedVariables, only: get_shared_variable
!
real, dimension (mx,my,mz,mfarray) :: f
!
! Fatal error if Particle_radius module not used.
!
if (.not.lparticles_radius) &
call fatal_error('initialize_particles_coag', &
'must use Particles_radius module for coagulation')
!
! Allocate neighbour array necessary for identifying collisions.
!
if (.not.allocated(kneighbour)) allocate(kneighbour(mpar_loc))
!
! Get the gravity constant from the Selfgravity module.
!
if (lselfgravity) then
call get_shared_variable('rhs_poisson_const',rhs_poisson_const)
GNewton=rhs_poisson_const/(4*pi)
endif
!
! Precalculate inverse of coagulation time-step parameter.
!
cdtpcoag1=1/cdtpcoag
!
! Short hand for any kernel test.
!
lkernel_test=lconstant_kernel_test.or.llinear_kernel_test.or. &
lproduct_kernel_test
!
! Calculate squared and inverse volume factor.
!
four_pi_rhopmat_over_three2=four_pi_rhopmat_over_three**2
three_over_four_pi_rhopmat=1/four_pi_rhopmat_over_three
!
! Define the minimum particle mass or radius that will be allowed. Small
! particles may form by erosion and fragmentation.
!
if (minimum_particle_mass/=0.0 .and. minimum_particle_radius/=0.0) then
call fatal_error('initialize_particles_coag', &
'not allowed to set both minimum_particle_mass and '// &
'minimum_particle_radius')
endif
if (minimum_particle_mass/=0.0) then
minimum_particle_radius = &
(minimum_particle_mass*three_over_four_pi_rhopmat)**(1.0/3.0)
elseif (minimum_particle_radius/=0.0) then
minimum_particle_mass = &
four_pi_rhopmat_over_three*minimum_particle_radius**3
endif
!
if (.not.(lcartesian_coords.and.(all(lequidist)))) call fatal_error( &
'initialize_particles_coagulation', 'Coagulation only '// &
'implemented for Cartesian equidistant grids.')
!
! If using the Zsom-Dullemond Monte Carlo method
!
if(lzsomdullemond.and.(.not.allocated(r_ik_mat))) then
allocate(r_ik_mat(mpar_loc,mpar_loc))
allocate(r_i_tot(mpar_loc))
allocate(cum_func_first_i(mpar_loc))
allocate(cum_func_sec_ik(mpar_loc,mpar_loc))
end if
!
call keep_compiler_quiet(f)
!
endsubroutine initialize_particles_coag
!***********************************************************************
subroutine particles_coagulation_timestep(fp,ineargrid)
!
! Time-step contribution from particle coagulation.
!
! 30-nov-10/anders: coded
! 15-nov-12/KWJ: modified
!
real, dimension (mpar_loc,mparray) :: fp
integer, dimension (mpar_loc,3) :: ineargrid
!
real, dimension (3) :: xpj, xpk, vpj, vpk
real :: deltavjk, dt1_coag_par, kernel
real :: npswarmj, npswarmk
integer :: j, k, l
!
! If using the Zsom and Dullemond Monte Carlo method
!
if(lzsomdullemond) then
call particles_coag_timestep_zd(fp)
return
end if
!
if (lfirst.and.ldt) then
!
! Create list of shepherd and neighbour particles for each grid cell in the
! current pencil.
!
call shepherd_neighbour_pencil(fp,ineargrid,kshepherd,kneighbour)
!
! Calculate coagulation time-step.
!
do l=l1,l2
!
! Start with shepherd particle.
!
k=kshepherd(l-nghost)
if (k>0) then
do while (k/=0)
if (lparticles_number) then
npswarmk=fp(k,inpswarm)
else
npswarmk=rhop_swarm/(four_pi_rhopmat_over_three*fp(k,iap)**3)
endif
dt1_coag_par=0.0
j=kshepherd(l-nghost)
!
! Move through neighbours.
!
do while (.true.)
if (lparticles_number) then
npswarmj=fp(j,inpswarm)
else
npswarmj=rhop_swarm/(four_pi_rhopmat_over_three*fp(j,iap)**3)
endif
!
! Calculate the relative speed of particles j and k.
!
xpk=fp(k,ixp:izp)
vpk=fp(k,ivpx:ivpz)
if (lshear .and. lshear_in_vp) vpk(2)=vpk(2)-qshear*Omega*xpk(1)
xpj=fp(j,ixp:izp)
vpj=fp(j,ivpx:ivpz)
if (lshear .and. lshear_in_vp) vpj(2)=vpj(2)-qshear*Omega*xpj(1)
!
! Special treatment for kernel tests.
!
if (lkernel_test) then
if (lconstant_kernel_test) then
kernel=kernel_cst
elseif (llinear_kernel_test) then
kernel=kernel_lin* &
four_pi_rhopmat_over_three*(fp(j,iap)**3+fp(k,iap)**3)
elseif (lproduct_kernel_test) then
kernel=kernel_pro* &
four_pi_rhopmat_over_three2*fp(j,iap)**3*fp(k,iap)**3
endif
if (j==k .and. lnoselfcollision) kernel=0.0
dt1_coag_par=dt1_coag_par+kernel*min(npswarmj,npswarmk)
else
!
! Time-step for physical kernel.
!
if (sum((vpk-vpj)*(xpk-xpj))<0.0) then
deltavjk=sqrt(sum((vpk-vpj)**2))
dt1_coag_par=dt1_coag_par+ &
pi*(fp(k,iap)+fp(k,iap))**2*deltavjk* &
min(npswarmj,npswarmk)
endif
!
endif
j=kneighbour(j)
if (j==0) exit
enddo
!
! Put particle's time-step into inverse time-step array.
!
dt1_max(l-nghost)=max(dt1_max(l-nghost),dt1_coag_par*cdtpcoag1)
!
! Move to next particle in the grid cell.
!
k=kneighbour(k)
!
enddo
endif
enddo
endif
!
endsubroutine particles_coagulation_timestep
!***********************************************************************
subroutine particles_coagulation_pencils(fp,ineargrid)
!
! Calculate outcome of superparticle collisions by comparing the collision
! time-scale to the time-step. A random number is used to determine
! whether two superparticles collide in this time-step.
!
! Collisions lead to coagulation, shattering, erosion, or bouncing.
!
! 24-nov-10/anders: coded
! 15-nov-12/KWJ: modified
!
use Diagnostics
use General, only: random_number_wrapper
!
real, dimension (mpar_loc,mparray) :: fp
integer, dimension (mpar_loc,3) :: ineargrid
!
real, dimension (3) :: xpj, xpk, vpj, vpk
real :: lambda_mfp1, deltavjk, tau_coll1, prob, rran, kernel
real :: npswarmj, npswarmk
integer :: l, j, k, ncoll, ncoll_par, npart_par
!
intent (in) :: ineargrid
intent (inout) :: fp
!
real, dimension (10,10) :: kernel_array
real, dimension (10) :: radius_all
real, dimension (10) :: radius_ratio
real :: rmin, rmax
integer :: ibin, ik, ij, ikernel, row
integer, parameter :: max_rows = 10
real, parameter :: radius_diff=5.e-6
real, dimension (rbin) :: adplus, adminus, ad
character(len=50) :: itn,filename
real :: r100p, r100m, r80p, r80m, r60p, r60m, r50p, r50m, r40p, r40m, r30p, r30m, &
r20p, r20m, r10p, r10m
integer :: k100_100, k100_80, k100_60, k100_50, k100_40, k100_30, k100_20, k100_10, &
k80_80, k80_60, k80_50, k80_40, k80_30, k80_20, k80_10, &
k60_60, k60_50, k60_40, k60_30, k60_20, k60_10, &
k50_50, k50_40, k50_30, k50_20, k50_10, &
k40_40, k40_30, k40_20, k40_10, &
k30_30, k30_20, k30_10, &
k20_20, k20_10, &
k10_10
!
! a flag to start collision on the fly
if (t >= tstart_droplet_coagulation) ldroplet_coagulation=.true.
!
! If using the Zsom & Dullemond Monte Carlo method (KWJ)
!
if(lzsomdullemond) then
call particles_coag_outcome_zd(fp)
return
end if
!
! Reset collision counter.
!
ncoll=0
!
! Pencil loop.
!
do imn=1,ny*nz
!
! Create list of shepherd and neighbour particles for each grid cell in the
! current pencil.
!
call shepherd_neighbour_pencil(fp,ineargrid,kshepherd,kneighbour)
!
do l=l1,l2
!
! Start with shepherd particle.
!
k=kshepherd(l-nghost)
if (k>0) then
do while (k/=0)
if (lcoag_simultaneous) then
j=k
else
j=kshepherd(l-nghost)
endif
npart_par=0
ncoll_par=0
!17-06-30:Xiang-Yu coded.
k100_100=0; k100_80=0; k100_60=0; k100_50=0; k100_40=0; k100_30=0; k100_20=0; k100_10=0; &
k80_80=0; k80_60=0; k80_50=0; k80_40=0; k80_30=0; k80_20=0; k80_10=0; &
k60_60=0; k60_50=0; k60_40=0; k60_30=0; k60_20=0; k60_10=0; &
k50_50=0; k50_40=0; k50_30=0; k50_20=0; k50_10=0; &
k40_40=0; k40_30=0; k40_20=0; k40_10=0; &
k30_30=0; k30_20=0; k30_10=0; &
k20_20=0; k20_10=0; &
k10_10=0
!Xiang-Yu
!
if (lparticles_number) then
npswarmk=fp(k,inpswarm)
else
npswarmk=rhop_swarm/(four_pi_rhopmat_over_three*fp(k,iap)**3)
endif
!
! Move through neighbours, excluding particles tagged with negative radius
! (they have coagulated with a sink particle and will be removed later).
!
do while (fp(k,iap)>=0.0)
if (lcoag_simultaneous) then
j=kneighbour(j)
if (j==0) exit
!
! Do not attempt to collide with particles tagged for removal.
!
if (lparticles_number) then
do while (fp(j,iap)<0.0)
j=kneighbour(j)
if (j==0) exit
enddo
if (j==0) exit
endif
endif
!
! Particle number density either from f array or from rhop_swarm.
!
if (lparticles_number) then
npswarmj=fp(j,inpswarm)
else
npswarmj=rhop_swarm/(four_pi_rhopmat_over_three*fp(j,iap)**3)
endif
!
! Calculate the relative speed of particles j and k.
!
xpk=fp(k,ixp:izp)
vpk=fp(k,ivpx:ivpz)
if (lshear .and. lshear_in_vp) vpk(2)=vpk(2)-qshear*Omega*xpk(1)
xpj=fp(j,ixp:izp)
vpj=fp(j,ivpx:ivpz)
if (lshear .and. lshear_in_vp) vpj(2)=vpj(2)-qshear*Omega*xpj(1)
!
! Only consider collisions between particles approaching each other.
!
if (ldroplet_coagulation .or. &
(sum((vpk-vpj)*(xpk-xpj))<0.0) .or.&
lkernel_test) then
!
! Relative particle speed.
!
if (sphericalKernel) then
deltavjk=2*sqrt(sum(((vpk-vpj)*(xpk-xpj))**2))/sqrt(sum((xpk-xpj)**2))
else
deltavjk=sqrt(sum((vpk-vpj)**2))
endif
!
! Check for consistensy
!
if (droplet_coagulation_model=='standard') then
if (lcoag_simultaneous) then
call fatal_error('particles_coagulation_pencils',&
'One can not set droplet_coagulation_model=standard for lcoag_simultaneous=T')
endif
else if (droplet_coagulation_model=='shima') then
if (.not. lcoag_simultaneous) then
call fatal_error('particles_coagulation_pencils',&
'One can not set droplet_coagulation_model=shima for lcoag_simultaneous=F')
endif
endif
!
! The time-scale for collisions between a representative particle from
! superparticle k and the particle swarm in superparticle j is
!
! tau_coll = 1/(n_j*sigma_jk*deltav_jk)
!
! where n_j is the number density of swarm j and sigma_jk is the collisional
! cross section [=pi*(a_j+a_k)^2].
!
! The collision time-scale above is actually the *mass collision time-scale*,
! i.e. the time for a particle to collide with its own mass. Since super-
! particles always represent the same mass, that means that
!
! - The interaction time-scale for a particle k to collide with a swarm
! of larger particles j is simply the collision time-scale (when one
! small particle has collid with a large, then they all have collided).
!
! - The interaction time-scale for a particle k to collide with a swarm
! of lighter particles is the collision time-scale times the mass ratio
! m_k/m_j, so that tau_coll = 1/(n_k*sigma_jk*deltav_jk).
!
! We use lcoag_simultaneous to specify that two colliding superparticles
! change internal properties at the same time. The default is that the
! swarms evolve separately.
!
if (lkernel_test) then
if (lconstant_kernel_test) then
kernel=kernel_cst
elseif (llinear_kernel_test) then
kernel=kernel_lin* &
four_pi_rhopmat_over_three*(fp(j,iap)**3+fp(k,iap)**3)
elseif (lproduct_kernel_test) then
kernel=kernel_pro* &
four_pi_rhopmat_over_three2*fp(j,iap)**3*fp(k,iap)**3
endif
if (j==k .and. lnoselfcollision) kernel=0.0
if (fp(k,iap)<fp(j,iap)) then
tau_coll1=kernel*npswarmj
else
tau_coll1=kernel*npswarmk
endif
else
if (ldroplet_coagulation .and. .not. lcoag_simultaneous) then
tau_coll1=deltavjk*pi*(fp(k,iap)+fp(j,iap))**2*npswarmj
elseif (ldroplet_coagulation .and. &
lcoag_simultaneous .and. &
droplet_coagulation_model=='shima') then
tau_coll1=deltavjk*pi*(fp(k,iap)+fp(j,iap))**2* &
max(npswarmj,npswarmk)
elseif (normal_coagulation) then
tau_coll1=deltavjk*pi*(fp(k,iap)+fp(j,iap))**2* &
min(npswarmj,npswarmk)
else
tau_coll1=0
if (lgravitational_cross_section) then
tau_coll1=tau_coll1*(1.0+2*GNewton* &
four_pi_rhopmat_over_three* &
(fp(j,iap)**3+fp(k,iap)**3)/((fp(k,iap)+fp(k,iap))* &
(deltavjk+deltav_grav_floor)**2))
endif
endif
endif
!
if (tau_coll1/=0.0) then
!
! The probability for a collision in this time-step is dt/tau_coll.
!
prob=dt*tau_coll1
!NILS: I found that when the particles_coagulation module was turned on,
!NILS: the fluid turbulence was affected. By inspecting the energy spectrum
!NILS: it became apparent that it was primarily the small scales that were
!NILS: affected. This happened even though all back reactions from the
!NILS: particles to the fluid were turned off. After some intense debugging
!NILS: I finally found that the problem was due to a call to
!NILS: the random_number_wrapper routine. When I replaced this call with a call
!NILS: to the intinsic random_number subroutine, the energy spectrum was fine.
!NILS:
!NILS: I am not able to understand why the call to the random_number_wrapper
!NILS: routine should cause any problems, but I have reproduce the problem
!NILS: both on our linux cluster (with pfg90) and on my laptop (with gfortran).
if (luser_random_number_wrapper) then
call random_number_wrapper(rran)
else
call random_number(rran)
endif
if (rran<=prob) then
!
call coagulation_fragmentation(fp,j,k,deltavjk, &
xpj,xpk,vpj,vpk)
!
if (lparticles_number) then
npswarmk=fp(k,inpswarm)
else
npswarmk=rhop_swarm/(four_pi_rhopmat_over_three* &
fp(k,iap)**3)
endif
!
ncoll=ncoll+1
ncoll_par=ncoll_par+1
!
!17-06-21: Xiang-Yu coded: kernel of ri rj, diagnostics as time series
if (radius_output) then
!logorithmic binning
do ibin = 1,rbin
ad(ibin) = a0*deltad**ibin
enddo
do ibin = 1,rbin
if (abs(ad(ibin)-r1) == minval(abs(ad-r1))) then
r100p = ad(ibin+1)
r100m = ad(ibin-1)
endif
if (abs(ad(ibin)-r2) == minval(abs(ad-r2))) then
r80p = ad(ibin+1)
r80m = ad(ibin-1)
endif
if (abs(ad(ibin)-r3) == minval(abs(ad-r3))) then
r60p = ad(ibin+1)
r60m = ad(ibin-1)
endif
if (abs(ad(ibin)-r4) == minval(abs(ad-r4))) then
r50p = ad(ibin+1)
r50m = ad(ibin-1)
endif
if (abs(ad(ibin)-r5) == minval(abs(ad-r5))) then
r40p = ad(ibin+1)
r40m = ad(ibin-1)
endif
if (abs(ad(ibin)-r6) == minval(abs(ad-r6))) then
r30p = ad(ibin+1)
r30m = ad(ibin-1)
endif
if (abs(ad(ibin)-r7) == minval(abs(ad-r7))) then
r20p = ad(ibin+1)
r20m = ad(ibin-1)
endif
if (abs(ad(ibin)-r8) == minval(abs(ad-r8))) then
r10p = ad(ibin+1)
r10m = ad(ibin-1)
endif
enddo
if (fp(j,iap)>=fp(k,iap)) then
!r1=100 mum
if (fp(j,iap)>=r100m .and. fp(j,iap)<=r100p .and. fp(k,iap)>=r100m &
.and. fp(k,iap)<=r100p) k100_100 = k100_100 + 1
if (fp(j,iap)>=r100m .and. fp(j,iap)<=r100p .and. fp(k,iap)>=r80m &
.and. fp(k,iap)<=r80p) k100_80 = k100_80 + 1
if (fp(j,iap)>=r100m .and. fp(j,iap)<=r100p .and. fp(k,iap)>=r60m &
.and. fp(k,iap)<=r60p) k100_60 = k100_60 + 1
if (fp(j,iap)>=r100m .and. fp(j,iap)<=r100p .and. fp(k,iap)>=r50m &
.and. fp(k,iap)<=r50p) k100_50 = k100_50 + 1
if (fp(j,iap)>=r100m .and. fp(j,iap)<=r100p .and. fp(k,iap)>=r40m &
.and. fp(k,iap)<=r40p) k100_40 = k100_40 + 1
if (fp(j,iap)>=r100m .and. fp(j,iap)<=r100p .and. fp(k,iap)>=r30m &
.and. fp(k,iap)<=r30p) k100_30 = k100_30 + 1
if (fp(j,iap)>=r100m .and. fp(j,iap)<=r100p .and. fp(k,iap)>=r20m &
.and. fp(k,iap)<=r20p) k100_20 = k100_20 + 1
if (fp(j,iap)>=r100m .and. fp(j,iap)<=r100p .and. fp(k,iap)>=r10m &
.and. fp(k,iap)<=r10p) k100_10 = k100_10 + 1
!r2=80 mum
if (fp(j,iap)>=r80m .and. fp(j,iap)<=r80p .and. fp(k,iap)>=r80m .and. fp(k,iap)<=r80p) k80_80 = k80_80 + 1
if (fp(j,iap)>=r80m .and. fp(j,iap)<=r80p .and. fp(k,iap)>=r60m .and. fp(k,iap)<=r60p) k80_60 = k80_60 + 1
if (fp(j,iap)>=r80m .and. fp(j,iap)<=r80p .and. fp(k,iap)>=r50m .and. fp(k,iap)<=r50p) k80_50 = k80_50 + 1
if (fp(j,iap)>=r80m .and. fp(j,iap)<=r80p .and. fp(k,iap)>=r40m .and. fp(k,iap)<=r40p) k80_40 = k80_40 + 1
if (fp(j,iap)>=r80m .and. fp(j,iap)<=r80p .and. fp(k,iap)>=r30m .and. fp(k,iap)<=r30p) k80_30 = k80_30 + 1
if (fp(j,iap)>=r80m .and. fp(j,iap)<=r80p .and. fp(k,iap)>=r20m .and. fp(k,iap)<=r20p) k80_20 = k80_20 + 1
if (fp(j,iap)>=r80m .and. fp(j,iap)<=r80p .and. fp(k,iap)>=r10m .and. fp(k,iap)<=r10p) k80_10 = k80_10 + 1
!r3=60 mum
if (fp(j,iap)>=r60m .and. fp(j,iap)<=r60p .and. fp(k,iap)>=r60m .and. fp(k,iap)<=r60p) k60_60 = k60_60 + 1
if (fp(j,iap)>=r60m .and. fp(j,iap)<=r60p .and. fp(k,iap)>=r50m .and. fp(k,iap)<=r50p) k60_50 = k60_50 + 1
if (fp(j,iap)>=r60m .and. fp(j,iap)<=r60p .and. fp(k,iap)>=r40m .and. fp(k,iap)<=r40p) k60_40 = k60_40 + 1
if (fp(j,iap)>=r60m .and. fp(j,iap)<=r60p .and. fp(k,iap)>=r30m .and. fp(k,iap)<=r30p) k60_30 = k60_30 + 1
if (fp(j,iap)>=r60m .and. fp(j,iap)<=r60p .and. fp(k,iap)>=r20m .and. fp(k,iap)<=r20p) k60_20 = k60_20 + 1
if (fp(j,iap)>=r60m .and. fp(j,iap)<=r60p .and. fp(k,iap)>=r10m .and. fp(k,iap)<=r10p) k60_10 = k60_10 + 1
!r4=50 mum
if (fp(j,iap)>=r50m .and. fp(j,iap)<=r50p .and. fp(k,iap)>=r50m .and. fp(k,iap)<=r50p) k50_50 = k50_50 + 1
if (fp(j,iap)>=r50m .and. fp(j,iap)<=r50p .and. fp(k,iap)>=r40m .and. fp(k,iap)<=r40p) k50_40 = k50_40 + 1
if (fp(j,iap)>=r50m .and. fp(j,iap)<=r50p .and. fp(k,iap)>=r30m .and. fp(k,iap)<=r30p) k50_30 = k50_30 + 1
if (fp(j,iap)>=r50m .and. fp(j,iap)<=r50p .and. fp(k,iap)>=r20m .and. fp(k,iap)<=r20p) k50_20 = k50_20 + 1
if (fp(j,iap)>=r50m .and. fp(j,iap)<=r50p .and. fp(k,iap)>=r10m .and. fp(k,iap)<=r10p) k50_10 = k50_10 + 1
!r5=40 mum
if (fp(j,iap)>=r40m .and. fp(j,iap)<=r40p .and. fp(k,iap)>=r40m .and. fp(k,iap)<=r40p) k40_40 = k40_40 + 1
if (fp(j,iap)>=r40m .and. fp(j,iap)<=r40p .and. fp(k,iap)>=r30m .and. fp(k,iap)<=r30p) k40_30 = k40_30 + 1
if (fp(j,iap)>=r40m .and. fp(j,iap)<=r40p .and. fp(k,iap)>=r20m .and. fp(k,iap)<=r20p) k40_20 = k40_20 + 1
if (fp(j,iap)>=r40m .and. fp(j,iap)<=r40p .and. fp(k,iap)>=r10m .and. fp(k,iap)<=r10p) k40_10 = k40_10 + 1
!r6=30 mum
if (fp(j,iap)>=r30m .and. fp(j,iap)<=r30p .and. fp(k,iap)>=r30m .and. fp(k,iap)<=r30p) k30_30 = k30_30 + 1
if (fp(j,iap)>=r30m .and. fp(j,iap)<=r30p .and. fp(k,iap)>=r20m .and. fp(k,iap)<=r20p) k30_20 = k30_20 + 1
if (fp(j,iap)>=r30m .and. fp(j,iap)<=r30p .and. fp(k,iap)>=r10m .and. fp(k,iap)<=r10p) k30_10 = k30_10 + 1
!r7=20 mum
if (fp(j,iap)>=r20m .and. fp(j,iap)<=r20p .and. fp(k,iap)>=r20m .and. fp(k,iap)<=r20p) k20_20 = k20_20 + 1
if (fp(j,iap)>=r20m .and. fp(j,iap)<=r20p .and. fp(k,iap)>=r10m .and. fp(k,iap)<=r10p) k20_10 = k20_10 + 1
!r7=10 mum
if (fp(j,iap)>=r10m .and. fp(j,iap)<=r10p .and. fp(k,iap)>=r10m .and. fp(k,iap)<=r10p) k10_10 = k10_10 + 1
else
!r1=100 mum
if (fp(k,iap)>=r100m .and. fp(k,iap)<=r100p .and. fp(j,iap)>=r100m &
.and. fp(j,iap)<=r100p) k100_100 = k100_100 + 1
if (fp(k,iap)>=r100m .and. fp(k,iap)<=r100p .and. fp(j,iap)>=r80m &
.and. fp(j,iap)<=r80p) k100_80 = k100_80 + 1
if (fp(k,iap)>=r100m .and. fp(k,iap)<=r100p .and. fp(j,iap)>=r60m &
.and. fp(j,iap)<=r60p) k100_60 = k100_60 + 1
if (fp(k,iap)>=r100m .and. fp(k,iap)<=r100p .and. fp(j,iap)>=r50m &
.and. fp(j,iap)<=r50p) k100_50 = k100_50 + 1
if (fp(k,iap)>=r100m .and. fp(k,iap)<=r100p .and. fp(j,iap)>=r40m &
.and. fp(j,iap)<=r40p) k100_40 = k100_40 + 1
if (fp(k,iap)>=r100m .and. fp(k,iap)<=r100p .and. fp(j,iap)>=r30m &
.and. fp(j,iap)<=r30p) k100_30 = k100_30 + 1
if (fp(k,iap)>=r100m .and. fp(k,iap)<=r100p .and. fp(j,iap)>=r20m &
.and. fp(j,iap)<=r20p) k100_20 = k100_20 + 1
if (fp(k,iap)>=r100m .and. fp(k,iap)<=r100p .and. fp(j,iap)>=r10m &
.and. fp(j,iap)<=r10p) k100_10 = k100_10 + 1
!r2=80 mum
if (fp(k,iap)>=r80m .and. fp(k,iap)<=r80p .and. fp(j,iap)>=r80m .and. fp(j,iap)<=r80p) k80_80 = k80_80 + 1
if (fp(k,iap)>=r80m .and. fp(k,iap)<=r80p .and. fp(j,iap)>=r60m .and. fp(j,iap)<=r60p) k80_60 = k80_60 + 1
if (fp(k,iap)>=r80m .and. fp(k,iap)<=r80p .and. fp(j,iap)>=r50m .and. fp(j,iap)<=r50p) k80_50 = k80_50 + 1
if (fp(k,iap)>=r80m .and. fp(k,iap)<=r80p .and. fp(j,iap)>=r40m .and. fp(j,iap)<=r40p) k80_40 = k80_40 + 1
if (fp(k,iap)>=r80m .and. fp(k,iap)<=r80p .and. fp(j,iap)>=r30m .and. fp(j,iap)<=r30p) k80_30 = k80_30 + 1
if (fp(k,iap)>=r80m .and. fp(k,iap)<=r80p .and. fp(j,iap)>=r20m .and. fp(j,iap)<=r20p) k80_20 = k80_20 + 1
if (fp(k,iap)>=r80m .and. fp(k,iap)<=r80p .and. fp(j,iap)>=r10m .and. fp(j,iap)<=r10p) k80_10 = k80_10 + 1
!r3=60 mum
if (fp(k,iap)>=r60m .and. fp(k,iap)<=r60p .and. fp(j,iap)>=r60m .and. fp(j,iap)<=r60p) k60_60 = k60_60 + 1
if (fp(k,iap)>=r60m .and. fp(k,iap)<=r60p .and. fp(j,iap)>=r50m .and. fp(j,iap)<=r50p) k60_50 = k60_50 + 1
if (fp(k,iap)>=r60m .and. fp(k,iap)<=r60p .and. fp(j,iap)>=r40m .and. fp(j,iap)<=r40p) k60_40 = k60_40 + 1
if (fp(k,iap)>=r60m .and. fp(k,iap)<=r60p .and. fp(j,iap)>=r30m .and. fp(j,iap)<=r30p) k60_30 = k60_30 + 1
if (fp(k,iap)>=r60m .and. fp(k,iap)<=r60p .and. fp(j,iap)>=r20m .and. fp(j,iap)<=r20p) k60_20 = k60_20 + 1
if (fp(k,iap)>=r60m .and. fp(k,iap)<=r60p .and. fp(j,iap)>=r10m .and. fp(j,iap)<=r10p) k60_10 = k60_10 + 1
!r4=50 mum
if (fp(k,iap)>=r50m .and. fp(k,iap)<=r50p .and. fp(j,iap)>=r50m .and. fp(j,iap)<=r50p) k50_50 = k50_50 + 1
if (fp(k,iap)>=r50m .and. fp(k,iap)<=r50p .and. fp(j,iap)>=r40m .and. fp(j,iap)<=r40p) k50_40 = k50_40 + 1
if (fp(k,iap)>=r50m .and. fp(k,iap)<=r50p .and. fp(j,iap)>=r30m .and. fp(j,iap)<=r30p) k50_30 = k50_30 + 1
if (fp(k,iap)>=r50m .and. fp(k,iap)<=r50p .and. fp(j,iap)>=r20m .and. fp(j,iap)<=r20p) k50_20 = k50_20 + 1
if (fp(k,iap)>=r50m .and. fp(k,iap)<=r50p .and. fp(j,iap)>=r10m .and. fp(j,iap)<=r10p) k50_10 = k50_10 + 1
!r5=40 mum
if (fp(k,iap)>=r40m .and. fp(k,iap)<=r40p .and. fp(j,iap)>=r40m .and. fp(j,iap)<=r40p) k40_40 = k40_40 + 1
if (fp(k,iap)>=r40m .and. fp(k,iap)<=r40p .and. fp(j,iap)>=r30m .and. fp(j,iap)<=r30p) k40_30 = k40_30 + 1
if (fp(k,iap)>=r40m .and. fp(k,iap)<=r40p .and. fp(j,iap)>=r20m .and. fp(j,iap)<=r20p) k40_20 = k40_20 + 1
if (fp(k,iap)>=r40m .and. fp(k,iap)<=r40p .and. fp(j,iap)>=r10m .and. fp(j,iap)<=r10p) k40_10 = k40_10 + 1
!r6=30 mum
if (fp(k,iap)>=r30m .and. fp(k,iap)<=r30p .and. fp(j,iap)>=r30m .and. fp(j,iap)<=r30p) k30_30 = k30_30 + 1
if (fp(k,iap)>=r30m .and. fp(k,iap)<=r30p .and. fp(j,iap)>=r20m .and. fp(j,iap)<=r20p) k30_20 = k30_20 + 1
if (fp(k,iap)>=r30m .and. fp(k,iap)<=r30p .and. fp(j,iap)>=r10m .and. fp(j,iap)<=r10p) k30_10 = k30_10 + 1
!r7=20 mum
if (fp(k,iap)>=r20m .and. fp(k,iap)<=r20p .and. fp(j,iap)>=r20m .and. fp(j,iap)<=r20p) k20_20 = k20_20 + 1
if (fp(k,iap)>=r20m .and. fp(k,iap)<=r20p .and. fp(j,iap)>=r10m .and. fp(j,iap)<=r10p) k20_10 = k20_10 + 1
!r7=10 mum
if (fp(k,iap)>=r10m .and. fp(k,iap)<=r10p .and. fp(j,iap)>=r10m .and. fp(j,iap)<=r10p) k10_10 = k10_10 + 1
endif
endif
!17-06-21: Xiang-Yu
!17-06-18: Xiang-Yu coded: kernel of ri rj, diagnostics as 2-D matrix
if (kernel_output) then
! read radius and radius ratio
open(unit=11,file="radius.txt")
do row = 1, max_rows
read(11,*) radius_all(row)
enddo
close(unit=11)
rmin = minval(radius_all)
rmax = maxval(radius_all)
!open(unit=12,file="ratio.txt")
!do row = 1, max_rows
! read(12,*) radius_ratio(row)
!enddo
!close(unit=12)
!logorithmic binning
do ibin = 1,rbin
adminus(ibin) = a0*deltad**(ibin-1)
adplus(ibin) = a0*deltad**ibin
ad(ibin) = (adminus(ibin)+adplus(ibin))*.5
enddo
! search for collector and collected particles and bin them in the kernel
if (fp(k,iap) >= rmin .and. fp(k,iap) <= rmax) then
! ikernel=0
do ibin=1,rbin
if (abs(fp(k,iap)-radius_all(ibin)) == minval(abs(fp(k,iap)-radius_all))) then
ik=ibin
endif
enddo
do ibin=1,rbin
if (abs(fp(j,iap)-radius_all(ibin)) == minval(abs(fp(j,iap)-radius_all))) then
ij=ibin
endif
enddo
! ikernel=ikernel+1
!kernel_array(ik,ij)=(kernel_array(ik,ij)+tau_coll1)/ikernel
kernel_array(ik,ij) = tau_coll1
! output the kernel
write(itn,'(I5)') it ! convert integer to char
filename = adjustl(trim(adjustr(itn)))
open(13, file=trim(directory_dist)//'/'//trim(filename)//'.dat')
write(13,"(10e11.3)") kernel_array
close(13)
! print*,'kernel_array=',kernel_array
! print*,'tau_coll1=',tau_coll1
endif
endif
!17-06-18:XY
endif
endif
endif
!
! Move to next particle neighbour.
!
if (.not.lcoag_simultaneous) then
j=kneighbour(j)
if (j==0) exit
endif
!
enddo
!
! Collision diagnostics. Since this subroutine is called in the last sub-
! time-step, we can not use ldiagnos. Therefore we calculate collision
! diagnostics in the preceding time-step. This has the side effect that
! collision diagnostics are not particle normalized in it==1 or if it1==1.
!
if (it==1 .or. mod(it,it1)==0) then
if (idiag_ncoagpm/=0) &
call sum_par_name((/real(ncoll_par)/),idiag_ncoagpm)
if (idiag_ncoagpartpm/=0) &
call sum_par_name((/real(npart_par)/),idiag_ncoagpartpm)
if (idiag_dt1_coag_par/=0) &
call sum_par_name((/real(ncoll_par)/),idiag_dt1_coag_par)
if (idiag_k100_100/=0) &
call sum_par_name((/real(k100_100)/),idiag_k100_100)
if (idiag_k100_80/=0) &
call sum_par_name((/real(k100_80)/),idiag_k100_80)
if (idiag_k100_60/=0) &
call sum_par_name((/real(k100_60)/),idiag_k100_60)
if (idiag_k100_50/=0) &
call sum_par_name((/real(k100_50)/),idiag_k100_50)
if (idiag_k100_40/=0) &
call sum_par_name((/real(k100_40)/),idiag_k100_40)
if (idiag_k100_30/=0) &
call sum_par_name((/real(k100_30)/),idiag_k100_30)
if (idiag_k100_20/=0) &
call sum_par_name((/real(k100_20)/),idiag_k100_20)
if (idiag_k100_10/=0) &
call sum_par_name((/real(k100_10)/),idiag_k100_10)
if (idiag_k80_80/=0) call sum_par_name((/real(k80_80)/),idiag_k80_80)
if (idiag_k80_60/=0) call sum_par_name((/real(k80_60)/),idiag_k80_60)
if (idiag_k80_50/=0) call sum_par_name((/real(k80_50)/),idiag_k80_50)
if (idiag_k80_40/=0) call sum_par_name((/real(k80_40)/),idiag_k80_40)
if (idiag_k80_30/=0) call sum_par_name((/real(k80_30)/),idiag_k80_30)
if (idiag_k80_20/=0) call sum_par_name((/real(k80_20)/),idiag_k80_20)
if (idiag_k80_10/=0) call sum_par_name((/real(k80_10)/),idiag_k80_10)
if (idiag_k60_60/=0) call sum_par_name((/real(k60_60)/),idiag_k60_60)
if (idiag_k60_50/=0) call sum_par_name((/real(k60_50)/),idiag_k60_50)
if (idiag_k60_40/=0) call sum_par_name((/real(k60_40)/),idiag_k60_40)
if (idiag_k60_30/=0) call sum_par_name((/real(k60_30)/),idiag_k60_30)
if (idiag_k60_20/=0) call sum_par_name((/real(k60_20)/),idiag_k60_20)
if (idiag_k60_10/=0) call sum_par_name((/real(k60_10)/),idiag_k60_10)
if (idiag_k50_50/=0) call sum_par_name((/real(k50_50)/),idiag_k50_50)
if (idiag_k50_40/=0) call sum_par_name((/real(k50_40)/),idiag_k50_40)
if (idiag_k50_30/=0) call sum_par_name((/real(k50_30)/),idiag_k50_30)
if (idiag_k50_20/=0) call sum_par_name((/real(k50_20)/),idiag_k50_20)
if (idiag_k50_10/=0) call sum_par_name((/real(k50_10)/),idiag_k50_10)
if (idiag_k40_40/=0) call sum_par_name((/real(k40_40)/),idiag_k40_40)
if (idiag_k40_30/=0) call sum_par_name((/real(k40_30)/),idiag_k40_30)
if (idiag_k40_20/=0) call sum_par_name((/real(k40_20)/),idiag_k40_20)
if (idiag_k40_10/=0) call sum_par_name((/real(k40_10)/),idiag_k40_10)
if (idiag_k30_30/=0) call sum_par_name((/real(k30_30)/),idiag_k30_30)
if (idiag_k30_20/=0) call sum_par_name((/real(k30_20)/),idiag_k30_20)
if (idiag_k30_10/=0) call sum_par_name((/real(k30_10)/),idiag_k30_10)
if (idiag_k20_20/=0) call sum_par_name((/real(k20_20)/),idiag_k20_20)
if (idiag_k20_10/=0) call sum_par_name((/real(k20_10)/),idiag_k20_10)
if (idiag_k10_10/=0) call sum_par_name((/real(k10_10)/),idiag_k10_10)
endif
!
! Move to next particle in the grid cell.
!
k=kneighbour(k)
!
enddo
endif
enddo
enddo
!
! We need to register the diagnostic type, even if there are no particles
! at the local processor. In the same spirit we calculate diagnostics even
! for it==1 on processors that do have particles (above). Otherwise a
! processor starting with N>0 particles, but arriving at mod(it,it1)==0 with
! zero particles, will not register collision diagnostics properly.
!
if (it==1) then
if (npar_loc==0) then
if (idiag_ncoagpm/=0) &
call sum_par_name(fp(1:npar_loc,ixp),idiag_ncoagpm)
if (idiag_ncoagpartpm/=0) &
call sum_par_name(fp(1:npar_loc,ixp),idiag_ncoagpartpm)
if (idiag_dt1_coag_par/=0) &
call sum_par_name(fp(1:npar_loc,ixp),idiag_dt1_coag_par)
endif
endif
!
! Remove the particles that have been tagged for removal. We have to remove
! them after the coagulation step, as we will otherwise confuse the
! shepherd-neighbour algorithm.
!
k=1
do while (.true.)
if (fp(k,iap)<0.0) then
fp(k,iap)=-fp(k,iap)
call remove_particle(fp,ipar,k)
else
k=k+1
if (k>npar_loc) exit
endif
enddo
!
endsubroutine particles_coagulation_pencils
!***********************************************************************
subroutine particles_coagulation_blocks(fp,ineargrid)
!
! Calculate outcome of superparticle collisions by comparing the collision
! time-scale to the time-step. A random number is used to determine
! whether two superparticles collide in this time-step.
!
! Collisions lead to coagulation, shattering, erosion, or bouncing.
!
! 24-nov-10/anders: coded
!
real, dimension (mpar_loc,mparray) :: fp
integer, dimension (mpar_loc,3) :: ineargrid
!
call fatal_error('particles_coagulation_blocks','not implemented yet')
!
call keep_compiler_quiet(fp)
call keep_compiler_quiet(ineargrid)
!
endsubroutine particles_coagulation_blocks
!***********************************************************************
subroutine coagulation_fragmentation(fp,j,k,deltavjk,xpj,xpk,vpj,vpk)
!
! Change the particle size to the new size.
! If droplet_coagulation_model does not equal 'shima' the total mass in the
! particle swarm is kept constant.
!
! A representative particle k colliding with a swarm of larger particles j
! simply obtains the new mass m_k -> m_k + m_j.
!
! A representative particle k colliding with a swarm of lighter particles j
! obtains the new mass m_k -> 2*m_k (this is true only when
! ldroplet_coagulation=F.)
!
! 24-nov-10/anders: coded
! 25-jan-16/Xiang-yu and nils: added the model of Shima et al. (2009)
!
real, dimension (mpar_loc,mparray) :: fp
integer :: j, k, swarm_index1, swarm_index2, iswap, ipar_j_
real :: deltavjk
real, dimension (3) :: xpj, xpk, vpj, vpk, vpnew
!
real, dimension (3) :: vvcm, vvkcm, vvkcmnew, vvkcm_normal, vvkcm_parall
real, dimension (3) :: nvec
real :: mpsma, mpbig, npsma, npbig, npnew, mpnew, apnew
real :: rhopsma, rhopbig, apsma, apbig
real :: coeff_restitution, deltav_recoil, escape_speed
real :: apj,mpj,npj,apk,mpk,npk
logical :: lswap
!
!
if (lparticles_number) then
!
! Physical particles in the two swarms collide simultaneously.
!
if (lcoag_simultaneous) then
!
! Define mpsma, npsma, rhopsma for the superparticle with lighter particles.
! Define mpbig, npbig, rhopbig for the superparticle with bigger particles.
!
if (fp(j,iap)<fp(k,iap)) then
apsma = fp(j,iap)
mpsma = four_pi_rhopmat_over_three*fp(j,iap)**3
npsma = fp(j,inpswarm)
apbig = fp(k,iap)
mpbig = four_pi_rhopmat_over_three*fp(k,iap)**3
npbig = fp(k,inpswarm)
else
apsma = fp(k,iap)
mpsma = four_pi_rhopmat_over_three*fp(k,iap)**3
npsma = fp(k,inpswarm)
apbig = fp(j,iap)
mpbig = four_pi_rhopmat_over_three*fp(j,iap)**3
npbig = fp(j,inpswarm)
endif
rhopsma=mpsma*npsma
rhopbig=mpbig*npbig
!
! If we are working on droploet coagulation all colliding droplets are
! assumed to coalesce. I.e. the mass of the new representative droplet equals
! the sum of the two initial droplet masses.
!
if (ldroplet_coagulation) then
if (droplet_coagulation_model=='standard') then
mpnew=mpsma+mpbig
apnew=(mpnew*three_over_four_pi_rhopmat)**(1.0/3.0)
npnew=(rhopsma+rhopbig)/(2.*mpnew)
vpnew=(&
fp(j,ivpx:ivpz)*four_pi_rhopmat_over_three*fp(j,iap)**3+&
fp(k,ivpx:ivpz)*four_pi_rhopmat_over_three*fp(k,iap)**3)/mpnew
fp(j,iap)=apnew
fp(k,iap)=apnew
fp(j,inpswarm)=npnew
fp(k,inpswarm)=npnew
fp(j,ivpx:ivpz)=vpnew
fp(k,ivpx:ivpz)=vpnew
else if (droplet_coagulation_model=='shima') then
mpj = four_pi_rhopmat_over_three*fp(j,iap)**3
npj = fp(j,inpswarm)
mpk = four_pi_rhopmat_over_three*fp(k,iap)**3
npk = fp(k,inpswarm)
!
! Identify the swarm with the highest particle number density
!
if (npk<npj) then
swarm_index1=k
swarm_index2=j
lswap=mpk<mpj
else
swarm_index1=j
swarm_index2=k
lswap=mpk>mpj
endif
if (lcollision_output) then
open(99,POSITION='append', &
FILE=trim(directory_dist)//'/collisions.dat')
write(99,"(f18.6,2i10,2f12.8,1p,2e11.3)") t,ipar(j),ipar(k),fp(j,iap),fp(k,iap),fp(j,inpswarm),fp(k,inpswarm)
close(99)
endif
!
! Check if we have the special case where the particle number
! densities are equal. This will typically be the case in the beginning of
! a simulation.
!
if (npk == npj) then
! 13-June-18/Xiang-Yu: coded
if (npk*dx*dy*dz<1.0 .and. lremove_particle) then
if (fp(swarm_index1,iap)>=fp(swarm_index2,iap)) then
fp(swarm_index1,ivpx:ivpz)=(rhopbig*fp(swarm_index1,ivpx:ivpz) + &
rhopsma*fp(swarm_index2,ivpx:ivpz))/(rhopsma+rhopbig)
! fp(swarm_index1,inpswarm)=1.0/(dx*dy*dz)
fp(swarm_index1,iap)=((rhopsma+rhopbig)/fp(swarm_index1,inpswarm)/ &
four_pi_rhopmat_over_three)**(1.0/3.0)
fp(swarm_index2,iap)=-fp(swarm_index2,iap)
else
fp(swarm_index2,ivpx:ivpz)=(rhopbig*fp(swarm_index2,ivpx:ivpz) + &
rhopsma*fp(swarm_index1,ivpx:ivpz))/(rhopsma+rhopbig)
! fp(swarm_index2,inpswarm)=1.0/(dx*dy*dz)
fp(swarm_index2,iap)=((rhopsma+rhopbig)/fp(swarm_index2,inpswarm)/ &
four_pi_rhopmat_over_three)**(1.0/3.0)
fp(swarm_index1,iap)=-fp(swarm_index1,iap)
endif
else
! 13-June-18/Xiang-Yu.
fp(swarm_index1,inpswarm)=fp(swarm_index1,inpswarm)/2.
fp(swarm_index2,inpswarm)=fp(swarm_index2,inpswarm)/2.
mpnew=mpj+mpk
fp(swarm_index1,iap)&
=(mpnew*three_over_four_pi_rhopmat)**(1.0/3.0)
fp(swarm_index2,iap)&
=(mpnew*three_over_four_pi_rhopmat)**(1.0/3.0)
vpnew=(fp(j,ivpx:ivpz)*mpj+fp(k,ivpx:ivpz)*mpk)/mpnew
fp(swarm_index1,ivpx:ivpz)=vpnew
fp(swarm_index2,ivpx:ivpz)=vpnew
endif
else
!
! Set particle radius. The radius of the swarm with the highes particle
! number density does not change since its mass doesn't change.
!
mpnew=mpj+mpk
fp(swarm_index1,iap)&
=(mpnew*three_over_four_pi_rhopmat)**(1.0/3.0)
!
! Set particle number densities. Only the the swarm with the highest
! particle number density, changes its particle number density.
!
fp(swarm_index2,inpswarm)=abs(npj-npk)
!
! Set particle velocities. Only the swarm with the lowest
! particle number density, changes its particle velcity.
!
vpnew=(fp(j,ivpx:ivpz)*mpj+fp(k,ivpx:ivpz)*mpk)/mpnew
fp(swarm_index1,ivpx:ivpz)=vpnew
endif
!
! 10-July-18/Xiang-Yu coded: remove superparticle containing less than one droplet.
if (npnew*dx*dy*dz<1.0 .and. lremove_particle2) then
if (fp(j,iap)<fp(k,iap)) then
fp(k,ivpx:ivpz)=(rhopsma*fp(j,ivpx:ivpz) + &
rhopbig*fp(k,ivpx:ivpz))/(rhopsma+rhopbig)
else
fp(k,ivpx:ivpz)=(rhopbig*fp(j,ivpx:ivpz) + &
rhopsma*fp(k,ivpx:ivpz))/(rhopsma+rhopbig)
endif
fp(j,iap)=-fp(j,iap)
fp(k,iap)=((rhopsma+rhopbig)*dx*dy*dz/ &
four_pi_rhopmat_over_three)**(1.0/3.0)
fp(k,inpswarm)=1/(dx*dy*dz)
endif
! 10-July-18/Xiang-Yu.
!
! swap names ipar to follow lucky droplets
!
if (lrelabelling) then
!if (lswap.and.(fp(j,iap)/=fp(k,iap))) then
if (lswap.and.(abs(fp(j,iap)-fp(k,iap)))>rdifference) then
ipar_j_=ipar(k)
ipar(k)=ipar(j)
ipar(j)=ipar_j_
iswap=1
else
iswap=0
endif
endif
!
if (lcollision_output_swapped) then
open(99,POSITION='append', &
FILE=trim(directory_dist)//'/collisions_swapped.dat')
write(99,"(f18.6,2i10,2f12.8,1p,2e11.3)") &
t,ipar(j),ipar(k),fp(j,iap),fp(k,iap), &
fp(j,inpswarm),fp(k,inpswarm),iswap
close(99)
endif
!
else
call fatal_error('','No such droplet_coagulation_model')
endif
else
!
! Calculate the coefficient of restitution. We base in here on the
! experimental fit of Higa et al. (1986) to ice at 100 K.
!
coeff_restitution=(deltavjk/1.8)**(-alog10(deltavjk/1.8))
!
! Particles stick:
!
! a) when the mass ratio between large and small particles is
! sufficiently high.
!
! b) when the recoil velocity is less than the escape speed
!
nvec=(xpj-xpk)
nvec=nvec/sqrt(sum(nvec**2))
vvcm=0.5*(vpj+vpk)
vvkcm=vpk-vvcm
vvkcm_normal=nvec*(sum(vvkcm*nvec))
vvkcm_parall=vvkcm-vvkcm_normal
deltav_recoil=sqrt(sum((vvkcm_parall - &
coeff_restitution*vvkcm_normal)**2))
escape_speed =sqrt(2*GNewton*mpbig/apbig)
!
if (mpbig/mpsma>critical_mass_ratio_sticking .or. &
deltav_recoil<escape_speed .or. &
fp(j,inpswarm)==1/(dx*dy*dz) .or. &
fp(k,inpswarm)==1/(dx*dy*dz)) then
mpnew=mpbig+rhopsma/npbig
else
if (mpsma>2*minimum_particle_mass) then
mpnew=0.5*mpsma
else
mpnew=minimum_particle_mass
endif
endif
!
! The new radius is defined from the new particle mass, while the new
! particle number in each superparticle comes from total mass conservation.
!
apnew=(mpnew*three_over_four_pi_rhopmat)**(1.0/3.0)
npnew=0.5*(rhopsma+rhopbig)/mpnew
!
! Turn into sink particle if number of physical particles is less than one.
!
if (npnew*dx*dy*dz<1.0) then
if (fp(j,iap)<fp(k,iap)) then
fp(k,ivpx:ivpz)=(rhopsma*fp(j,ivpx:ivpz) + &
rhopbig*fp(k,ivpx:ivpz))/(rhopsma+rhopbig)
else
fp(k,ivpx:ivpz)=(rhopbig*fp(j,ivpx:ivpz) + &
rhopsma*fp(k,ivpx:ivpz))/(rhopsma+rhopbig)
endif
!
! Tag particle for removal by making the radius negative.
!
fp(j,iap)=-fp(j,iap)
fp(k,iap)=((rhopsma+rhopbig)*dx*dy*dz/ &
four_pi_rhopmat_over_three)**(1.0/3.0)
fp(k,inpswarm)=1/(dx*dy*dz)
else
fp(j,iap)=apnew
fp(k,iap)=apnew
fp(j,inpswarm)=npnew
fp(k,inpswarm)=npnew
endif
endif ! if (ldroplet_coagulation)
else
!
! If we are working on droploet coagulation all colliding droplets are
! assumed to coalesce. I.e. the mass of the new representative droplet equals
! the sum of the two initial droplet masses.
!
if (ldroplet_coagulation) then
apj = fp(j,iap)
mpj = four_pi_rhopmat_over_three*fp(j,iap)**3
npj = fp(j,inpswarm)
apk = fp(k,iap)
mpk = four_pi_rhopmat_over_three*fp(k,iap)**3
npk = fp(k,inpswarm)
if (droplet_coagulation_model=='standard') then
if (lcollision_output) then
open(99,POSITION='append', &
FILE=trim(directory_dist)//'/collisions.dat')
write(99,"(f18.6,2i10,2f12.8,1p,2e11.3)") t,ipar(j),ipar(k),fp(j,iap),fp(k,iap),fp(j,inpswarm),fp(k,inpswarm)
close(99)
endif
mpnew=mpj+mpk
apnew=(mpnew*three_over_four_pi_rhopmat)**(1.0/3.0)
! npnew=(mpj*npk+mpk*npk)/(2.*mpnew)
npnew=mpk*npk/mpnew
vpnew=(fp(j,ivpx:ivpz)*mpj+fp(k,ivpx:ivpz)*mpk)/mpnew
fp(k,iap)=apnew
fp(k,inpswarm)=npnew
fp(k,ivpx:ivpz)=vpnew
else if (droplet_coagulation_model=='shima') then
!
! Identify the swarm with the highest particle number density
!
if (npk<npj) then
swarm_index1=k
swarm_index2=j
else
swarm_index1=j
swarm_index2=k
endif
!
! Check if we have the special case where the particle number
! densities are equal. This will typically be the case in the beginning of
! a simulation.
!
if (npk == npj) then
fp(swarm_index1,inpswarm)=fp(swarm_index1,inpswarm)/2.
fp(swarm_index2,inpswarm)=fp(swarm_index2,inpswarm)/2.
mpnew=mpj+mpk
fp(swarm_index1,iap)&
=(mpnew*three_over_four_pi_rhopmat)**(1.0/3.0)
fp(swarm_index2,iap)&
=(mpnew*three_over_four_pi_rhopmat)**(1.0/3.0)
vpnew=(fp(j,ivpx:ivpz)*mpj+fp(k,ivpx:ivpz)*mpk)/mpnew
fp(swarm_index1,ivpx:ivpz)=vpnew
fp(swarm_index2,ivpx:ivpz)=vpnew
else
!
! Set particle radius. The radius of the swarm with the highes particle
! number density does not change since its mass doesn't change.
!
mpnew=mpj+mpk
fp(swarm_index1,iap)&
=(mpnew*three_over_four_pi_rhopmat)**(1.0/3.0)
!
! Set particle number densities. Only the the swarm with the highest
! particle number density, changes its particle number density.
!
fp(swarm_index2,inpswarm)=abs(npj-npk)
!
! Set particle velocities. Only the swarm with the lowest
! particle number density, changes its particle velcity.
!
vpnew=(fp(j,ivpx:ivpz)*mpj+fp(k,ivpx:ivpz)*mpk)/mpnew
fp(swarm_index1,ivpx:ivpz)=vpnew
endif
else
call fatal_error('','No such droplet_coagulation_model')
endif
else
!
! Physical particles in the two swarms collide asymmetrically.
!
if (fp(k,iap)<fp(j,iap)) then
fp(k,inpswarm)=fp(k,inpswarm)* &
(1/(1.0+(fp(j,iap)/fp(k,iap))**3))
fp(k,iap)=(fp(k,iap)**3+fp(j,iap)**3)**(1.0/3.0)
else
fp(k,inpswarm)=0.5*fp(k,inpswarm)
fp(k,iap)=2**(1.0/3.0)*fp(k,iap)
endif
endif
endif
else
!
! In case we only evolve the particle radius, number their number.
! We assume that the mass represented by each superparticle is the same.
!
if (lcoag_simultaneous) then
fp(j,iap)=2**(1.0/3.0)*max(fp(j,iap),fp(k,iap))
fp(k,iap)=fp(j,iap)
else
if (fp(k,iap)<fp(j,iap)) then
fp(k,iap)=(fp(k,iap)**3+fp(j,iap)**3)**(1.0/3.0)
else
fp(k,iap)=2**(1.0/3.0)*fp(k,iap)
endif
endif
endif
!
endsubroutine coagulation_fragmentation
!***********************************************************************
subroutine particles_coag_timestep_zd(fp)
!
! Calculate the next timestep with the same scheme as in Zsom & Dullemond
! (2008). Only use for kernel tests. Self-collision included.
!
! 15-nov-12/KWJ: coded
! 08-jan-13/KWJ: modified
!
use General, only: random_number_wrapper
! Get all variables
!NILS: Shouldn't npar be mpar_loc?
! real, dimension (npar,mparray) :: fp
real, dimension (mpar_loc,mparray) :: fp
!
real :: r, kernel, dt1_coag_par
integer :: j, l
!
! If first timestep calculate the r_ik-matrix and the total rates.
! Loop over all particle pairs.
!
! We also need to get the cumulative distributions to pick the particles.
! N.B. Non-normalized cumulative distributions for easier handle of coagulation
! outcome.
!
!
if (it==1) then
j=1
l=1
! Initialize all rates and cumulative functions
r_i_tot = 0.
r_total = 0.
cum_func_first_i = 0.
cum_func_sec_ik = 0.
! Loop over all particles
do while (.true.)
do while (.true.)
! Get kernel for jl-pair. Rep. particle j and phys. particle l
! Theoretical kernels with analytic solutions
if (lconstant_kernel_test) then
kernel=kernel_cst
elseif (llinear_kernel_test) then
kernel=kernel_lin* &
four_pi_rhopmat_over_three*(fp(j,iap)**3+fp(l,iap)**3)
elseif (lproduct_kernel_test) then
kernel=kernel_pro* &
four_pi_rhopmat_over_three2*fp(j,iap)**3*fp(l,iap)**3
! Physical kernels
elseif (lconstant_deltav) then ! constant relative velocity
kernel=deltav*pi*(fp(j,iap)+fp(l,iap))**2
elseif (lmaxwell_deltav) then ! maxwellian relative velocity
call particles_coag_maxwell(deltav, maxwell_param)
kernel=deltav*pi*(fp(j,iap)+fp(l,iap))**2
endif
! Calculate rates and cum. functions
r_ik_mat(j,l) = fp(l,inpswarm)*kernel ! rate matrix element
r_i_tot(j) = r_i_tot(j) + r_ik_mat(j,l) ! total rate for part. j
cum_func_sec_ik(j,l) = r_i_tot(j) ! cum. func. for phys. part.
!
l = l + 1 ! go to next particle
if(l == npar + 1) exit ! exit if last particle
enddo
r_total = r_total + r_i_tot(j) ! total rate for all particles
cum_func_first_i(j) = r_total ! cum. func. for rep. part.
!
l = 1
j = j + 1 ! go to next particle
if(j == npar + 1) exit
enddo
endif
!
! Now get the timestep (time until next collision) and put in
! inverse time-step array. Use eq. 5 i Zsom & Dullemond 2008
!
call random_number_wrapper(r)
dt1_coag_par = -r_total/log(r)
dt1_max = dt1_coag_par
!
endsubroutine particles_coag_timestep_zd
!***********************************************************************
subroutine particles_coag_outcome_zd(fp)
!
! Calculate which two particles collide and the outcome with the same
! scheme as in Zsom & Dullemond (2008). Only use for kernel tests.
! Assume perfect sticking for now.
!
! 15-nov-12/KWJ: coded
!
use General, only: random_number_wrapper
! Get all variables
!NILS: Shouldn't npar be mpar_loc?
! real, dimension (npar,mparray) :: fp
real, dimension (mpar_loc,mparray) :: fp
!
real :: r, r_i_old, kernel
integer :: i, j, k, l
real :: delta_r ! change in r_i rate
!
! Get the particles that is involved in the collision with
! random numbers and cumulative distribution functions.
! N.B. Using un-normalized cumulative functions.
! Use the bisection method to find the right particle.
!
! Find the representative particle, j.
call random_number_wrapper(r)
r = r*r_total
call particles_coagulation_bisection(cum_func_first_i, r, j)
! Find the physical particle, l.
call random_number_wrapper(r)
r = r*r_i_tot(j)
call particles_coagulation_bisection(cum_func_sec_ik(j,:), r, l)
!
! Change properties of the representative particle (particle j). Change
! mass (radius) and particle density in superparticle j.
!
fp(j,inpswarm)=fp(j,inpswarm)* &
(1.0/(1.0+(fp(l,iap)/fp(j,iap))**3))
fp(j,iap)=(fp(j,iap)**3+fp(l,iap)**3)**(1.0/3.0)
!
! Update r_ik matrix. Only need to update elements with representative
! particle i = j (column) and physical particle k = j (row).
!
! Loop over i and update row with k = j. Skip k = i = j
! (value in i = j column).
!
i = 1
do while (.true.)
if (i /= j) then
! Theoretical kernels with analytic solutions
if (lconstant_kernel_test) then
kernel=kernel_cst
elseif (llinear_kernel_test) then
kernel=kernel_lin* &
four_pi_rhopmat_over_three*(fp(i,iap)**3+fp(j,iap)**3)
elseif (lproduct_kernel_test) then
kernel=kernel_pro* &
four_pi_rhopmat_over_three2*fp(i,iap)**3*fp(j,iap)**3
! Physical kernels
elseif (lconstant_deltav) then ! constant relative velocity
kernel=deltav*pi*(fp(j,iap)+fp(l,iap))**2
elseif (lmaxwell_deltav) then ! maxwellian relative velocity
call particles_coag_maxwell(deltav, maxwell_param)
kernel=deltav*pi*(fp(j,iap)+fp(l,iap))**2
endif
!
delta_r = fp(j,inpswarm)*kernel - &
r_ik_mat(i,j) ! change in matrix element
r_ik_mat(i,j) = fp(j,inpswarm)*kernel ! rate matrix element
! Update rate and cumulative function
r_i_tot(i) = r_i_tot(i) + delta_r ! total rate for rep. particle i
! Only need to update cumulative function from element j and onwards
cum_func_sec_ik(i, j:) = cum_func_sec_ik(i, j:) + &
delta_r
endif
!
i = i + 1 ! go to next superparticle
if(i == npar + 1) exit
enddo
!
! Update i = j column.
!
k = 1
do while (.true.)
if (lconstant_kernel_test) then
kernel=kernel_cst
elseif (llinear_kernel_test) then
kernel=kernel_lin* &
four_pi_rhopmat_over_three*(fp(j,iap)**3+fp(k,iap)**3)
elseif (lproduct_kernel_test) then
kernel=kernel_pro* &
four_pi_rhopmat_over_three2*fp(j,iap)**3*fp(k,iap)**3
! Physical kernels
elseif (lconstant_deltav) then ! constant relative velocity
kernel=deltav*pi*(fp(j,iap)+fp(l,iap))**2
elseif (lmaxwell_deltav) then ! maxwellian relative velocity
call particles_coag_maxwell(deltav, maxwell_param)
kernel=deltav*pi*(fp(j,iap)+fp(l,iap))**2
endif
!
r_ik_mat(j,k) = fp(k,inpswarm)*kernel ! rate matrix element
! Update rate for representative particle j
if (k == 1) then
r_i_tot(j) = r_ik_mat(j,k)
else
r_i_tot(j) = r_i_tot(j) + r_ik_mat(j,k)
endif
cum_func_sec_ik(j,k) = r_i_tot(j)
!
k = k + 1
if (k == npar + 1) exit
enddo
!
! Update total rate and cumulative function for rep. particles
!
i = 1
do while (.true.)
if (i == 1) then
r_total = r_i_tot(i)
else
r_total = r_total + r_i_tot(i)
end if
cum_func_first_i(i) = r_total
!
i = i + 1
if (i == npar + 1) exit
end do
!
endsubroutine particles_coag_outcome_zd
!***********************************************************************
subroutine particles_coagulation_bisection(qArr, qVal, iPart)
!
! Given random value qVal [0,max_rate[ find particle iPart through the bisection
! method given cumulative function qArr.
!
! 15-nov-12/KWJ: coded
!
real, dimension (npar) :: qArr
real :: qVal
integer :: iPart, jl, ju, jm
!
intent (in) :: qArr,qVal
intent (out) :: iPart
!
jl = 1 ! lower index limit
ju = size(qArr) ! upper index limit
if (qVal <= qArr(1)) then ! qVal is in first bin
iPart = 1
else
do while ((ju-jl) > 1) ! not yet done
jm = (ju+jl)/2 ! midpoint
if (qVal >= qArr(jm)) then
jl = jm ! new lower limit
else
ju = jm ! or new upper limit
endif
enddo
iPart = ju ! qVal is in the ju-bin
endif
!
endsubroutine particles_coagulation_bisection
!***********************************************************************
subroutine particles_coag_maxwell_johnk(dv, alpha)
!
! Generate a Maxwell-Boltzmann distributed relative velocity dv with
! parameter alpha through the Johnk's algorithm.
!
! For particles following the gas: alpha = sqrt(k_B*T/m)
! Mean velocity = 2*sqrt(2/pi)*alpha
!
! 14-jan-13/KWJ: coded
!
use General, only: random_number_wrapper
!
real :: dv, alpha
real :: r, w, w1, w2
!
intent (in) :: alpha
intent (out) :: dv
!
call random_number_wrapper(r)
dv = -log(r)
!
do while (.true.)
call random_number_wrapper(r)
w1 = r*r
call random_number_wrapper(r)
w2 = r*r
w = w1 + w2
if(w <= 1.) exit
enddo
!
call random_number_wrapper(r)
dv = dv - w1/w*log(r)
!
dv = alpha*sqrt(2.*dv)
!
endsubroutine particles_coag_maxwell_johnk
!***********************************************************************
subroutine particles_coag_maxwell(dv, alpha)
!
! Generate a Maxwell-Boltzmann distributed relative velocity dv with
! parameter alpha through the algorithm described in Mohamed (2011).
!
! For particles following the gas: alpha = sqrt(k_B*T/m)
! Mean velocity = 2*sqrt(2/pi)*alpha
!
! 14-jan-13/KWJ: coded
!
use General, only: random_number_wrapper
!
real :: dv, alpha
real :: y, r, tmp1, tmp2
! For findig random nr.
! C = 0.5, y0 = 1, k = 2/(C*sqrt(pi))*y0^(1/2)*e^(-C*y0)
! k = 4/sqrt(e*pi) ~ 1.37
! g = 2/(k*C*sqrt(pi)) = sqrt(e) ~ 1.649
! real :: g = 1.648721271
real :: e = 2.718281828
! Only really interested in g^2 = e
!
intent (in) :: alpha
intent (out) :: dv
!
!
do while (.true.)
call random_number_wrapper(r)
y = -2.*log(r)
tmp1 = e*y*r*r
call random_number_wrapper(r)
tmp2 = r*r
if(tmp1.ge.tmp2) exit
enddo
!
dv = alpha*sqrt(2.*y)
!
endsubroutine particles_coag_maxwell
!***********************************************************************
subroutine read_particles_coag_run_pars(iostat)
!
use File_io, only: parallel_unit
!
integer, intent(out) :: iostat
!
read(parallel_unit, NML=particles_coag_run_pars, IOSTAT=iostat)
!
endsubroutine read_particles_coag_run_pars
!***********************************************************************
subroutine write_particles_coag_run_pars(unit)
!
integer, intent(in) :: unit
!
write(unit, NML=particles_coag_run_pars)
!
endsubroutine write_particles_coag_run_pars
!***********************************************************************
subroutine rprint_particles_coagulation(lreset,lwrite)
!
! Read and register diagnostic parameters.
!
! 28-mar-09/anders: adapted
!
use Diagnostics
!
logical :: lreset
logical, optional :: lwrite
!
integer :: iname
!
if (lreset) then
idiag_ncoagpm=0; idiag_ncoagpartpm=0; idiag_dt1_coag_par=0
idiag_k100_100=0; idiag_k100_80=0; idiag_k100_60=0; idiag_k100_50=0
idiag_k100_40=0; idiag_k100_30=0; idiag_k100_20=0; idiag_k100_10=0
idiag_k80_80=0; idiag_k80_60=0; idiag_k80_50=0; idiag_k80_40=0
idiag_k80_30=0; idiag_k80_20=0; idiag_k80_10=0
idiag_k60_60=0; idiag_k60_50=0; idiag_k60_40=0
idiag_k60_30=0; idiag_k60_20=0; idiag_k60_10=0
idiag_k50_50=0; idiag_k50_40=0; idiag_k50_30=0
idiag_k50_20=0; idiag_k50_10=0; idiag_k40_40=0
idiag_k40_30=0; idiag_k40_20=0; idiag_k40_10=0
idiag_k30_30=0; idiag_k30_20=0; idiag_k30_10=0
idiag_k20_20=0; idiag_k20_10=0; idiag_k10_10=0
endif
!
do iname=1,nname
call parse_name(iname,cname(iname),cform(iname),'ncoagpm',idiag_ncoagpm)
call parse_name(iname,cname(iname),cform(iname), &
'ncoagpartpm',idiag_ncoagpartpm)
call parse_name(iname,cname(iname),cform(iname),'dt1_coag_par',idiag_dt1_coag_par)
call parse_name(iname,cname(iname),cform(iname),'k100_100',idiag_k100_100)
call parse_name(iname,cname(iname),cform(iname),'k100_80',idiag_k100_80)
call parse_name(iname,cname(iname),cform(iname),'k100_60',idiag_k100_60)
call parse_name(iname,cname(iname),cform(iname),'k100_50',idiag_k100_50)
call parse_name(iname,cname(iname),cform(iname),'k100_40',idiag_k100_40)
call parse_name(iname,cname(iname),cform(iname),'k100_30',idiag_k100_30)
call parse_name(iname,cname(iname),cform(iname),'k100_20',idiag_k100_20)
call parse_name(iname,cname(iname),cform(iname),'k100_10',idiag_k100_10)
call parse_name(iname,cname(iname),cform(iname),'k80_80',idiag_k80_80)
call parse_name(iname,cname(iname),cform(iname),'k80_60',idiag_k80_60)
call parse_name(iname,cname(iname),cform(iname),'k80_50',idiag_k80_50)
call parse_name(iname,cname(iname),cform(iname),'k80_40',idiag_k80_40)
call parse_name(iname,cname(iname),cform(iname),'k80_30',idiag_k80_30)
call parse_name(iname,cname(iname),cform(iname),'k80_20',idiag_k80_20)
call parse_name(iname,cname(iname),cform(iname),'k80_10',idiag_k80_10)
call parse_name(iname,cname(iname),cform(iname),'k60_60',idiag_k60_60)
call parse_name(iname,cname(iname),cform(iname),'k60_50',idiag_k60_50)
call parse_name(iname,cname(iname),cform(iname),'k60_40',idiag_k60_40)
call parse_name(iname,cname(iname),cform(iname),'k60_30',idiag_k60_30)
call parse_name(iname,cname(iname),cform(iname),'k60_20',idiag_k60_20)
call parse_name(iname,cname(iname),cform(iname),'k60_10',idiag_k60_10)
call parse_name(iname,cname(iname),cform(iname),'k60_10',idiag_k60_10)
call parse_name(iname,cname(iname),cform(iname),'k50_50',idiag_k50_50)
call parse_name(iname,cname(iname),cform(iname),'k50_40',idiag_k50_40)
call parse_name(iname,cname(iname),cform(iname),'k50_30',idiag_k50_30)
call parse_name(iname,cname(iname),cform(iname),'k50_20',idiag_k50_20)
call parse_name(iname,cname(iname),cform(iname),'k50_10',idiag_k50_10)
call parse_name(iname,cname(iname),cform(iname),'k40_40',idiag_k40_40)
call parse_name(iname,cname(iname),cform(iname),'k40_30',idiag_k40_30)
call parse_name(iname,cname(iname),cform(iname),'k40_20',idiag_k40_20)
call parse_name(iname,cname(iname),cform(iname),'k40_10',idiag_k40_10)
call parse_name(iname,cname(iname),cform(iname),'k30_30',idiag_k30_30)
call parse_name(iname,cname(iname),cform(iname),'k30_20',idiag_k30_20)
call parse_name(iname,cname(iname),cform(iname),'k30_10',idiag_k30_10)
call parse_name(iname,cname(iname),cform(iname),'k20_20',idiag_k20_20)
call parse_name(iname,cname(iname),cform(iname),'k20_10',idiag_k20_10)
call parse_name(iname,cname(iname),cform(iname),'k10_10',idiag_k10_10)
enddo
!
if (present(lwrite)) call keep_compiler_quiet(lwrite)
!
endsubroutine rprint_particles_coagulation
!***********************************************************************
endmodule Particles_coagulation