subroutine bcmprsn(ierb)
 
  use nbi_com
  use nbi_mpi_com 
  use logmod
  !
  ! dmc nov 2001 -- compression operator
  ! rewritten, taken outside particle orbit loop, for speed and to
  ! remove dependence on d/dt{b,g} data from driver code
  !
  ! method:
  !   (toroidal) flux surface(wb) is saved from previous step;
  !   (r,z) position is saved from previous step;
  !   mu and canonical momentum (pphi) are samed from previous step.
  !   so...
  !
  !     (a) get new xion -- normalized flux coordinate -- from tor. flux
  !     (b) find new th -- closest point on xion surface to (r,z) from
  !         previous step
  !     (c) based on mu and pphi conservation, update the velocity vector
  !
  !
  ! only "old" particles of specie lsbeam are considered
 
  !-----------------------------------------------------------
  !  local...
 
!============
! idecl:  explicitize implicit real declarations:
  implicit none
!============
  integer iptcl,j,jz,jtr,ind,ierr,ioff,ierb
  real(fp) zxi,zxip,zxion,zth,zoff,zdiff,zpphi
  real(fp) zrmj,zz,zdum(2,2),zgdum
  real(fp) zf,zvion,zksid,zvperp,zvpll
  real(fp) zepp,zep,zmo,zxhw
 
  real(fp) xjadep,zxista
  integer iuncon,idum,itrap
 
  integer icount(mtbl),imax,imax2,ivec,ipts
  integer, dimension(:,:), allocatable :: ilists
  real(fp), dimension(:,:), allocatable :: zrlist,zzlist,zphlist
  real(fp), dimension(:), allocatable :: zrvec,zzvec,zphvec
  real(fp), dimension(:), allocatable :: zdist,zthnear,zphnear,zthprev
  real(fp), dimension(:), allocatable :: ztha,zthb
  real(fp) :: tdepay_ind
  integer find_ptcl_ind      ! integer function to find ptcl.id#

  !  replace all __mcdebug with __DEBUG and debug compile to get extra plotting
  !  of pre- and post- compression densities.

#ifdef __mcdebug
  real(fp) :: zdenpre(mj,2),zdenpost(mj,2),zdeln(mj,2)
  character*80 :: zlabloc
#endif
  !-----------------------------------------------------------
#ifndef __MPI
      integer myid_debug, numproc_debug
!
! 
      
      if(debug_mpi.and.nltest_output) then ! this will be a path for nubeam_test, 
                                           ! when debugging mpi nubeam 
                                           ! for one ptcl.sublistby serial code 
                                           ! nubeam_test 

         myid_debug=myidd
         numproc_debug=numprocd
      endif
#endif

      call nbi_mpi_com_init

#ifndef __MPI
      if(debug_mpi.and.nltest_output) then
         myidd=myid_debug       ! do not redefine  myidd and nuprocd for 
         numprocd=numproc_debug ! debugging of mpi version by serial code 
                                ! nubeam_test
      endif
#endif

#ifdef __mcdebug

  zlabloc=' '
  write(zlabloc,1001) runid,lsbeam,abeams(lsbeam),xzbeams(lsbeam)
1001 format(1x,a,' species ',i1,'  a=',1pe10.3,' z=',1pe10.3)

  call denchk(zdenpre,'before compression')
