subroutine xstraln_init(inphi,inrho,intheta,inr,inz,ierr)

  use xstraln_mod
  use xplasma_calls

  implicit NONE

  !
  !  initialize the xstraln (private) module
  !
  !----------------------------------------------------------
  !  passed:  grid sizes
  !    modification (dmc Jan 2006 -- if inr=inz=0 are passed, then
  !    xstraln is put into "no f(R,Z) mode".  This means, no bilinear
  !    inverse maps {rho(R,Z),theta(R,Z)}; it will make calculation
  !    of individual tracks more expensive, but it means that xplasma
  !    doesn't have to set up the inverse map which is itself somewhat
  !    expensive in compute time.
  !
  !    Bottom line on inr=inz=0:  if for a given geometry only a few
  !    tracks are going to be computed, 
  !

  integer, intent(in) :: inphi,inrho,intheta,inR,inZ

  !  completion code is returned >0 if a passed grid size is too small
  !  ierr=-1 means xplasma poloidal angle grid not found
  !  ierr=-2 means rhomax exceeds xplasma mapped region
  !  ierr=-3 means in xplasma something needed is not defined.
  !  ierr=-4 means some or all of {BR,BZ,Bphi} are undefined
  !          (this is allowed-- fast lookups of BR,BZ,Bphi are suppressed).

  integer, intent(out) :: ierr   ! completion code

  !----------------------------------------------------------
  !  local:

  real*8, dimension(:), allocatable :: zdata,zrho,zth,zphi,zr,zz
  integer, dimension(:), allocatable :: nregion
  real*8 bvec(3,1),zrmin,zrmax,zzmin,zzmax,zaxis

  real*8 ztol,zrho_rmin,zrho_rbrk1,zrho_rbrk2,zrho_rmax

  integer :: id_theta=0,id_phi=0,itype,isize,igot,i,id=0
  integer, parameter :: iscan=200, iscan2=20

  !----------------------------------------------------------
  !  executable code.
  !
  !  free previously allocated memory...
  !

  call free

  !
  !  grid sizes
  !

  ierr=0

  nphi=max(1,inphi)            ! toroidal angle
  if(inphi.lt.0) ierr=ierr+1

  nrho=max(20,inrho)            ! radial flux coordinate
  if(inrho.lt.nrho) ierr=ierr+1

  ntheta=max(20,intheta)        ! poloidal angle
  if(intheta.lt.ntheta) ierr=ierr+1

  if(inr.eq.0) then

     no_fRZ = .TRUE.
     nr = 0
     nz = 0

  else
     no_fRZ = .FALSE.

     nr=max(20,inr)                ! R
     if(inr.lt.nr) ierr=ierr+1

     nz=max(20,inz)                ! Z
     if(inz.lt.nz) ierr=ierr+1
  endif

  !
  !  grid arrays
  !

  if(allocated(phigrid)) deallocate(phigrid)
  if(nphi.gt.1) allocate(phigrid(nphi))

  if(allocated(thgrid)) deallocate(thgrid)
  allocate(thgrid(ntheta))

  if(allocated(rhogrid)) deallocate(rhogrid)
  allocate(rhogrid(nrho))

  if(allocated(rgrid)) deallocate(rgrid)
  if(.not.no_fRZ) allocate(rgrid(nr))

  if(allocated(zgrid)) deallocate(zgrid)
  if(.not.no_fRZ) allocate(zgrid(nz))

  !
  !  angle grid ids
  !

  call eq_ganum('__CHI',id_theta)

  call eq_ganum('__PHI',id_phi)

  if(id_theta.lt.0) ierr=-1

  if(ierr.ne.0) return

  !
  !  grid limits & grids
  !
  !
  !  ...rho grid

  call eq_rholim(rhomin,rhobrk)
  if(no_fRZ) then
     rhomax=rhobrk
     nrho_brk=nrho
  else
     rhomax=1.35d0*rhobrk             ! ** temporary hack ** (want bmzini grid)
     nrho_brk=(rhobrk-rhomin)*(nrho+1)/(rhomax-rhomin)
  endif

  delrho_hinv1=(nrho_brk-1)/(rhobrk-rhomin)
  if(no_fRZ) then
     delrho_hinv2 = 0
  else
     delrho_hinv2=(nrho-nrho_brk)/(rhomax-rhobrk)
  endif
  
  do i=1,nrho_brk
     rhogrid(i)=(rhomin*(nrho_brk-i)+rhobrk*(i-1))/(nrho_brk-1)
  enddo
  if(.not.no_fRZ) then
     do i=nrho_brk+1,nrho
        rhogrid(i)=(rhobrk*(nrho-i)+rhomax*(i-nrho_brk))/(nrho-nrho_brk)
     enddo
  endif

  !
  !  ...theta grid

  call eq_ngrid(id_theta,isize)

  allocate(zdata(isize))

  call eq_grid(id_theta,zdata,isize,igot,ierr)
  if(ierr.ne.0) call errmsg_exit(' ?xstraln_init fatal code error no.1')

  thmin=-xpi
  thmax= xpi
  th_period=zdata(isize)-zdata(1)
  if(abs(th_period-twopi).gt.0.001d0) then
     call errmsg_exit(' ?xstraln_init: xplasma chi grid range error.')
  endif
  delth_hinv=(ntheta-1)/th_period

  deallocate(zdata)
  
  do i=1,ntheta
     thgrid(i)=(thmin*(ntheta-i)+thmax*(i-1))/(ntheta-1)
  enddo

  !
  !  ...phi grid

  if(id_phi.gt.0) then

     call eq_ngrid(id_phi,isize)

     allocate(zdata(isize))

     call eq_grid(id_phi,zdata,isize,igot,ierr)
     if(ierr.ne.0) call errmsg_exit(' ?xstraln_init fatal code error no.2')

     phimin=-xpi
     phimax= xpi
     phi_period=zdata(isize)-zdata(1)
     if(abs(phi_period-twopi).gt.0.001d0) then
        call errmsg_exit(' ?xstraln_init: xplasma phi grid range error.')
     endif
     delphi_hinv=(nphi-1)/phi_period

     deallocate(zdata)

  else

     phimin=thmin
     phimax=thmax
     phi_period=th_period
     delphi_hinv=delth_hinv

  endif
  
  if(nphi.gt.1) then
     do i=1,nphi
        phigrid(i)=(phimin*(nphi-i)+phimax*(i-1))/(nphi-1)
     enddo
  endif

  !
  !  ...R,Z limits:  base on min/max R,Z
  !

  if(no_fRZ) then

     ! rmin/rmax/zmin/zmax based on limiter specifications
     ! even though the xplasma grid may not be extrapolated, the limiter
     ! locations must be supplied as this determines scrape-off of straight
     ! line tracks...

     call eq_bdlims(itype,zrmin,zrmax,zzmin,zzmax,ierr)

     ! impose a small expansion
     rmin = zrmax - 1.005d0*(zrmax-zrmin)
     rmax = zrmin + 1.005d0*(zrmax-zrmin)
     zmin = zzmax - 1.005d0*(zzmax-zzmin)
     zmax = zzmin + 1.005d0*(zzmax-zzmin)

     if(ierr.ne.0) then
        ierr=-5
        return
     endif

     no_BRZ = .TRUE.

     delz_hinv = 0

     delr_hinv1 = 0
     delr_hinv2 = 0
     delr_hinv3 = 0

     nbrk1 = 0
     nbrk2 = 0

     dlgrid_min = 0
     dlgrid_max = 0

  else
     call eq_rzlims(rmin,rmax,zmin,zmax,ierr)   ! xplasma (R,Z) grid
     if(ierr.ne.0) then
        ierr=-2
        return
     endif

     call eq_glimrz(rhomax,zrmin,zrmax,zzmin,zzmax,ierr) ! rhomax (R,Z) extent
     if(ierr.ne.0) then
        ierr=-2
        return
     endif

     !
     !  ...Z grid

     zmin=max(zmin,zzmin)
     zmax=min(zmax,zzmax)
     delz_hinv=(nz-1)/(zmax-zmin)

     do i=1,nz
        zgrid(i)=(zmin*(nz-i)+zmax*(i-1))/(nz-1)
     enddo

     !
     !  ...R grid

     rmin=max(rmin,zrmin)      ! choose intersection of above regions
     rmax=min(rmax,zrmax)

     !
     !  ...decide on break points based on axis and outer bdy position
     !
     !  the R grid is partitioned into three equally spaced subregions
     !  spanning:  [Rmin,Raxis], [Raxis,Rbdy_outer], [Rbdy_outer,Rmax]
     !    where Rmin is the minimum R in the map;
     !          Rmax is the maximum R in the map;
     !          Raxis is the mag. axis location
     !          Rbdy_outer is the plasma bdy midplane intercept, large R side
     !

     allocate(zrho(iscan),zth(iscan),zphi(iscan),zr(iscan),zz(iscan))
     allocate(nregion(iscan))

     do i=1,iscan
        zrho(i)=rhomin
        zth(i)=(thmin+thmax)/2
        zphi(i)=(phimin*(iscan-i)+phimax*(i-1))/(iscan-1)
     enddo

     call eq_rzget(iscan,zrho,zth,zphi,zr,zz,nregion,ierr)
     if(ierr.ne.0) then
        ierr=-3
        rbrk1=(rmin+rmax)/2
        zaxis=(zmin+zmax)/2
     else
        rbrk1=sum(zr)/iscan
        zaxis=sum(zz)/iscan
     endif

     call eq_glimrz(rhobrk,zrmin,zrmax,zzmin,zzmax,ierr)
     if(ierr.ne.0) then
        ierr=-3
        rbrk2=(rbrk1+rmax)/2
     else
        rbrk2=zrmax
     endif

     !  ...get rho values at R brk pts; decide how to partition R grid

     zz(1:iscan)=zaxis
     ztol=1.0d-2  ! polar map OK

     zr(1:iscan)=rmin
     call eq_inv(iscan,zr,zz,zphi,ztol,zrho,zth,nregion,ierr)
     if(ierr.eq.0) then
        zrho_rmin=-sum(zrho)/iscan
     else
        zrho_rmin=-rhomax
     endif

     zrho_rbrk1=rhomin

     zr(1:iscan)=rbrk2
     call eq_inv(iscan,zr,zz,zphi,ztol,zrho,zth,nregion,ierr)
     if(ierr.eq.0) then
        zrho_rbrk2=sum(zrho)/iscan
     else
        zrho_rbrk2=rhobrk
     endif

     zr(1:iscan)=rmax
     call eq_inv(iscan,zr,zz,zphi,ztol,zrho,zth,nregion,ierr)
     if(ierr.eq.0) then
        zrho_rmax=sum(zrho)/iscan
     else
        zrho_rmax=rhomax
     endif

     nbrk1=(zrho_rbrk1-zrho_rmin)*(nr+1)/(zrho_rmax-zrho_rmin)
     nbrk2=(zrho_rbrk2-zrho_rmin)*(nr+1)/(zrho_rmax-zrho_rmin)
     nbrk2=min(nbrk2,(nr-2))

     delr_hinv1=(nbrk1-1)/(rbrk1-rmin)
     delr_hinv2=(nbrk2-nbrk1)/(rbrk2-rbrk1)
     delr_hinv3=(nr-nbrk2)/(rmax-rbrk2)

     do i=1,nr
        if(i.lt.nbrk1) then
           rgrid(i)=(rmin*(nbrk1-i)+rbrk1*(i-1))/(nbrk1-1)
        else if(i.lt.nbrk2) then
           rgrid(i)=(rbrk1*(nbrk2-i)+rbrk2*(i-nbrk1))/(nbrk2-nbrk1)
        else
           rgrid(i)=(rbrk2*(nr-i)+rmax*(i-nbrk2))/(nr-nbrk2)
        endif
     enddo

     deallocate(zrho,zth,zphi,zr,zz,nregion)

     dlgrid_min=min(1/delr_hinv1,1/delr_hinv2,1/delr_hinv3,1/delz_hinv)
     dlgrid_max=max(1/delr_hinv1,1/delr_hinv2,1/delr_hinv3,1/delz_hinv)

     !
     ! some final checks:  certain xplasma things must be defined...
     !

     no_brz = .FALSE.

     call eq_gfnum('BPHI',id)
     if(id.eq.0) then
        if(ierr.eq.0) ierr=-4
        no_brz=.TRUE.
     endif

     call eq_gfnum('BR',id)
     if(id.eq.0) then
        if(ierr.eq.0) ierr=-4
        no_brz=.TRUE.
     endif

     call eq_gfnum('BZ',id)
     if(id.eq.0) then
        if(ierr.eq.0) ierr=-4
        no_brz=.TRUE.
     endif

  endif

end subroutine xstraln_init
