! $Id$
!
!  This modules takes care of instantaneous collisions between
!  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_collisions = .true.
!
! MVAR CONTRIBUTION 0
! MAUX CONTRIBUTION 0
!
!***************************************************************
module Particles_collisions
!
  use Cparam
  use Cdata
  use General, only: keep_compiler_quiet
  use Messages
  use Particles_cdata
  use Particles_map
  use Particles_mpicomm
  use Particles_sub
!
  implicit none
!
  include 'particles_collisions.h'
!
  real, pointer, dimension (:) :: tausp_species, tausp1_species
  real, pointer :: gravr
  real :: lambda_mfp_single=1.0, coeff_restitution=1.0
  real :: energy_gain_inelastic=0.0
  integer :: ncoll_max_par=-1, npart_max_par=-1, npart_max_par_2=0
  logical :: lcollision_random_angle=.false., lcollision_big_ball=.false.
  logical :: lshear_in_vp=.false., lkeplerian_flat=.false.
  logical :: ltauc_from_tauf=.false., lapproaching_collisions=.false.
  logical :: lstop_at_first_collision=.false.
  character (len=labellen) :: icoll='big-ball'
!
  integer :: idiag_ncollpm=0, idiag_npartpm=0, idiag_decollpm=0
!
  namelist /particles_coll_run_pars/ &
      lambda_mfp_single, coeff_restitution, icoll, lshear_in_vp, &
      ncoll_max_par, npart_max_par, lkeplerian_flat, ltauc_from_tauf, &
      lapproaching_collisions, lstop_at_first_collision
!
  contains
!***********************************************************************
    subroutine initialize_particles_collisions(f)
!
!  Perform any post-parameter-read initialization i.e. calculate derived
!  parameters.
!
!  07-oct-08/anders: coded
!
      use SharedVariables, only: get_shared_variable
!
      real, dimension (mx,my,mz,mfarray) :: f
!
      select case (icoll)
      case ('random-angle'); lcollision_random_angle=.true.
      case ('big-ball');     lcollision_big_ball=.true.
      case default
        if (lroot) print*, 'No such value for icoll: ', trim(icoll)
        call fatal_error('initialize_particles_collisions','')
      endselect
!
!  Allocate neighbour array necessary for identifying collisions.
!
      if (.not.allocated(kneighbour)) allocate(kneighbour(mpar_loc))
!
      if (lparticles_blocks) then
        if (npart_max_par/=-1 .and. (.not.lrandom_particle_blocks)) then
          if (lroot) then
            print*, 'initialize_particles_collisions: '// &
                'with npart_max_par/=-1 one should set'
            print*, '    lrandom_particle_blocks=T in particles_run_pars'
            print*, '    to avoid artificial sub domains in the collisions'
          endif
          call fatal_error('initialize_particles_collisions','')
        endif
      else
        if (npart_max_par/=-1 .and. (.not.lrandom_particle_pencils)) then
          if (lroot) then
            print*, 'initialize_particles_collisions: '// &
                'with npart_max_par/=-1 one should set'
            print*, '    lrandom_particle_pencils=T in particles_run_pars'
            print*, '    to avoid artificial sub domains in the collisions'
          endif
          call fatal_error('initialize_particles_collisions','')
        endif
      endif
!
!  The maximum number of collision partners must be an even number.
!
      if (npart_max_par/=-1) then
        if (mod(npart_max_par,2)/=0) &
            call fatal_error('initialize_particles_collisions', &
            'npart_par_max must be an even number')
        npart_max_par_2=npart_max_par/2
      endif
!
!  Get radial gravity from gravity module.
!
      if (lkeplerian_flat) call get_shared_variable('gravr',gravr)
!
!  Get friction time from dust particle module.
!
      if (ltauc_from_tauf) then
        call get_shared_variable( 'tausp_species', tausp_species)
        call get_shared_variable('tausp1_species',tausp1_species)
      endif