#endif

  ierb=0    ! error check
  icount=0
  ipts=nznblo+1
 
  if(nuse(lsbeam).eq.0) return   ! no particles...
 
  ! ...get new xion for each ptcl and find max no. in a single "bin"
 
  do iptcl=1,nuse(lsbeam)
     ind=nbndex(iptcl,lsbeam)
     ind=find_ptcl_ind(numprocd,myidd,ind)
 
     ! ...only do "old" particles
     
     tdepay_ind=max(0.0_fp,tdepay(ind,lsbeam)-tbm1) !mg 11/2013

     if(tdepay_ind.ne.0.0_fp) cycle   ! tdepay.gt.0.0 for newly deposited
 
     zxi=sqrt(tfluxay(ind,lsbeam)/tflux)
     xiay(ind,lsbeam)=zxi
     j = min(nznblo, int(zxi/dxi(lcentr,2)) + 1 )
     icount(j)=icount(j)+1
  enddo
 
  ! ...find max list sizes
 
  imax=icount(1)
  imax2=icount(1)
  do j=2,nznblo
     imax=max(imax,icount(j))  ! max in any one zone
     imax2=max(imax2,(icount(j)+icount(j-1)))  ! max next to any one bdy
  enddo
 
  if(imax.eq.0) return   ! no "old" particles...
 
  allocate(ilists(imax,nznblo))
  allocate(zrlist(imax,nznblo))
  allocate(zzlist(imax,nznblo))
  allocate(zphlist(imax,nznblo))
 
  ilists=0
 
  ! ...sort particles into bin lists; get nearest (r,z,phi) for each one.
 
  icount=0
  do iptcl=1,nuse(lsbeam)
     ind=nbndex(iptcl,lsbeam)
     ind=find_ptcl_ind(numprocd,myidd,ind)
 
     ! ...only do "old" particles
     tdepay_ind=max(0.0_fp,tdepay(ind,lsbeam)-tbm1) !mg 11/2013
     if(tdepay_ind.ne.0.0_fp) cycle   ! tdepay.gt.0.0 for newly deposited
 
     zxi=xiay(ind,lsbeam)
     j = min(nznblo, int(zxi/dxi(lcentr,2)) + 1 )
     icount(j)=icount(j)+1
     ilists(icount(j),j)=ind
     zrlist(icount(j),j)=rmjionay(ind,lsbeam)
     zzlist(icount(j),j)=xzionay(ind,lsbeam)
     zphlist(icount(j),j)=0.0_fp
  enddo
 
  ! input vectors for search routine
 
  allocate(zrvec(imax2),zzvec(imax2),zphvec(imax2))
 
  ! output vectors for search routine -- clear to zero
 
  allocate(zdist(imax2),zthnear(imax2),zphnear(imax2))
  zdist=0.0_fp
  zthnear=0.0_fp
  zphnear=0.0_fp
 
  ! vectors for binned ptcl theta interpolation
 
  allocate(ztha(imax),zthb(imax))
! (following initialization by kumar-- needs to be checked)
  ztha=0.0_fp
  zthb=0.0_fp
 
  ! ...loop over "interior" boundaries
 
  do j=2,ipts
 
     jtr=lcentr+j-1
     zxi=xiblo(jtr)
     zxip=xiblo(jtr-1)
 
     jz=j-1
 
     zrvec=0.0_fp
     zzvec=0.0_fp
     zphvec=0.0_fp
 
     ivec=icount(jz)+icount(j)
 
     zrvec(1:icount(jz)) = zrlist(1:icount(jz),jz)*0.01_fp  ! m
     zzvec(1:icount(jz)) = zzlist(1:icount(jz),jz)*0.01_fp  ! m
     zphvec(1:icount(jz)) = zphlist(1:icount(jz),jz)
     if(j.le.nznblo) then
        zrvec(icount(jz)+1:ivec) = zrlist(1:icount(j),j)*0.01_fp  ! m
        zzvec(icount(jz)+1:ivec) = zzlist(1:icount(j),j)*0.01_fp  ! m
        zphvec(icount(jz)+1:ivec) = zphlist(1:icount(j),j)
     endif
 
     ! ...find nearest points on target surface
 
