subroutine xstraln_setup_r8(ztol,ierr)
 
  use xstraln_mod
  use xplasma_calls
 
  implicit NONE
 
  REAL*8, intent(in) :: ztol     ! relative accuracy tolerance for inverse maps
 
  integer, intent(out) :: ierr     ! completion code, 0=OK
 
  !
  !  setup xstraln fields:  f(phi,theta,rho), g(phi,R,Z)
  ! *** edit "real" version & run xs_convert.csh to gen "real*8" version!
  !
  !---------------------------------
  !  modifications -- dmc Jan 2006 -- for "standalone" straight line tracker
  !  without fast inverse map support
  !
  !  if xstraln_mod'no_fRZ is set -- skip construction of all functions...
  !
  !---------------------------------
  !  local...
 
  integer irho,itheta,iphi,ir,iz,isize,iadr,iret
 
  REAL*8 zrhomin,zrhoval,zthval,zrval,zval,zb
  REAL*8, dimension(:), allocatable :: zphi,ztheta,zrho,zr,zz,zf
  REAL*8, dimension(:,:), allocatable :: bvec
  integer, dimension(:), allocatable :: nregion
  real*8 ztemp8
 
  integer irsave,izsave
 
  !  debug plotting...
 
  REAL*8, dimension(:), allocatable :: rhopl,thpl,rpl,zpl,wk
 
  !---------------------------------
 
  if( max( ffx_r8%init, frz_r8%init ) .ne. 0 ) then
     ierr=99
     return
  else
     ierr=0
  endif
 
  call alloc_r8
 
  if(no_fRZ) return
 
  !  for debug plotting...
 
  allocate(rpl(nr),zpl(nz),rhopl(nrho),thpl(ntheta))
  allocate(wk(max(nr,nz,nrho,ntheta)))
  wk=0
 
  !---------------------------------
  !  get R(phi,theta,rho) and Z(phi,theta,rho) piecewise linear maps
 
  isize=nphi*ntheta
 
  allocate(zphi(isize),ztheta(isize),zrho(isize),zr(isize),zz(isize))
  allocate(nregion(isize))
 
  zr=0
  zz=0
  nregion=0
 
  do irho=1,nrho
 
     zrhoval=rhogrid(irho)
     rhopl(irho)=zrhoval
 
     if(nphi.eq.1) then
        zphi(1)=(phimin+phimax)/2
     else
        zphi(1:nphi)=phigrid
     endif
 
     iadr=-nphi
     do itheta=1,ntheta
        iadr=iadr+nphi
 
        zthval = thgrid(itheta)
        thpl(itheta)=zthval
 
        if(iadr.gt.0) zphi(iadr+1:iadr+nphi)=zphi(1:nphi)
        zrho(iadr+1:iadr+nphi)=zrhoval
        ztheta(iadr+1:iadr+nphi)=zthval
     enddo
 
     call eq_rzget(isize,zrho,ztheta,zphi,zr,zz,nregion,iret)
     ierr=max(ierr,iret)
 
     if(ierr.eq.0) then
 
        iadr=-nphi
        do itheta=1,ntheta
           iadr=iadr+nphi
           if(irho.eq.1) iadr=0  ! suppress theta variation on axis
 
           ffx_r8%r(1:nphi,itheta,irho)=zr(iadr+1:iadr+nphi)
           ffx_r8%z(1:nphi,itheta,irho)=zz(iadr+1:iadr+nphi)
        enddo
 
     else
 
        exit    !  exit the loop, an error has occurred.
 
     endif
 
  enddo
 
  deallocate(zphi,ztheta,zrho,zr,zz)
  deallocate(nregion)
 
  if(ierr.ne.0) go to 9999
 
  !---------------------------------
  ! compute all f(phi,R,Z) fast maps
 
  isize=nphi*nr
 
  allocate(zphi(isize),ztheta(isize),zrho(isize),zr(isize),zz(isize))
  allocate(zf(isize),bvec(3,isize))
  allocate(nregion(isize))
 
  zrho=0
  ztheta=0
  zf=0
  bvec=0
  nregion=0
 
  do iz=1,nz
 
     zval = zgrid(iz)
     zpl(iz)=zval
 
     if(nphi.eq.1) then
        zphi(1)=(phimin+phimax)/2
     else
        zphi=phigrid
     endif
 
     iadr=-nphi
     do ir=1,nr
        iadr=iadr+nphi
 
        zrval=rgrid(ir)
        rpl(ir)=zrval
 
        if(iadr.gt.0) zphi(iadr+1:iadr+nphi)=zphi(1:nphi)
        zr(iadr+1:iadr+nphi)=zrval
        zz(iadr+1:iadr+nphi)=zval
     enddo
 
     !  inverse map...
 
     call eq_inv(isize,zr,zz,zphi,ztol,zrho,ztheta,nregion,iret)
     ierr=max(ierr,iret)
 
     if(ierr.eq.0) then
 
        iadr=-nphi
        do ir=1,nr
           iadr=iadr+nphi
 
           frz_r8%rho(1:nphi,ir,iz)=zrho(iadr+1:iadr+nphi)
           do iphi=1,nphi
              ztemp8=ztheta(iadr+iphi)
              do while(ztemp8 < -xpi)
                 ztemp8 = ztemp8 + twopi
              enddo
              do while(ztemp8 > xpi)
                 ztemp8 = ztemp8 - twopi
              enddo
              frz_r8%theta(iphi,ir,iz)=ztemp8
           enddo
        enddo
 
     else
 
        exit    !  exit the loop, an error has occurred.
 
     endif
 
     !  additional functions...
 
     !  distance to plasma boundary...
 
     call eq_dfast(isize,zr,zz,zphi,zf,iret)
     ierr=max(ierr,iret)
 
     if(ierr.eq.0) then
 
        iadr=-nphi
        do ir=1,nr
           iadr=iadr+nphi
 
           frz_r8%dbdy(1:nphi,ir,iz)=zf(iadr+1:iadr+nphi)
        enddo
 
     else
 
        exit    !  exit the loop, an error has occurred.
 
     endif
 
     !  distance to limiter...
 
     call eq_rzdist(isize,zr,zz,zphi,zf,0,zrho,zrho,zrho,iret)
     ierr=max(ierr,iret)
 
     if(ierr.eq.0) then
 
        iadr=-nphi
        do ir=1,nr
           iadr=iadr+nphi
 
           frz_r8%dlim(1:nphi,ir,iz)=zf(iadr+1:iadr+nphi)
        enddo
 
     else
 
        exit    !  exit the loop, an error has occurred.
 
     endif
 
     !  field strength
 
     if(.not.NO_BRZ) then
        call eqrz_bget(isize,zr,zz,zphi,ztol,bvec,nregion,iret)
        ierr=max(ierr,iret)
 
        if(ierr.eq.0) then
 
           iadr=-nphi
           do ir=1,nr
              iadr=iadr+nphi
 
              do iphi=1,nphi
                 zb=sqrt(bvec(1,iadr+iphi)**2 + &
                      bvec(2,iadr+iphi)**2 + &
                      bvec(3,iadr+iphi)**2)
                 if(zb.eq.0) then
                    zb=1
                    bvec(3,iadr+iphi)=1
                 endif
                 zf(iphi)=zb
              enddo
 
              frz_r8%modb(1:nphi,ir,iz)=zf(1:nphi)
 
              zf(1:nphi)=1/zf(1:nphi)
              frz_r8%binv(1:nphi,ir,iz)=zf(1:nphi)
 
              frz_r8%br(1:nphi,ir,iz)=bvec(1,iadr+1:iadr+nphi)
              frz_r8%bz(1:nphi,ir,iz)=bvec(2,iadr+1:iadr+nphi)
              frz_r8%bphi(1:nphi,ir,iz)=bvec(3,iadr+1:iadr+nphi)
 
           enddo
 
        else
 
           exit    !  exit the loop, an error has occurred.
 
        endif
     endif
 
  enddo
 
  !  compute some stuff related to handling inv. maps near the axis
 
  invmap_tol = ztol
 
  rho_small=0
 
  do iphi=1,nphi
 
     zrhomin=rhomax
 
     do iz=1,nz
        do ir=1,nr
 
           if(frz_r8%rho(iphi,ir,iz).lt.zrhomin) then
              zrhomin=frz_r8%rho(iphi,ir,iz)
              irsave=ir
              izsave=iz
           endif
 
        enddo
     enddo
 
     ir=irsave
     iz=izsave
 
     ! take largest rho 1.5 zone widths from minimum for this phi
 
     ztemp8=(frz_r8%rho(iphi,ir-2,iz)+frz_r8%rho(iphi,ir-1,iz))/2
     rho_small=max(rho_small,ztemp8)
     ztemp8=(frz_r8%rho(iphi,ir+2,iz)+frz_r8%rho(iphi,ir+1,iz))/2
     rho_small=max(rho_small,ztemp8)
     ztemp8=(frz_r8%rho(iphi,ir,iz-2)+frz_r8%rho(iphi,ir,iz-1))/2
     rho_small=max(rho_small,ztemp8)
     ztemp8=(frz_r8%rho(iphi,ir,iz+2)+frz_r8%rho(iphi,ir,iz+1))/2
     rho_small=max(rho_small,ztemp8)
 
  enddo
 
  deallocate(zphi,ztheta,zrho,zr,zz)
  deallocate(zf,bvec)
  deallocate(nregion)
 