!
!  Friction time must be set when calculating collision time from friction
!  time.
!
      if (ltauc_from_tauf) then
        if (npar_species==1 .and. tausp1_species(1)==0.0) call fatal_error( &
            'initialize_particles_collisions', 'tausp must be set when '// &
            'calculating collision time from friction time')
      endif
!
      if (.not.(lcartesian_coords.and.(all(lequidist)))) call fatal_error( &
           'initialize_particles_collisions','collisions only implemented '// &
           'for Cartesian equidistant grids.')
!
      call keep_compiler_quiet(f)
!
    endsubroutine initialize_particles_collisions
!***********************************************************************
    subroutine particles_collisions_timestep(fp,ineargrid)
!
!  Time-step contribution from particle collisions.
!
!  30-nov-10/anders: dummy
!
      real, dimension (mpar_loc,mparray) :: fp
      integer, dimension (mpar_loc,3) :: ineargrid
!
      call keep_compiler_quiet(fp)
      call keep_compiler_quiet(ineargrid)
!
    endsubroutine particles_collisions_timestep
!***********************************************************************
    subroutine particles_collisions_pencils(fp,ineargrid)
!
!  Calculate collisions between superparticles 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 change the velocity vectors of the colliding particles
!  instantaneously.
!
!  23-mar-09/anders: coded
!
      use Diagnostics
      use EquationOfState, only: cs0, rho0
      use General
!
      real, dimension (mpar_loc,mparray) :: fp
      integer, dimension (mpar_loc,3) :: ineargrid
!
      real, dimension (3) :: xpj, xpk, vpj, vpk
      real :: deltavjk, tau_coll1, prob, r
      real :: espec, asemi, omega_orbit
      real :: tausp_j, tausp_k, tausp1_j, tausp1_k
      integer, dimension (nx) :: np_pencil
      integer :: l, j, k, npart_point, ncoll, ncoll_par, npart_par
      integer :: jspec, kspec
!
      intent (in) :: ineargrid
      intent (inout) :: fp
!
!  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.
!
!  Note: with npart_max_par>0, it is safest to also set
!  lrandom_particle_pencils=T. Otherwise there is a risk that a particle
!  always interacts with the same subset of other particles.
!
        call shepherd_neighbour_pencil(fp,ineargrid,kshepherd,kneighbour)
!
!  Calculate number of particles per grid point. This is only needed in order
!  to limit the number of collision partners. Note that f(:,:,:,inp) is not
!  up-to-date at this time.
!
        if (npart_max_par/=-1) then
          np_pencil=0
          if (npar_imn(imn)/=0) then
            do k=k1_imn(imn),k2_imn(imn)
              np_pencil(ineargrid(k,1)-nghost)= &
                  np_pencil(ineargrid(k,1)-nghost)+1
            enddo
          endif
        endif
!
        do l=l1,l2
          k=kshepherd(l-nghost)
          if (k>0) then
            if (npart_max_par/=-1) npart_point=np_pencil(l-nghost)-1
            do while (k/=0)
              j=k
              if (ltauc_from_tauf) then
                if (lparticles_radius) then
                  tausp_k=fp(k,iap)*rhopmat
                  tausp1_k=1/tausp_k
                else
                  kspec=npar_species*(ipar(k)-1)/npar+1
                  tausp_k=tausp_species(kspec)
                  tausp1_k=tausp1_species(kspec)
                endif
              endif
              npart_par=0
              ncoll_par=0
              do while (.true.)
                j=kneighbour(j)
                if (j==0) then
                  if (npart_max_par/=-1 .and. npart_max_par<npart_point) then
                    j=kshepherd(l-nghost)
                  else
                    exit
                  endif
                endif
                if (ltauc_from_tauf) then
                  if (lparticles_radius) then
                    tausp_j=fp(j,iap)*rhopmat
                    tausp1_j=1/tausp_j
                  else
                    jspec=npar_species*(ipar(j)-1)/npar+1
                    tausp_j=tausp_species(jspec)
                    tausp1_j=tausp1_species(jspec)
                  endif
                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 ((.not.lapproaching_collisions) .or. &
                    sum((vpk-vpj)*(xpk-xpj))<0.0) then