!    if((ivec.gt.0).and.(j.lt.nznblo)) then
     if((ivec.gt.0).and.(j.lt.ipts)) then
        call eq_bdfind(ivec,zrvec,zzvec,zphvec, zxi,  &
             zthnear, zphnear, zdist, ierr)
        if(ierr.ne.0) exit
        call ckrange(zthnear,ivec)
     endif
 
     ! ...nearest theta, outer surface, zone jz
     ! ...inner surface data from previous loop iteration (except 1st one)
 
     if(j.eq.2) then
        zthb(1:icount(jz))=zthnear(1:icount(jz))
        ztha(1:icount(jz))=zthb(1:icount(jz))
     else if(j.eq.ipts) then
        zthb(1:icount(jz))=ztha(1:icount(jz))
     else
        zthb(1:icount(jz))=zthnear(1:icount(jz))
     endif
 
     if(icount(jz).gt.0) then
 
        do iptcl=1,icount(jz)
 
           ind=ilists(iptcl,jz)
           zxion=xiay(ind,lsbeam)
           nbsce=nbscay(ind,lsbeam)
 
           jtr = min(lep1, (int(zxion/dxi(lcentr,2)) + lcentr))
           if(zxion.lt.1.0_fp) jtr=min(ledge,jtr)
 
           zf=(zxion-zxip)/(zxi-zxip)
 
           zdiff=(ztha(iptcl)-zthb(iptcl))
           if(zdiff.gt.xpi) then
              zoff=-twopi
              ioff=-1
           else if(zdiff.lt.-xpi) then
              zoff=twopi
              ioff=1
           else
              zoff=0
              ioff=0
           endif
           zth=(1.0_fp-zf)*(ztha(iptcl)+zoff)+zf*zthb(iptcl)
           if(ioff.ne.0) call ckrang0(zth)
           thay(ind,lsbeam)=zth

           !  if the particle is outside the plasma, we don't trust
           !  the fields:  so, leave the velocities unchanged.
 
           if(jtr.le.ledge) then
 
              if(sawflag) then
                 zpphi = pmechay(ind,lsbeam)  ! just conserve mech. momentum
              else
                 zpphi = pphiay(ind,lsbeam)   ! conserve canonical momentum
              endif

              call ptclrgen(zxion,zth,xmuay(ind,lsbeam),zpphi,  &
                   sawflag,zvion,zksid,zvperp,zvpll,zmo)
 
              zepp=vay(ind,lsbeam)**2
              vay(ind,lsbeam)=zvion
              xksidy(ind,lsbeam)=zksid
 
              zep=zvion**2
              pbcprss(jtr,lsbeam)=pbcprss(jtr,lsbeam)+ &
                   (zep-zepp)*v2toev*wghtay(ind,lsbeam)
 
              if(v2toev*zep .gt. edbfac*einjs(lsbeam)) then
                 if(orbrzv_rf_option .ne. 1.and.(ntest_max.eq.0.or.abs(ntest_max).eq.10)) then
                   ! error flag setting removed by kjm
                   ! ierb = ierb + 1
                    call errorlog(' bcmprsn: -- compression acceleration:')
                    call errorlog(' bcmprsn: iptcl ',iptcl)
                    call errorlog(' bcmprsn: zone ',jz)
                    call errorlog(' bcmprsn: index ',ind)
                    call errorlog(' bcmprsn: old energy: ',v2toev*zepp)
                    call errorlog(' bcmprsn: new energy: ', v2toev*zep)
 
                    write(nonlin,*) ' ?bcmprsn -- compression acceleration:'
                    write(nonlin,*) '  iptcl ',iptcl,' zone ',jz,' index ',ind
                    write(nonlin,*) '  old energy: ',v2toev*zepp, &
                         & ' new energy: ', v2toev*zep
                    if(ierb.eq.1) write(nonlin,*) '  max energy: ', &
                         &edbfac*einjs(lsbeam)
                 endif
              endif

              !  count mech. momentum difference as jxb torque
 
              zgdum=0.0_fp
              call nb0jxb(wghtay(ind,lsbeam),zep,zep,pmechay(ind,lsbeam),zmo, &
                   zgdum,zgdum,jtr,jtr,zxion,zxion,4,.false.)
 
              !  new <x> for averaged flow calculation
 
              if(nfast_map.eq.1) then
                 call nbmomry(zxion,zth,0,zrmj,zz,zdum) ! use bilinear maps
              else
                 call nbmomry_xs(zxion,zth,0,zrmj,zz,zdum) !use xplasma spline maps
              endif
              zxista = xjadep(zxion,zth,zrmj,zz,zvion,zvpll,  &
                   itrap,iuncon,zdum(1,1),zdum(2,1),zdum(1,2),idum,zxhw)
 
              if(zxista.ge.1.0_fp) zxista=max(1.01_fp,zxista) ! bind away from bdy
 
              xistay(ind,lsbeam)=zxista
 
           endif  ! ptcl inside plasma
 
        enddo  ! loop over ptcls
 
     endif  ! test if there are ptcls
 