!  if(ierr.eq.0) call copy_r8
 
9999 continue
 
  !  debug plotting... assumes nphi=1 & real*4 (for now)
 
  !dbg call grf3f1(thpl,ffx_r4%r,rhopl,wk,ntheta,nrho,ntheta,  &
  !dbg      'theta','R','rho',' ','m',' ','xstraln_setup','debug plot',0)
 
  !dbg call grf3f1(thpl,ffx_r4%z,rhopl,wk,ntheta,nrho,ntheta,  &
  !dbg      'theta','Z','rho',' ','m',' ','xstraln_setup','debug plot',0)
 
  !dbg call grf3f1(rpl,frz_r4%rho,zpl,wk,nr,nz,nr, &
  !dbg      'R','rho','Z','m',' ','m','xstraln_setup','debug plot',0)
 
  !dbg call grf3f1(rpl,frz_r4%theta,zpl,wk,nr,nz,nr, &
  !dbg      'R','theta','Z','m','rad','m','xstraln_setup','debug plot',0)
 
  !dbg call grf3f1(rpl,frz_r4%dbdy,zpl,wk,nr,nz,nr, &
  !dbg      'R','dbdy','Z','m','m','m','xstraln_setup','debug plot',0)
 
  !dbg call grf3f1(rpl,frz_r4%dlim,zpl,wk,nr,nz,nr, &
  !dbg      'R','dlim','Z','m','m','m','xstraln_setup','debug plot',0)
 
  !dbg call grf3f1(rpl,frz_r4%modb,zpl,wk,nr,nz,nr, &
  !dbg      'R','mod(B)','Z','m','T','m','xstraln_setup','debug plot',0)
 
  !dbg call grf3f1(rpl,frz_r4%binv,zpl,wk,nr,nz,nr, &
  !dbg      'R','1/mod(B)','Z','m','T^-1','m','xstraln_setup','debug plot',0)
 
  deallocate(rpl,zpl,rhopl,thpl,wk)
  return
 
end subroutine xstraln_setup_r8