!
!  For Keplerian particle discs, where the scale height of the particles is
!  given by their velocity dispersion Hp~vrms/OmegaK, we can use the 2-D
!  approach suggested by Lithwick & Chiang (2007). The collision time scale is
!
!    tcol=1/(n*sigma*vrms)=lambda/vrms
!
!  The particle density is n~Sigma/Hp, giving
!
!    tcol=1/(Sigma*sigma*OmegaK)=lambda0/(Omega*dx)
!
!  Here we used that lambda0=1/(npswarm*sigma) has been calculated as if
!  the particle scale height was dx.
!
                  if (lkeplerian_flat) then
                    espec=sum(vpk**2)/2-gravr/sqrt(sum(xpk**2))
                    asemi=-gravr/(2*espec)
                    omega_orbit=sqrt(gravr/asemi**3)
                    deltavjk=omega_orbit*dx
                  else
                    deltavjk=sqrt(sum((vpk-vpj)**2))
                  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*sigma*dv) = lambda/dv
!
!  where lambda is the mean free path of a particle relative to a single
!  superparticle and sigma is the collisional cross section.
!
!  We can further write the collision time of a representative particle from
!  swarm k colliding with the particles of swarm j in terms of their friction
!  times.
!
!    tau_coll_k = 1/(nj*sigma_jk*dv_jk) = mj/(rhoj*sigma_jk*dv_jk)
!               = (4/3)*pi*rhopmat*aj^3/(rhoj*pi*(aj+ak)^2*dv_jk)
!               = (4/3)*rhopmat*aj^3/(rhoj*(aj+ak)^2*dv_jk)
!               = (4/3)*tau_fric_j*(rhog/rhoj)*(cs/dv_jk)*
!                     tau_fric_j^2/(tau_fric_j+tau_fric_k)^2
!
                  if (ltauc_from_tauf) then
                    if (lparticles_radius.and.lparticles_density) then
                      tau_coll1=0.75*min(tausp1_j,tausp1_k)*deltavjk/cs0* &
                          fp(j,irhopswarm)/rho0*(tausp_j+tausp_k)**2* &
                          min(tausp1_j,tausp1_k)**2
                    elseif (npar_species>1) then
                      tau_coll1=0.75*min(tausp1_j,tausp1_k)*deltavjk/cs0* &
                          rhop_swarm/rho0*(tausp_j+tausp_k)**2* &
                          min(tausp1_j,tausp1_k)**2
                    else
                      tau_coll1=3*tausp1_j*deltavjk/cs0*rhop_swarm/rho0
                    endif
                  else
!
!  If the mean free path is a user supplied constant, then we can readily
!  calculate the collision time-scale.
!
                    tau_coll1=deltavjk/lambda_mfp_single
                  endif
!
!  Exclude sink particles from collisions.
!
                  if (lsinkparticle_1 .and. (j==1 .or. k==1)) &
                      tau_coll1=0.0
                  if (lparticles_sink) then
                    if (fp(j,iaps)/=0.0 .or. fp(k,iaps)/=0.0) &
                        tau_coll1=0.0
                  endif
!
!  Increase collision rate artificially for fewer collisions.
!
                  if (npart_max_par/=-1 .and. npart_max_par<npart_point) then
                    tau_coll1=tau_coll1*npart_point/npart_max_par
                  endif
!
                  if (tau_coll1/=0.0) then
!
!  The probability for a collision in this time-step is dt/tau_coll.
!
                    prob=dt*tau_coll1
                    call random_number_wrapper(r)
                    if (r<=prob) then
                      call particle_collision(xpj,xpk,vpj,vpk,j,k)
                      if (lshear .and. lshear_in_vp) then
                        vpk(2)=vpk(2)+qshear*Omega*xpk(1)
                        vpj(2)=vpj(2)+qshear*Omega*xpj(1)
                      endif
                      fp(k,ivpx:ivpz)=vpk
                      fp(j,ivpx:ivpz)=vpj
                      ncoll=ncoll+1
                      ncoll_par=ncoll_par+1
                    endif
                  endif
                endif
                npart_par=npart_par+1
                if (ncoll_max_par/=-1 .and. ncoll_par==ncoll_max_par) exit
                if (npart_max_par/=-1 .and. npart_par==npart_max_par_2) exit
              enddo
              k=kneighbour(k)