! following change by kumar (replaced .lt. by .le.)
!    if(j.lt.nznblo) then
     if(j.le.nznblo) then
 
        ! save data for next iteration
 
        ztha(1:icount(j))=zthnear(icount(jz)+1:ivec)
 
     endif
 
  enddo  ! end of loop over bins
  !
 
  if(ierr.ne.0) then
     call errorlog(' bcmprsn: eq_bdfind error, exiting before completion.')
     write(nonlin,*) ' ?bcmprsn: eq_bdfind error, exiting before completion.'
  endif
 
  !  cleanup...............
 
  deallocate(ilists,zrlist,zzlist,zphlist)
  deallocate(zrvec,zzvec,zphvec)
  deallocate(zdist,zthnear,zphnear)
  deallocate(ztha,zthb)
 
#ifdef __mcdebug
  call denchk(zdenpost,'after compression')
  zdeln = zdenpost - zdenpre

  call r8_grafx3(xizc(lcentr),zdenpre(lcentr,1),zdenpost(lcentr,1), &
       zdeln(lcentr,1),nzones, &
       ' ','#/cm3',zlabloc, &
       'bcmprsn.f90 debug plot: n from xiay: pre, post, delta',' ')

  call r8_grafx3(xizc(lcentr),zdenpre(lcentr,2),zdenpost(lcentr,2), &
       zdeln(lcentr,2),nzones, &
       ' ','#/cm3',zlabloc, &
       'bcmprsn.f90 debug plot: n from xistay: pre, post, delta',' ')
  
#endif

contains
    subroutine ckrange(zthvec,ivec)
      !check/enforce range [-pi,pi] on vector of theta values
 
      integer, intent(in) :: ivec
      real(fp), intent(inout) :: zthvec(ivec)
      integer i
 
      do i=1,ivec
         call ckrang0(zthvec(i))
      enddo
    end subroutine ckrange
 
    subroutine ckrang0(zth)
      !check/enforce range [-pi,pi] on a single theta value
 
      real(fp), intent(inout) :: zth
 
      do while (zth < -xpi)
         zth = zth + twopi
      enddo
      do while (zth > xpi )
         zth = zth - twopi
      enddo
    end subroutine ckrang0
 
#ifdef __mcdebug
    subroutine denchk(zden,zlabel)
      
      real(fp) :: zden(mj,2)
      character*(*) :: zlabel

      zden = 0.0_fp

      do iptcl=1,nuse(lsbeam)
         ind=nbndex(iptcl,lsbeam)
         ind=find_ptcl_ind(numprocd,myidd,ind)
 
         ! ...only do "old" particles
         tdepay_ind=max(0.0_fp,tdepay(ind,lsbeam)-tbm1) !mg 11/2013
         if(tdepay_ind.ne.0.0_fp) cycle 
 
         zxi = xiay(ind,lsbeam)
         j = min(nznblo, int(zxi/dxi(lcentr,2)) + 1 )

         if((j.ge.lcentr).and.(j.le.ledge)) then
            zden(j,1) = zden(j,1) + wghtay(ind,lsbeam)
         endif

         zxi = xistay(ind,lsbeam)
         j = min(nznblo, int(zxi/dxi(lcentr,2)) + 1 )

         if((j.ge.lcentr).and.(j.le.ledge)) then
            zden(j,2) = zden(j,2) + wghtay(ind,lsbeam)
         endif

      enddo

      do j=lcentr,ledge
         zden(j,1)=zden(j,1)/dvolzc(j)
         zden(j,2)=zden(j,2)/dvolzc(j)
      enddo

      call r8_grafx2(xizc(lcentr),zden(lcentr,1),zden(lcentr,2),nzones, &
           ' ','#/cm3',zlabloc, &
           'bcmprsn.f90 debug plot: n from xiay, n from xistay',zlabel)
    end subroutine denchk
#endif
end subroutine bcmprsn
