integer function izmap(id_mcgrid,rho,theta)

  implicit NONE

  ! MCgrid zone map, with assumptions:
  !   1. axisymmetry
  !   2. ignore zones outside plasma
  !   3. updown symmetric theta range: [0,pi]
  !   4. updown asymmetric theta range: [-pi,pi]

  integer, intent(in) :: id_mcgrid  ! grid id (axisymmetry assumed)
  real*8, intent(in) :: rho,theta   ! (rho,theta in)

  !------------------------------------
  !  info on zone rows

  logical, save :: init = .FALSE.
  integer, dimension(:), allocatable, save :: izns,iztot
  integer, save :: inznr,inznri,iudsym,inth0
  integer, save :: id_mcgrid_int = 0
  integer, save :: istep_int = 0

  !------------------------------------
  integer i,izr,inth,ith,idum
  real*8 ztheta
  real*8, parameter :: cpi = 3.1415926535897931D+00
  real*8, parameter :: c2pi = 6.2831853071795862D+00

  !  check if id matches the one stored; if not, re-initialize...

  if(id_mcgrid.ne.id_mcgrid_int) init = .FALSE.
  if(istep_int.ne.0) init = .FALSE.

  if(.not.init) then
     ! ** initialize **
     ! get no. of zones per row; add up cumulative number

     istep_int = 1

     call mcgrid_getnumr(id_mcgrid,inznr,inznri)
     if(inznr.eq.0) then
        izmap=0   ! invalid MCgrid, evidently.
        return
     endif
     call mcgrid_getsym(id_mcgrid,iudsym)
     if(allocated(izns)) deallocate(izns,iztot)
     allocate(izns(inznr),iztot(0:inznr))
     call mcgrid_getnumzns(id_mcgrid,izns)
     iztot(0)=0
     do i=1,inznr
        iztot(i)=iztot(i-1)+izns(i)
     enddo
     init = .TRUE.
     id_mcgrid_int=id_mcgrid
  endif

  izr = 1 + rho*inznri           ! zone row index -- inznri are inside bdy
  izr = max(1,min(inznri,izr))   ! forced inside plasma

  inth = izns(izr)               ! no. of poloidal zones, this row

  ztheta = theta                 ! force theta in range [-pi,pi
  do while ( ztheta < -cpi )
     ztheta = ztheta + c2pi
  enddo
  do while ( ztheta > cpi )
     ztheta = ztheta - c2pi
  enddo
  
  if(iudsym.eq.1) then
     ztheta = abs(ztheta)  ! updown symmetry: map to upper half-plane
     ith=1 + ztheta*inth/cpi
     ith = max(1,min(inth,ith))
  else
     ith=1 + (ztheta+cpi)*inth/c2pi
     ith = max(1,min(inth,ith))
  endif

  izmap = iztot(izr-1)+ith

end function izmap