!
!  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_ncollpm/=0) &
                    call sum_par_name((/float(ncoll_par)/),idiag_ncollpm)
                if (idiag_npartpm/=0) &
                    call sum_par_name((/float(npart_par)/),idiag_npartpm)
                if (idiag_decollpm/=0) &
                    call sum_name(energy_gain_inelastic/npar,idiag_decollpm)
              endif
!
            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_ncollpm/=0) &
              call sum_par_name(fp(1:npar_loc,ixp),idiag_ncollpm)
          if (idiag_npartpm/=0) &
              call sum_par_name(fp(1:npar_loc,ixp),idiag_npartpm)
          call save_name(0.0,idiag_decollpm)
        endif
      endif
!
    endsubroutine particles_collisions_pencils
!***********************************************************************
    subroutine particles_collisions_blocks(fp,ineargrid)
!
!  Calculate collisions between superparticles 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 change the velocity vectors of the colliding particles
!  instantaneously.
!
!  23-mar-09/anders: coded
!
      use Diagnostics
      use EquationOfState, only: cs0, rho0
      use General
!
      real, dimension (mpar_loc,mparray) :: fp
      integer, dimension (mpar_loc,3) :: ineargrid
!
      real, dimension (3) :: xpj, xpk, vpj, vpk
      real :: deltavjk, tau_coll1, prob, r
      real :: espec, asemi, omega_orbit
      real :: tausp_j, tausp_k, tausp1_j, tausp1_k
      integer, dimension (nxb,nyb,nzb) :: kshepherdb, np_block
      integer :: ix, iy, iz, iblock, ix0, iy0, iz0
      integer :: j, k, npart_point, ncoll, ncoll_par, npart_par
      integer :: jspec, kspec
!
      intent (in) :: ineargrid
      intent (inout) :: fp
!
!  Reset collision counter.
!
      ncoll=0
!
!  Block loop.
!
      do iblock=0,nblock_loc-1
!
!  Create list of shepherd and neighbour particles for each grid cell in the
!  current block.
!
!  Note: with npart_max_par>0, it is safest to also set
!  lrandom_particle_pencils=T. Otherwise there is a risk that a particle
!  always interacts with the same subset of other particles.
!
        call shepherd_neighbour_block(fp,ineargrid,kshepherdb,kneighbour,iblock)
!
!  Calculate number of particles per grid point. This is only needed in order
!  to limit the number of collision partners. Note that f(:,:,:,inp) is not
!  up-to-date at this time.
!
        if (npart_max_par/=-1) then
          np_block=0
          if (npar_iblock(iblock)/=0) then
            do k=k1_iblock(iblock),k2_iblock(iblock)
              ix0=ineargrid(k,1); iy0=ineargrid(k,2); iz0=ineargrid(k,3)
              np_block(ix0-nghostb,iy0-nghostb,iz0-nghostb)= &
                  np_block(ix0-nghostb,iy0-nghostb,iz0-nghostb)+1
            enddo
          endif
        endif
!
        do iz=n1b,n2b; do iy=m1b,m2b; do ix=l1b,l2b
          k=kshepherdb(ix-nghostb,iy-nghostb,iz-nghostb)
          if (k>0) then
            if (npart_max_par/=-1) &
                npart_point=np_block(ix-nghostb,iy-nghostb,iz-nghostb)-1
            do while (k/=0)
              j=k
              if (ltauc_from_tauf) then
                kspec=npar_species*(ipar(k)-1)/npar+1
                tausp_k=tausp_species(kspec)
                tausp1_k=tausp1_species(kspec)
              endif
              npart_par=0
              ncoll_par=0
              do while (.true.)
                j=kneighbour(j)
                if (j==0) then
                  if (npart_max_par/=-1 .and. npart_max_par<npart_point) then
                    j=kshepherdb(ix-nghostb,iy-nghostb,iz-nghostb)
                  else
                    exit
                  endif
                endif
                if (ltauc_from_tauf) then
                  jspec=npar_species*(ipar(j)-1)/npar+1
                  tausp_j=tausp_species(jspec)
                  tausp1_j=tausp1_species(jspec)
                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 ((.not.lapproaching_collisions) .or. &
                    sum((vpk-vpj)*(xpk-xpj))<0.0) then
