subroutine nbi_distrsample(irsize,ierr)
  use nbi_mpi_com
  use nbi_random
  use nbi_random_pll
  use nbi_fbm_random !drndm is there
  use nbi_com, only: lsbeam,nznbmea,nznbmr,nthzsm,nznbma,bmvol,fbm_cx,lcentr,ncx0_use,xncx0_orb,wght_cx,fp
  use nbi_com, only: efbmb
  use nbi_com, only: snbx1d,nxskpb,lep1,nznbmri
  implicit none
!
  integer, intent(in) :: irsize  ! sample size
  integer, intent(out) :: ierr
  real(fp) :: dsum,dsum1,wtot,znorm_snb,zwrk,xzid,dsum_wght
  integer :: i,j,j1,ipos,ie,ix,ith1,ith2,ith,ia,indx
  integer :: idum,i_count
  real(fp), dimension(:,:), allocatable :: wrk
  integer, dimension(:), allocatable :: ix_zone
  logical :: il_wrk
  real(fp) :: de,da

!
!  not normalize beam distribution
!
  ierr=0
  
  znorm_snb=sum(snbx1d(:,lsbeam))
  if(znorm_snb.le.0) then
     ierr=ierr+1
     call errorlog( ' nbi_cxdistr: normalization <= 0.0')
     return
  endif


  il_wrk=.false.

  if(nxskpb.gt.1)then
     if(allocated(wrk)) deallocate(wrk)
     allocate(wrk(nznbmr,nxskpb))
     il_wrk=.true.
     wrk=0.0_fp
  endif

  if(allocated(ix_zone)) deallocate(ix_zone)
  allocate(ix_zone(irsize))
  ix_zone=0
 
!
!  get ntot
!
!
  
  da = 1.0_fp/nznbma
  dsum=0.0_fp
  dsum_wght=0.0_fp
  do ie=1,nznbmea(lsbeam)
     de=efbmb(ie+1,lsbeam)-efbmb(ie,lsbeam)
     do ix=1,nznbmr
        ith1=nthzsm(lcentr+ix-2)+1
        ith2=nthzsm(lcentr+ix-1)
        do ith=ith1,ith2
           do ia=1,nznbma
              dsum = dsum + fbm_cx(ie,ia,ith,lsbeam)*de*da*bmvol(ith)
              dsum_wght = dsum_wght + fbm_cx(ie,ia,ith,lsbeam)
           enddo
        enddo
     enddo
  enddo
!
! weight used for improving cx statistics
!
  wght_cx=dsum_wght/real(ncx0_use,fp)

!
! ascending sequence btw 0 and dsum
!
  drndm(:)=drndm(:)*dsum 
 
  if(il_wrk) then
     do i=1,nznbmr
        wrk(i,:)=0.0_fp
        j=min(lep1,lcentr+(i-1)*nxskpb) !analogy to ngc
        if(j.gt.lep1-nxskpb+1) exit !for exterior use default value for xzid=nxskpb/2._fp
        do ix=1,nxskpb
           j1=j+ix-1
           wrk(i,ix)= sum(snbx1d(j:j1,lsbeam))
        enddo
     enddo
  endif
!
!  compute distribution
!
  dsum1=0.0_fp
  ipos=1
  xz_id(:)=nxskpb/2._fp
  i_count=0
  do ie=1,nznbmea(lsbeam)
     de=efbmb(ie+1,lsbeam)-efbmb(ie,lsbeam)
     do i=1,nznbmr
        xzid=nxskpb/2._fp ! default value
        if(il_wrk) then
           i_count=i_count+1
           zwrk=nbrandm_pll(idum)
           if(wrk(i,nxskpb).gt.0.0_fp) then
              wrk(i,:)=wrk(i,:)/wrk(i,nxskpb)
              !now find zone id
              ix=1 
              do
                 if(zwrk.gt.wrk(i,ix))then
                    ix=ix+1
                    if(ix.gt.nxskpb)then
                       ierr=ierr+1
                       call errorlog( ' nbi_cxdistr: sampling error')
                       return
                    endif
                    cycle
                 else
                    exit
                 endif
              enddo
!              xzid=ix-0.5_fp !at zone center
              xzid=ix-1 
           endif
        endif
        
        ith1=nthzsm(lcentr+i-2)+1
        ith2=nthzsm(lcentr+i-1)
        do ith=ith1,ith2
           do ia=1,nznbma
              dsum1 = dsum1 + fbm_cx(ie,ia,ith,lsbeam)*de*da*bmvol(ith)
              do
                 if(ipos.gt.irsize) exit
                 if(drndm(ipos).gt.dsum1) exit
                 !
                 indx = indxs(ipos)
                 jb_a(indx)=ia
                 jb_e(indx)=ie
                 n_gc2(indx)=ith
                 xz_id(indx)=xzid
                 ix_zone(indx)=i
                 !
                 ipos=ipos+1
                 !
              enddo
              !
           enddo
        enddo
     enddo
  enddo
  !
!  xz_id(:)=nxskpb/2._fp
  if(allocated(wrk)) deallocate(wrk)
  if(ipos.ne.(irsize+1)) then
     ierr=ierr+1
     call errorlog( ' nbi_cxdistr: sampling error, #found .ne. #requested:')
  endif
  
end subroutine nbi_distrsample
