  subroutine nbi_zone_ind(zr,zphi,zz,n_j,n_iz)
    ! return transp 1d zone index and beam 2d zone index
    ! for a partilce posinition in tokamak cylindrical coord. system
    use iso_c_binding, only: fp => c_double
    use nbi_com, only : xminbm,lep1,lcentr,nzones,nznbmr,nthzsm,nlsym2b,&
         thbdy0,xpi
    use xstraln_calls, only : xs_fxget
    implicit none
    real(fp), intent(in) :: zr,zphi,zz !ptcl. coordinates in tokamak cylindrical coord. system
    integer, intent(out) :: n_j,n_iz! transp 1d zone and beam 2d zone
    real(fp) :: xi_flr,th_flr !rho and theta in (mag.) coordinates returned
    integer :: iwarn            ! out-of-range warning flag
    integer :: ibr,inz0,inz
    ! find (mag.) coordinates
    n_j=0
    n_iz=0
    call xs_fxget(1.0_fp,zphi,zr,zz,xi_flr,th_flr,iwarn)! zr and zz are given in [m]
    if(xi_flr.gt.xminbm) then ! boundary of xi-space is defined in bmzini
       n_iz=0
       n_j=lep1
    else
       n_j=min(lep1,int(lcentr+xi_flr*nzones))
       ibr=int(lcentr+nznbmr*(xi_flr/xminbm))
       inz0=nthzsm(ibr-1)
       inz=nthzsm(ibr)
       if(nlsym2b) then
!     midplane symmetry assumed; theta ranging from -pi to pi
          n_iz=min(inz,int(inz0+1+(abs(th_flr)/xpi)*(inz-inz0)))
       else
          !     no midplane symmetry
          n_iz=min(inz,int(inz0+1+((th_flr-thbdy0)*0.5_fp/xpi)& ! thbdy0 is defined in bmzini
               &             *(inz-inz0)))
       endif
    endif
    
  end subroutine nbi_zone_ind
