!-----------------------------------------------------------------------
!  bbnpik -- discretize the fast ion distribution function
!
subroutine bbnpik(izone,is,inum,iord,zrand,zf,ivv,ich,igot)
  !
  use nbi_com
  !
  !-----------------------------------------------------------------------
  !  mod dmc 24 sep 1991 -- energy grid of fast ion dist. fcns now varies
  !  with fast ion species, and it may not be divided up in constant steps
  !  of d(energy) -- take this into account in the code
  !    the output array zf now includes a d(energy) term
  !-----------------------------------------------------------------------
  !
  !
  !  input
  !
  !============
  implicit none
  integer :: itest,itry,iran,ia,ie,ict,ipik
  !============
  real(fp) :: zmax,zsmall,zdf,znorm,zdfnorm,zfgot,zzf,zdelebd,zerr
  real(fp) :: znormf
  !============
  integer :: izone,is,inum
  !
  !    izone -- 2d beam zone index
  !    is    -- fast ion species index
  !    inum  -- number to discretize to (target)
  !
  !  common input
  !    fbm   -- distribution fcn being discretized
  !    bdens2 -- dist. fcn normalization
  !              sum fbm(*,*,izone,is)*delebd/nzmbma = bdens2(izone,is)
  !
  !  workspaces
  !
  integer :: iord(inum)    ! order index
  real(fp) :: zrand(inum)      ! ordered set of random no.s btw 0 and 1
  !
  !  output
  !
  real(fp) :: zf(inum)         ! fbm discretized, normalized, and packed
  integer :: ivv(inum)     ! energy zone index of packed fbm entry
  integer :: ich(inum)     ! pitch angle index of packed fbm entry
  !
  integer :: igot          ! no. of pts actually returned, .le.inum
  !
  !
  real(fp) :: zfint,zdfint,ztest
  real(fp) :: fbm_flr !function to choose distribution function
  ! depending on nlfbmflr, fbm if nlfbmflr=.false.
  ! fbm_ptcl if nlfbmflr=.true.
  !c
  !-----------------------------------------------------------------------
  !  executable code.
  !
  !  get sorted random number sequence
  !
  call franay(zrand,iord,inum,1)
  !
  !  bug fix dmc oct 1992 -- make sure all rn's are bound away from 1.0
  !
  zmax=1.0000000_fp
  zsmall=5.0e-7_fp
  itest=inum+1
  do 
     itest=itest-1
     zmax=zmax-zsmall
     if(zrand(iord(itest)).le.zmax) exit
     zrand(iord(itest))=zmax
  enddo
  !
  zdf=1.000_fp/real(nznbma)
  !
  znorm=bdens2(izone,is)/bmvol(izone)
  zdfnorm=zdf/znorm
  !
  !  re-entry point to correct slight normalization error
  !
  itry=0
  !
  do 
     itry=itry+1
     !
     igot=0
     !
     zfgot=0.0_fp
     !
     zfint=0.0_fp
     !
     iran=1
     ztest=zrand(iord(iran))
     !
     do ia=1,nznbma
        do ie=1,nznbmea(is)
           zzf=fbm_flr(ie,ia,izone,is,0)
           if(zzf.eq.0.0_fp) cycle
           zdelebd=efbmb(ie+1,is)-efbmb(ie,is)
           zdfint=zzf*zdelebd*zdfnorm
           zfint=zfint+zdfint
           !
           ict=0
           do
              if(zfint.lt.ztest) exit
              ict=ict+1
              if(ict.eq.1) then
                 igot=igot+1
                 zf(igot)=zdelebd
                 ivv(igot)=ie
                 ich(igot)=ia
              else
                 zf(igot)=zf(igot)+zdelebd
              end if
              zfgot=zfgot+zdf*zdelebd
              iran=iran+1
              if(iran.gt.inum) then ! normal exit
                 znormf=znorm/zfgot
                 do ipik=1,igot
                    zf(ipik)=znormf*zf(ipik)
                 end do
                 return
              endif
              ztest=zrand(iord(iran))
           enddo
        end do
     end do
     !
     !  possibly in trouble here -- cut some slack if real close ...
     !
     if(itry.ne.1) exit
     zerr=1.0_fp-zfint
     if(abs(zerr).ge.1.e-5_fp) exit
     zdfnorm=zdfnorm/zfint
  enddo
  !
  !  cannot exit loop this way -- it would mean that bdens2 does not
  !  give the summed fbm as required
  !
  igot=-1
  write(nonlin,*) ' bbnpik -- bdens2/fbm normalization error.'
  return
end subroutine bbnpik
! 11sep2004 fgtok