!
!  For Keplerian particle discs, where the scale height of the particles is
!  given by their velocity dispersion Hp~vrms/OmegaK, we can use the 2-D
!  approach suggested by Lithwick & Chiang (2007). The collision time scale is
!
!    tcol=1/(n*sigma*vrms)=lambda/vrms
!
!  The particle density is n~Sigma/Hp, giving
!
!    tcol=1/(Sigma*sigma*OmegaK)=lambda0/(Omega*dx)
!
!  Here we used that lambda0=1/(npswarm*sigma) has been calculated as if
!  the particle scale height was dx.
!
                  if (lkeplerian_flat) then
                    espec=sum(vpk**2)/2-gravr/sqrt(sum(xpk**2))
                    asemi=-gravr/(2*espec)
                    omega_orbit=sqrt(gravr/asemi**3)
                    deltavjk=omega_orbit*dx
                  else
                    deltavjk=sqrt(sum((vpk-vpj)**2))
                  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*sigma*dv) = lambda/dv
!
!  where lambda is the mean free path of a particle relative to a single
!  superparticle and sigma is the collisional cross section.
!
!  We can further write the collision time of a representative particle from
!  swarm k colliding with the particles of swarm j in terms of their friction
!  times.
!
!    tau_coll_k = 1/(nj*sigma_jk*dv_jk) = mj/(rhoj*sigma_jk*dv_jk)
!               = (4/3)*pi*rhopmat*aj^3/(rhoj*pi*(aj+ak)^2*dv_jk)
!               = (4/3)*rhopmat*aj^3/(rhoj*(aj+ak)^2*dv_jk)
!               = (4/3)*tau_fric_j*(rhog/rhoj)*(cs/dv_jk)*
!                     tau_fric_j^2/(tau_fric_j+tau_fric_k)^2
!
                  if (ltauc_from_tauf) then
                    if (lparticles_radius.or.lparticles_number) then
                      if (lroot) print*, 'particles_collisions_blocks: ', &
                          'not implemented for variable particle radius '// &
                          'or particle number'
                      call fatal_error('particles_collisions_blocks','')
                    endif
                    if (npar_species>1) then
                      tau_coll1=0.75*min(tausp1_j,tausp1_k)*deltavjk/cs0* &
                          rhop_swarm/rho0*(tausp_j+tausp_k)**2* &
                          min(tausp1_j,tausp1_k)**2
                    else
                      tau_coll1=3*tausp1_j*deltavjk/cs0*rhop_swarm/rho0
                    endif
                  else
!
!  If the mean free path is a user supplied constant, then we can readily
!  calculate the collision time-scale.
!
                    tau_coll1=deltavjk/lambda_mfp_single
                  endif
!
!  Exclude sink particles from collisions.
!
                  if (lsinkparticle_1 .and. (j==1 .or. k==1)) &
                      tau_coll1=0.0
                  if (lparticles_sink) then
                    if (fp(j,iaps)/=0.0 .or. fp(k,iaps)/=0.0) &
                        tau_coll1=0.0
                  endif
!
!  Increase collision rate artificially for fewer collisions.
!
                  if (npart_max_par/=-1 .and. npart_max_par<npart_point) then
                    tau_coll1=tau_coll1*npart_point/npart_max_par
                  endif
!
                  if (tau_coll1/=0.0) then
