!******************** start file nbkill.for ; group fimain ******************
!--------------------------------------------------------------
!  nbkill
!
!  russian roulette old ptcls to make room for new freshly deposited
!  ptcls -- probability of death = inkill/nuse(lsbeam) < 1.0
!
subroutine nbkill_pll(iorbrr,inkill)
  !  inkill-- number of ptcls to kill (approx)
  !
  use nbi_com
  use nbi_mpi_com
  use nbi_random
  use nbi_random_nbkill
  use logmod
  !
  implicit none
  !
  !  russian roulette bookkeeping arrays
  !
  integer, intent(in) :: inkill ! approximate number to kill (dep rr)
  logical, intent(in) :: iorbrr ! option to apply rr to underweight ions
  !
  integer :: imc,idum,ilp,jloc,indnew,iscreen1,iseed,&
       idt,inuseyd,mpierr

  real(fp) :: zkill,znkil,zwght,zwnew,zsurvive,zsurvi
  real(fp) :: zrandom
  integer :: find_ptcl_ind

  !--------------------------------------------------------------
  !  1.  option to apply russian roulette to underweight ions.
  !      each ion is subject to a 50% chance of removal; survivors
  !      have their weight doubled.
  !

  iscreen1=1600+myidd
  zrandom=nbrandm(idum)
  iseed=abs(15432000*zrandom)*lsbeam
  
  !
  ! diagnostics (kumar 09-22-08)
  !
  !     call infolog('nbkill_pll: 1 - random seed= ',iseed)

  zkill = 0.5_fp
  if(iorbrr) then

     ilp=0
     do
        ilp=ilp+1
        if(ilp.gt.nuse(lsbeam)) exit
        imc=nbndex(ilp,lsbeam)
        imc=find_ptcl_ind(numprocd,myidd,imc)
        !
        zwght=wghtay(imc,lsbeam)
        indnew=ilp+myidd*1.0e+6_fp
        call nbi_random_init_nbkill(iseed,indnew,iscreen1)
        jloc=njavay(imc,lsbeam)
        if(zwght.le.1.0e-4_fp*wbav_kill(jloc,lsbeam)) then
           ! mg 07/10/2013            
           ! if(zwght.le.0.3333333333333_fp*wbav_kill(jloc,lsbeam)) then
           ! mg must be consistent with orball.for
           !
           !  underweight ions only...
           !
           if(nbrandm_nbkill(idum).le.zkill) then
              call idelet(nbndex(1,lsbeam),minb,&
                   nuse(lsbeam),ilp)
              wnkill(lsbeam)=wnkill(lsbeam)+1.0_fp
              sforbrr(lsbeam)=sforbrr(lsbeam)+zwght ! weight killed
           else
              wghtay(imc,lsbeam)=2*zwght
              sforbrr(lsbeam)=sforbrr(lsbeam)-zwght ! weight added
           end if
           !
        end if
     end do
     !
  elseif((.not.iorbrr).and.(inkill.le.0)) then
     return
  end if

  !--------------------------------------------------------------
  !  2.  fraction of particles to kill -- make room for newly deposited
  !      ions
  !

  !   first, split inkill particles among processors
  !   participating in mpi

#ifdef __MPI

  stime41=mpi_wtime()

  mpierr=0
  inuseyd=0
  call mpi_allreduce(nuse(lsbeam),inuseyd,1,MPI_INTEGER,&
       MPI_SUM,nbi_comm,mpierr)

  puseyp(lsbeam)=inuseyd

  stime42=mpi_wtime()
  stime4=stime4+(stime42-stime41)

#else
  puseyp(lsbeam)=nuse(lsbeam)

#endif

  zrandom=nbrandm(idum)
  iseed=abs(18111100*zrandom)*lsbeam

  !
  ! diagnostics (kumar 09-22-08)
  !
  !     call infolog('nbkill_pll: 2 - random seed= ',iseed)


  if(inkill.le.0) return

  !
  zkill=inkill
  zkill=zkill/puseyp(lsbeam)
  zsurvive=1.0_fp-zkill
  zsurvi=1.0_fp/zsurvive
  !
  !  total weight killed
  !  define number of particles to be killed in each zone
  !  choose particle from population at zone j
  !
  ilp=0
  do
     ilp=ilp+1
     if(ilp.gt.nuse(lsbeam)) exit
     imc=nbndex(ilp,lsbeam)
     imc=find_ptcl_ind(numprocd,myidd,imc)
     !
     zwght=wghtay(imc,lsbeam)
     indnew=ilp+myidd*1.0e+6_fp
     call nbi_random_init_nbkill(iseed,indnew,iscreen1)
     !  search for particle
     if(nbrandm_nbkill(idum).le.zkill) then
        !           imc=nbndex(ilp,lsbeam)
        !           zwght=wghtay(imc,lsbeam)
        call idelet(nbndex(1,lsbeam),minb,&
             nuse(lsbeam),ilp) ! ilp modified
        sfdeprr(lsbeam)=sfdeprr(lsbeam)-zwght ! weight removed
        wkill(lsbeam)=wkill(lsbeam)+zwght
        wnkill(lsbeam)=wnkill(lsbeam)+1.0_fp
     end if
     !  keep looking
  end do

#ifdef __MPI

  stime51=mpi_wtime()

  mpierr=0
  inuseyd=0
  call mpi_allreduce(nuse(lsbeam),inuseyd,1,MPI_INTEGER,&
       MPI_SUM,nbi_comm,mpierr)

  puseyp(lsbeam)=inuseyd

  stime52=mpi_wtime()
  stime5=stime5+(stime52-stime51)

#else
  puseyp(lsbeam)=nuse(lsbeam)
#endif

#ifdef __DEBUG
  call infolog('nbkill_pll: nbkill(ii loop) lsbeam = ',lsbeam)
  call infolog('nbkill_pll: nbkill(ii loop) puseyp(lsbeam)', &
       int(puseyp(lsbeam)))
#endif

  return
end subroutine nbkill_pll
!******************** end file nbkill.for ; group fimain ******************
