
!
! --------------------- xstraln_2d_zone ----------------------
! return the beam zone containing this point for the monte carlo
! grid initialized in the xstraln_2d module -- returns 0 if the
! point is outside the rho=1 boundary
!
function xstraln_2d_zone_r8(zrho, ztheta_arg) result (ibzone)
  use xstraln_2d

  implicit none

  real*8, intent(in) :: zrho       ! flux zone coordinate
  real*8, intent(in) :: ztheta_arg ! poloidal angle coordinate (radians)

  integer :: ibzone       ! beam zone index returned by function
  integer :: i            ! flux zone index
  integer :: in           ! number of poloidal zones in row i
  integer :: ib           ! beam zone offset in this row

  real*8  :: zfac         ! interpolation factor
  real*8  :: ztheta       ! poloidal angle coordinate in -pi,pi range
  real*8  :: ztheta_min   ! lower bound on poloidal angle

  real*8, parameter :: zpi  = 3.1415926535897931D0
  real*8, parameter :: zpi2 = 2.D0*zpi
  real*8, parameter :: zpin = -zpi

  ibzone=0

  zfac = (zrho-rhobr(1))/(rhobr(inbr+1)-rhobr(1))
  i    = max(1,int(real(zfac*inbr+1)))       ! force >0

  if (i>inbr) return    ! past rho=1 boundary, return 0

  ! bring theta in -pi,pi range
  ztheta = ztheta_arg
  do while(ztheta>zpi) 
     ztheta=ztheta-zpi2  ! too large
  end do
  do while(ztheta<zpin) 
     ztheta=ztheta+zpi2  ! too small
  end do

  if (ilsym_grid) then
     ztheta_min = 0.D0
     ztheta     = abs(ztheta)  ! up down symmetric, beam zones over top half only
  else
     ztheta_min = zpin
  end if

  zfac = (ztheta-ztheta_min)/(zpi-ztheta_min)
  in   = inthzsm(i)-inthzsm(i-1)             ! inthzsm(0:iznbmr)
  ib   = max(1,min(in, int(real(in*zfac+1))))

  ibzone = inthzsm(i-1) + ib

end function xstraln_2d_zone_r8