!
!  The probability for a collision in this time-step is dt/tau_coll.
!
                    prob=dt*tau_coll1
                    call random_number_wrapper(r)
                    if (r<=prob) then
                      call particle_collision(xpj,xpk,vpj,vpk,j,k)
                      if (lshear .and. lshear_in_vp) then
                        vpk(2)=vpk(2)+qshear*Omega*xpk(1)
                        vpj(2)=vpj(2)+qshear*Omega*xpj(1)
                      endif
                      fp(k,ivpx:ivpz)=vpk
                      fp(j,ivpx:ivpz)=vpj
                      ncoll=ncoll+1
                      ncoll_par=ncoll_par+1
                    endif
                  endif
                endif
                npart_par=npart_par+1
                if (ncoll_max_par/=-1 .and. ncoll_par==ncoll_max_par) exit
                if (npart_max_par/=-1 .and. npart_par==npart_max_par_2) exit
              enddo
              k=kneighbour(k)
!
!  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_ncollpm/=0) &
                    call sum_par_name((/float(ncoll_par)/),idiag_ncollpm)
                if (idiag_npartpm/=0) &
                    call sum_par_name((/float(npart_par)/),idiag_npartpm)
                call save_name(energy_gain_inelastic/npar,idiag_decollpm)
              endif
!
            enddo
          endif
        enddo; enddo; 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_ncollpm/=0) &
              call sum_par_name(fp(1:npar_loc,ixp),idiag_ncollpm)
          if (idiag_npartpm/=0) &
              call sum_par_name(fp(1:npar_loc,ixp),idiag_npartpm)
          call save_name(0.0,idiag_decollpm)
        endif
      endif
!
    endsubroutine particles_collisions_blocks
!***********************************************************************
    subroutine particle_collision(xpj,xpk,vpj,vpk,j,k)
!
!  Calculate collision between two particles.
!
!  13-nov-09/anders: coded
!
      use General
      use Mpicomm
      use Sub
!
      real, dimension (3) :: xpj, xpk, vpj, vpk
      integer :: j,k
!
      real, dimension (3) :: vvcm, vvkcm, vvkcmnew, vvkcm_normal, vvkcm_parall
      real, dimension (3) :: tmp1, tmp2, nvec
      real :: theta_rot, phi_rot
!
      intent (inout) :: xpj, xpk, vpj,vpk
      intent (in) ::j, k
!
      if (ip<=6) then
        print*, 'particle_collision: collision between'// &
            ' superparticles ', ipar(k), ' and ', ipar(j)
        print*, 'particle_collision: xj, xk='
        print*, '  ', xpj, sqrt(sum(xpj**2))
        print*, '  ', xpk, sqrt(sum(xpk**2))
        print*, 'particle_collision: **before**'
        print*, 'particle_collision: vj, vk, vj+vk='
        print*, '  ', vpj, sqrt(sum(vpj**2))
        print*, '  ', vpk, sqrt(sum(vpk**2))
        print*, '  ', vpj+vpk
        print*, 'particle_collision: ej, ek, ej+ek='
        print*, '  ', 0.5*(sum(vpj**2)), 0.5*(sum(vpk**2)), &
            0.5*(sum(vpj**2))+0.5*(sum(vpk**2))
        print*, 'particle_collision: lj, lk, lj+lk='
        call cross(xpj,vpj,tmp1)
        print*, '  ', tmp1
        call cross(xpk,vpk,tmp2)
        print*, '  ', tmp2
        print*, '  ', tmp1+tmp2
      endif
