subroutine xstraln_2d_init(izp1,id_mcgrid,ierr)

  use xstraln_mod
  use xstraln_2d

  implicit NONE

  ! (dmc Feb 2002)
  ! initialize for 2d straight line tracker calculations

  integer, intent(in) :: izp1   ! number of rho surfaces
  integer, intent(in) :: id_mcgrid   ! xplasma id for 2d grid information
  integer, intent(out) :: ierr       ! completion code (0=OK)

  ! -------------local-----------------

  character*50 zmodel
  integer iaxisymm,iextend,iudsym,iznbmri,iznbmr,i,j,ichk,imax,inum
  integer, dimension(:), allocatable :: iregion
  real*8, dimension(:), allocatable :: zrhoi,zchii,zphii,zri,zzi
  real*8, parameter :: zpi = 3.1415926535897931D+00
  !-----------------------------------------------------------------

  call eq_get_lunerr(lunmsgs)

  ierr=0

  call eq_ident(zmodel,iaxisymm,iextend)

  if(iaxisymm.eq.0) then
     write(lunmsgs,*) ' ?xstraln_2d_init: xplasma model not axisymmetric.'
     write(lunmsgs,*) '  xplasma ident = "',zmodel,'".'
     ierr=1
  endif

  if(.not.no_fRZ) then
     if(iextend.eq.0) then
        write(lunmsgs,*) &
             ' ?xstraln_2d_init: xplasma model lacks scrape-off region.'
        write(lunmsgs,*) '  xplasma ident = "',zmodel,'".'
        ierr=1
     endif
  endif

  if(ffx_r8%init.ne.1) then
     write(lunmsgs,*) &
          ' ?xstraln_2d_init:  call xstraln_init and xstraln_setup first!'
     ierr=1
  endif

  if(ierr.ne.0) return

  !---------------------------------------------------
  ! axis location

  if(no_fRZ) then
     inum=10
     allocate(zrhoi(inum),zchii(inum),zphii(inum),zri(inum),zzi(inum))
     allocate(iregion(inum))
     do i = 1,inum
        zrhoi(i)=rhomin
        zphii(i)=0
        zchii(i)=-zpi+2*(i-1)*zpi/(inum) ! -pi to (inum-2)*pi/inum evenly...
     enddo
     call eq_rzget(inum,zrhoi,zchii,zphii,zri,zzi,iregion,ierr)
     if(ierr.ne.0) then
        write(lunmsgs,*) &
             ' ?xstraln_2d_init: unexpected xplama eq_rzget error.'
     endif

     Raxis = sum(zri)/inum
     Zaxis = sum(zzi)/inum

     deallocate(zrhoi,zchii,zphii,zri,zzi,iregion)
  else
     Raxis = sum(ffx_r8%R(1,1:(ntheta-1),1))/(ntheta-1)
     Zaxis = sum(ffx_r8%Z(1,1:(ntheta-1),1))/(ntheta-1)
  endif

  ! reference stepsizes for tracker

  dsref=Raxis/50
  dsmin=Raxis/1000

  !---------------------------------------------------
  ! radial grid (rho) -- evenly spaced rhomin to rhobrk

  nzones=izp1-1
  nrho_target=izp1
  if(allocated(rho_target)) deallocate(rho_target)
  allocate(rho_target(nrho_target))
  do i=1,izp1
     rho_target(i)=((i-1)*rhobrk+(izp1-i)*rhomin)/nzones
  enddo

  ! beam zone rows -- coarsely spaced radial grid

  call mcgrid_getnumr(id_mcgrid,iznbmr,iznbmri)
  ichk=nzones/iznbmri
  if(ichk*iznbmri.ne.nzones) then
     write(lunmsgs,*) &
          ' ?xstraln_2d_init: radial grid / beam zone grid mismatch:'
     write(lunmsgs,*) &
          '  nzones = ',nzones,'  iznbmri = ',iznbmri
     ierr=2
     return
  endif

  inbr=iznbmri

  if(allocated(rhobr)) deallocate(rhobr)
  allocate(rhobr(inbr+1))
  rhobr(1)=rhomin
  do i=1,iznbmri
     rhobr(1+i)=rho_target(1+ichk*i)
  enddo

  ! check updown symmetry

  call mcgrid_getsym(id_mcgrid,iudsym)
  ilsym_grid = (iudsym.eq.1)

  ! get the numbers of zones (sums)

  if(allocated(inthbr)) deallocate(inthbr)
  allocate(inthbr(inbr))

  if(allocated(inthzsm)) deallocate(inthzsm)
  allocate(inthzsm(0:iznbmr))
  inthzsm(0)=0
  call mcgrid_getnumzns(id_mcgrid,inthzsm(1:iznbmr))

  imax=0
  do i=1,inbr
     inthzsm(i)=inthzsm(i-1)+inthzsm(i)
     if(ilsym_grid) then
        inthbr(i)=2*(inthzsm(i)-inthzsm(i-1))+1
     else
        inthbr(i)=(inthzsm(i)-inthzsm(i-1))+1
     endif
     imax=max(imax,inthbr(i))
  enddo

  ! set up arrays of theta boundary locations and zone indices within

  if(allocated(thbr)) deallocate(thbr)
  if(allocated(ibzns)) deallocate(ibzns)

  allocate(thbr(imax,inbr)); thbr=0
  allocate(ibzns(imax-1,inbr)); ibzns=0

  do i=1,inbr
     do j=1,inthbr(i)-1
        thbr(j,i)= (j-1)*twopi/(inthbr(i)-1) - xpi
     enddo
     thbr(inthbr(i),i)=xpi
     inum=inthzsm(i)-inthzsm(i-1)
     do j=1,inum
        if(ilsym_grid) then
           ibzns(inum+j,i)=inthzsm(i-1)+j
           ibzns(inum+1-j,i)=ibzns(inum+j,i)
        else
           ibzns(j,i)=inthzsm(i-1)+j
        endif
     enddo
  enddo

  if(allocated(izonrb)) deallocate(izonrb)
  allocate(izonrb(0:nrho_target))
  izonrb(0)=0
  izonrb(nrho_target)=0
  do i=1,inbr
     izonrb(ichk*(i-1)+1:ichk*i)=i   ! zone row -> beam zone row mapping
  enddo

  ! at the outer boundary of the first zone, get the (R,Z) coordinates
  ! at the theta boundaries

  inum=inthbr(1)

  allocate(zrhoi(inum),zchii(inum),zphii(inum),iregion(inum))
  zrhoi = rhobr(2)
  zchii = thbr(1:inum,1)
  zphii = 0

  if(allocated(rdel1)) deallocate(rdel1)
  if(allocated(zdel1)) deallocate(zdel1)
  allocate(rdel1(inum),zdel1(inum))

  call eq_rzget(inum,zrhoi,zchii,zphii,rdel1,zdel1,iregion,ierr)

  deallocate(zrhoi,zchii,zphii,iregion)

  rdel1=rdel1-Raxis
  zdel1=zdel1-Zaxis

end subroutine xstraln_2d_init