!
!  Choose random unit vector direction in COM frame. Here theta_rot is the
!  pole angle and phi_rot is the azimuthal angle. While we can choose phi_rot
!  randomly, theta_rot must be chosen with care so that equal areas on the
!  unit sphere are equally probable.
!    (see http://mathworld.wolfram.com/SpherePointPicking.html)
!
      if (lcollision_random_angle) then
        call random_number_wrapper(theta_rot)
        theta_rot=acos(2*theta_rot-1)
        call random_number_wrapper(phi_rot)
        phi_rot  =phi_rot*2*pi
        vvcm=0.5*(vpj+vpk)
        vvkcm=vpk-vvcm
        vvkcmnew(1)=+sin(theta_rot)*cos(phi_rot)
        vvkcmnew(2)=+sin(theta_rot)*sin(phi_rot)
        vvkcmnew(3)=+cos(theta_rot)
!
!  Multiply unit vector by length of velocity vector.
!
        vvkcmnew=vvkcmnew* &
            sqrt(vvkcm(1)**2+vvkcm(2)**2+vvkcm(3)**2)
!
!  Inelastic collision leads to energy loss of both particles.
!
        if (coeff_restitution/=1.0) then
          if (energy_gain_inelastic/=impossible) &
              energy_gain_inelastic=energy_gain_inelastic- &
              (1.0-coeff_restitution**2)*sum(vvkcm_normal**2)
          vvkcmnew=vvkcmnew*coeff_restitution
        endif
!
!  Change velocity vectors in normal frame.
!
        vpk=+vvkcmnew+vvcm
        vpj=-vvkcmnew+vvcm
!
      endif
!
!  Alternative method for collisions where each particle is considered to be
!  a big ball stretching to the surface of the other particle. We can then
!  solve the collision problem with perfect conservation of momentum, energy
!  and angular momentum.
!
      if (lcollision_big_ball) then
        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
        if (ip<=6) then
          print*, 'particle_collision: nvec, vvkcm, vvkcm_normal, vvkcm_parall='
          print*, '  ', nvec
          print*, '  ', vvkcm
          print*, '  ', vvkcm_normal
          print*, '  ', vvkcm_parall
        endif
!
!  Inelastic collision leads to energy loss of both particles.
!
        if (coeff_restitution/=1.0) then
          if (energy_gain_inelastic/=impossible) &
              energy_gain_inelastic=energy_gain_inelastic- &
              (1.0-coeff_restitution**2)*sum(vvkcm_normal**2)
          vvkcm_normal=vvkcm_normal*coeff_restitution
        endif
        vpk=vvcm+vvkcm_parall-vvkcm_normal
        vpj=vvcm-vvkcm_parall+vvkcm_normal
      endif
      if (ip<=6) then
        print*, 'particle_collision: **after**'
        print*, 'particle_collision: vj, vk, vj+vk='
        print*, '  ', vpj, sqrt(sum(vpj**2))
        print*, '  ', vpk, sqrt(sum(vpk**2))
        print*, '  ', vpj+vpk
        print*, 'particle_collision: ej, ek, ej+ek='
        print*, '  ', 0.5*(sum(vpj**2)), 0.5*(sum(vpk**2)), &
            0.5*(sum(vpj**2))+0.5*(sum(vpk**2))
        print*, 'particle_collision: lj, lk, lj+lk='
        call cross(xpj,vpj,tmp1)
        print*, '  ', tmp1
        call cross(xpk,vpk,tmp2)
        print*, '  ', tmp2
        print*, '  ', tmp1+tmp2
      endif
!
!  Stop after one collision (for testing purposes).
!
      if (lstop_at_first_collision) call fatal_error('particle_collision', &
          'stopping after first collision')
!
    endsubroutine particle_collision
!***********************************************************************
    subroutine read_particles_coll_run_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=particles_coll_run_pars, IOSTAT=iostat)
!
    endsubroutine read_particles_coll_run_pars
!***********************************************************************
    subroutine write_particles_coll_run_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=particles_coll_run_pars)
!
    endsubroutine write_particles_coll_run_pars
!***********************************************************************
    subroutine rprint_particles_collisions(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_ncollpm=0; idiag_npartpm=0; idiag_decollpm=0
      endif
!
      do iname=1,nname
        call parse_name(iname,cname(iname),cform(iname),'ncollpm',idiag_ncollpm)
        call parse_name(iname,cname(iname),cform(iname),'npartpm',idiag_npartpm)
        call parse_name(iname,cname(iname),cform(iname), &
            'decollpm',idiag_decollpm)
      enddo
!
      if (present(lwrite)) call keep_compiler_quiet(lwrite)
!
    endsubroutine rprint_particles_collisions
!***********************************************************************
endmodule Particles_collisions