subroutine xs_fxget1_r8(zscalrz,zphi,zri,zzi,zrho,zth,iwarn)
 
  use xstraln_mod
  use xstraln_params
 
  implicit NONE
 
  !  (phi,R,Z) -> ([phi],rho,theta) fast linear map
 
  !  mod DMC: near axis, make continuous transition to accurate map, instead
  !  of a discontinuous transition at rho_small
  !    -- transition btw 1*rho_small and (1/2)*rho_small
  !       decided against (3/2)*rho_small on the high end after tests,
  !       due to cpu cost of eq_inv calls...
 
  REAL*8, intent(in) :: zscalrz              ! factor to convert (R,Z) -> meters
                                           ! use 0.01 if R,Z given in cm
 
  REAL*8, intent(in) :: zphi,zri,zzi         ! (cyl.) coordinates
 
  REAL*8, intent(out) :: zrho,zth            ! (mag.) coordinates returned
  integer, intent(out) :: iwarn            ! out-of-range warning flag
 
  !-----------------------------------------------------
  !  local items:
 
  REAL*8 :: buf2(2,2)       ! R,Z lookup
  REAL*8 :: buf3(2,2,2)     ! phi,R,Z lookup
 
  REAL*8 zr,zz  ! scaled
 
  integer ir1,iz1,iphi1
  REAL*8 zfr,zfz,zfphi
 
  integer ir,iz,iphi,iregion,ierr,i,j,k,iadj
  REAL*8, parameter :: zero = 0
  REAL*8, parameter :: one = 1
  REAL*8 :: ztol,zthmin,zthmax,zxblend,zblend,zrhoinv,zthinv
 
  integer, save :: loc_init = 0
  REAL*8, save :: half, z3ov2 ! 0.5,1.5
 
  !-----------------------------------------------------
  if(loc_init.eq.0) then
     loc_init=1
     half = one/2
     z3ov2 = one + half
  endif
 
  iwarn=0
 
  zr=zscalrz*zri
  zz=zscalrz*zzi
 
  call rzgrid1_r8(zphi,zr,zz, &
       iphi1,zfphi, ir1,zfr, iz1,zfz, iwarn)
 
  if(nphi.eq.1) then
 
     ! bilinear rho
 
     buf2(1:2,1:2)=frz_r8%rho(1,ir1:ir1+1,iz1:iz1+1)
 
     zrho = ( &
             (one-zfz) * &
               ((one-zfr)*buf2(1,1) + &
                      zfr*buf2(2,1)) + &
              zfz * &
               ((one-zfr)*buf2(1,2) + &
                      zfr*buf2(2,2)) )
 
     ! bilinear theta
 
     buf2(1:2,1:2)=frz_r8%theta(1,ir1:ir1+1,iz1:iz1+1)
 
     iadj=0
     zthmin=min(buf2(1,1),buf2(1,2),buf2(2,1),buf2(2,2))
     zthmax=max(buf2(1,1),buf2(1,2),buf2(2,1),buf2(2,2))
 
     if(2*(zthmax-zthmin).gt.th_period) then
 
        ! deal with branchcut
 
        do j=1,2
           do i=1,2
              if(2*(zthmax-buf2(i,j)).gt.th_period) then
                 iadj=iadj+1
                 buf2(i,j)=buf2(i,j)+th_period
              endif
           enddo
        enddo
 
     endif
 
     zth = ( &
         (one-zfz) * &
           ((one-zfr)*buf2(1,1) + &
                  zfr*buf2(2,1)) + &
          zfz * &
           ((one-zfr)*buf2(1,2) + &
                  zfr*buf2(2,2)) )
 
     if(iadj.gt.0) then
        if(zth.gt.thmax) zth=zth-th_period
     endif
 
     if(zrho.lt.one*rho_small) then
 
        ! revert to xplasma iterative map if too close to axis
 
        if(zrho.le.half*rho_small) then
           zxblend= one
           zblend = one
        else
           zxblend = min(one,max(zero, 2*(one - zrho/rho_small) ))
           zblend = zxblend*zxblend*(3 - 2*zxblend)
           ! this a cubic satisfying f(0)=f'(0)=f'(1)=0 and f(1)=1
        endif
 
        ztol = invmap_tol
        call eq_inv(1,zr,zz,zphi,ztol,zrhoinv,zthinv,iregion,ierr)
        if(ierr.ne.0) then
           zrhoinv=rhomin
           zthinv=(thmin+thmax)/2
        endif
 
        if(zblend.ge.one) zth=zero
 
        ! if blending, bring zthinv into same branch as zth
 
        do while ( zthinv < zth-xpi )
           zthinv = zthinv + twopi
        enddo
        do while ( zthinv > zth+xpi )
           zthinv = zthinv - twopi
        enddo
 
        if(zblend.ge.one) then
           zrho=zrhoinv
           zth=zthinv
        else
           zrho=(one-zblend)*zrho + zblend*zrhoinv
           zth =(one-zblend)*zth  + zblend*zthinv
 
           do while ( zth < -xpi )
              zth = zth + twopi
           enddo
           do while ( zth > xpi )
              zth = zth - twopi
           enddo
        endif
 
     endif   ! (zrho.ge.rho_small)
 
  else
 
     ! trilinear rho
 
     buf3(1:2,1:2,1:2)=frz_r8%rho(iphi1:iphi1+1,ir1:ir1+1,iz1:iz1+1)
 
     zrho =( &
            (one-zfz) * &
               ((one-zfr) * &
                 ((one-zfphi)*buf3(1,1,1) + &
                        zfphi*buf3(2,1,1)) + &
                      zfr * &
                 ((one-zfphi)*buf3(1,2,1) + &
                        zfphi*buf3(2,2,1)) ) + &
                  zfz * &
               ((one-zfr) * &
                 ((one-zfphi)*buf3(1,1,2) + &
                        zfphi*buf3(2,1,2)) + &
                      zfr * &
                 ((one-zfphi)*buf3(1,2,2) + &
                        zfphi*buf3(2,2,2)) ) )
 
     ! trilinear theta
 
     buf3(1:2,1:2,1:2)=frz_r8%theta(iphi1:iphi1+1,ir1:ir1+1,iz1:iz1+1)
 
     iadj=0
     zthmin=min(buf3(1,1,1),buf3(1,1,2),buf3(1,2,1),buf3(1,2,2), &
          buf3(2,1,1),buf3(2,1,2),buf3(2,2,1),buf3(2,2,2))
     zthmax=max(buf3(1,1,1),buf3(1,1,2),buf3(1,2,1),buf3(1,2,2), &
          buf3(2,1,1),buf3(2,1,2),buf3(2,2,1),buf3(2,2,2))
 
     if(2*(zthmax-zthmin).gt.th_period) then
 
        ! deal with branchcut
 
        do k=1,2
           do j=1,2
              do i=1,2
                 if(2*(zthmax-buf3(i,j,k)).gt.th_period) then
                    iadj=iadj+1
                    buf3(i,j,k)=buf3(i,j,k)+th_period
                 endif
              enddo
           enddo
        enddo
 
     endif
 
     zth =( &
            (one-zfz) * &
               ((one-zfr) * &
                 ((one-zfphi)*buf3(1,1,1) + &
                        zfphi*buf3(2,1,1)) + &
                      zfr * &
                 ((one-zfphi)*buf3(1,2,1) + &
                        zfphi*buf3(2,2,1)) ) + &
                  zfz * &
               ((one-zfr) * &
                 ((one-zfphi)*buf3(1,1,2) + &
                        zfphi*buf3(2,1,2)) + &
                      zfr * &
                 ((one-zfphi)*buf3(1,2,2) + &
                        zfphi*buf3(2,2,2)) ) )
 
     if(iadj.gt.0) then
        if(zth.gt.thmax) zth=zth-th_period
     endif
 
     if(zrho.lt.one*rho_small) then
 
        ! revert to xplasma iterative map if too close to axis
 
        if(zrho.le.half*rho_small) then
           zxblend= one
           zblend = one
        else
           zxblend = min(one,max(zero, 2*(one - zrho/rho_small) ))
           zblend = zxblend*zxblend*(3 - 2*zxblend)
           ! this a cubic satisfying f(0)=f'(0)=f'(1)=0 and f(1)=1
        endif
 
        ztol = invmap_tol
        call eq_inv(1,zr,zz,zphi,ztol,zrhoinv,zthinv,iregion,ierr)
        if(ierr.ne.0) then
           zrhoinv=rhomin
           zthinv=(thmin+thmax)/2
        endif
 
        if(zblend.ge.one) zth=zero
 
        ! if blending, bring zthinv into same branch as zth
 
        do while ( zthinv < zth-xpi )
           zthinv = zthinv + twopi
        enddo
        do while ( zthinv > zth+xpi )
           zthinv = zthinv - twopi
        enddo
 
        if(zblend.ge.one) then
           zrho=zrhoinv
           zth=zthinv
        else
           zrho=(one-zblend)*zrho + zblend*zrhoinv
           zth =(one-zblend)*zth  + zblend*zthinv
 
           do while ( zth < -xpi )
              zth = zth + twopi
           enddo
           do while ( zth > xpi )
              zth = zth - twopi
           enddo
        endif
 
     endif   ! (zrho.ge.rho_small)
 
  endif   ! (nphi.eq.1 test)
 
  return
 
end subroutine xs_fxget1_r8
 
subroutine xs_fxget_r8(ivec,zscalrz,zphi,zri,zzi,zrho,zth,iwarn)
 
  !  mod DMC: near axis, make continuous transition to accurate map, instead
  !  of a discontinuous transition at rho_small
  !    -- transition btw 1*rho_small and (1/2)*rho_small
  !       decided against (3/2)*rho_small on the high end after tests,
  !       due to cpu cost of eq_inv calls...
 
  use xstraln_mod
  use xstraln_params
 
  implicit NONE
 
  !  (phi,R,Z) -> ([phi],rho,theta) fast linear map  -- VECTORIZED
 
  integer, intent(in) :: ivec              ! vector dimension
  REAL*8, intent(in) :: zscalrz              ! factor to convert (R,Z) -> meters
                                           ! use 0.01 if R,Z given in cm
 
  REAL*8, intent(in) :: zphi(ivec),zri(ivec),zzi(ivec)    ! (cyl.) coordinates
 
  REAL*8, intent(out) :: zrho(ivec),zth(ivec)   ! (mag.) coordinates returned
  integer, intent(out) :: iwarn(ivec)      ! out-of-range warning flag
 
  !-----------------------------------------------------
  !  local items:
 
  REAL*8 :: buf2(2,2)       ! R,Z lookup
  REAL*8 :: buf3(2,2,2)     ! phi,R,Z lookup
 
  REAL*8 zr(ivec),zz(ivec)  ! R,Z scaled...
 
  integer ir1(ivec),iz1(ivec),iphi1(ivec)
  REAL*8 zfr(ivec),zfz(ivec),zfphi(ivec)
 
  integer ir,iz,iphi,ii,j,k,iadj,iregion,ierr
  REAL*8, parameter :: zero = 0
  REAL*8, parameter :: one = 1
  REAL*8 :: ztol,zthmin,zthmax,zxblend,zblend,zrhoinv,zthinv
  integer i
 
  integer, save :: loc_init = 0
  REAL*8, save :: half, z3ov2 ! 0.5,1.5
 
  !-----------------------------------------------------
  if(loc_init.eq.0) then
     loc_init=1
     half = one/2
     z3ov2 = one + half
  endif
 
  iwarn=0
 
  zr=zri*zscalrz
  zz=zzi*zscalrz
 
  call rzgrid_r8(ivec,zphi,zr,zz, &
       iphi1,zfphi, ir1,zfr, iz1,zfz, iwarn)
 
  if(nphi.eq.1) then
 
     do i=1,ivec
 
        ir=ir1(i)
        iz=iz1(i)
 
        ! bilinear rho
 
        buf2(1:2,1:2)=frz_r8%rho(1,ir:ir+1,iz:iz+1)
 
        zrho(i) = ( &
                (one-zfz(i)) * &
                  ((one-zfr(i))*buf2(1,1) + &
                         zfr(i)*buf2(2,1)) + &
                      zfz(i) * &
                  ((one-zfr(i))*buf2(1,2) + &
                         zfr(i)*buf2(2,2)) )
 
        ! bilinear theta
 
        buf2(1:2,1:2)=frz_r8%theta(1,ir:ir+1,iz:iz+1)
 
        iadj=0
        zthmin=min(buf2(1,1),buf2(1,2),buf2(2,1),buf2(2,2))
        zthmax=max(buf2(1,1),buf2(1,2),buf2(2,1),buf2(2,2))
 
        if(2*(zthmax-zthmin).gt.th_period) then
 
           ! deal with branchcut
 
           do j=1,2
              do ii=1,2
                 if(2*(zthmax-buf2(ii,j)).gt.th_period) then
                    iadj=iadj+1
                    buf2(ii,j)=buf2(ii,j)+th_period
                 endif
              enddo
           enddo
 
        endif
 
        zth(i) = ( &
               (one-zfz(i)) * &
                 ((one-zfr(i))*buf2(1,1) + &
                        zfr(i)*buf2(2,1)) + &
                zfz(i) * &
                 ((one-zfr(i))*buf2(1,2) + &
                        zfr(i)*buf2(2,2)) )
 
        if(iadj.gt.0) then
           if(zth(i).gt.thmax) zth(i)=zth(i)-th_period
        endif
 
        if(zrho(i).lt.one*rho_small) then
 
           ! revert to xplasma iterative map if too close to axis
 
           if(zrho(i).le.half*rho_small) then
              zxblend= one
              zblend = one
           else
              zxblend = min(one,max(zero, 2*(one - zrho(i)/rho_small) ))
              zblend = zxblend*zxblend*(3 - 2*zxblend)
              ! this a cubic satisfying f(0)=f'(0)=f'(1)=0 and f(1)=1
           endif
 
           ztol = invmap_tol
           call eq_inv(1,zr(i),zz(i),zphi(i),ztol,zrhoinv,zthinv, &
                iregion,ierr)
           if(ierr.ne.0) then
              zrhoinv=rhomin
              zthinv=(thmin+thmax)/2
           endif
 
           if(zblend.ge.one) zth(i)=zero
 
           ! if blending, bring zthinv into same branch as zth(i)
 
           do while ( zthinv < zth(i)-xpi )
              zthinv = zthinv + twopi
           enddo
           do while ( zthinv > zth(i)+xpi )
              zthinv = zthinv - twopi
           enddo
 
           if(zblend.ge.one) then
              zrho(i)=zrhoinv
              zth(i)=zthinv
           else
              zrho(i)=(one-zblend)*zrho(i) + zblend*zrhoinv
              zth(i) =(one-zblend)*zth(i)  + zblend*zthinv
 
              do while ( zth(i) < -xpi )
                 zth(i) = zth(i) + twopi
              enddo
              do while ( zth(i) > xpi )
                 zth(i) = zth(i) - twopi
              enddo
           endif
 
        endif   ! (zrho.ge.rho_small)
 
     enddo
 
  else
 
     do i=1,ivec
 
        ir=ir1(i)
        iz=iz1(i)
        iphi=iphi1(i)
 
        ! trilinear rho
 
        buf3(1:2,1:2,1:2)=frz_r8%rho(iphi:iphi+1,ir:ir+1,iz:iz+1)
 
        zrho(i) =( &
            (one-zfz(i)) * &
               ((one-zfr(i)) * &
                 ((one-zfphi(i))*buf3(1,1,1) + &
                        zfphi(i)*buf3(2,1,1)) + &
                      zfr(i) * &
                 ((one-zfphi(i))*buf3(1,2,1) + &
                        zfphi(i)*buf3(2,2,1)) ) + &
                  zfz(i) * &
               ((one-zfr(i)) * &
                 ((one-zfphi(i))*buf3(1,1,2) + &
                        zfphi(i)*buf3(2,1,2)) + &
                      zfr(i) * &
                 ((one-zfphi(i))*buf3(1,2,2) + &
                        zfphi(i)*buf3(2,2,2)) ) )
 
        ! trilinear theta
 
        buf3(1:2,1:2,1:2)=frz_r8%theta(iphi:iphi+1,ir:ir+1,iz:iz+1)
 
        iadj=0
        zthmin=min(buf3(1,1,1),buf3(1,1,2),buf3(1,2,1),buf3(1,2,2), &
             buf3(2,1,1),buf3(2,1,2),buf3(2,2,1),buf3(2,2,2))
        zthmax=max(buf3(1,1,1),buf3(1,1,2),buf3(1,2,1),buf3(1,2,2), &
             buf3(2,1,1),buf3(2,1,2),buf3(2,2,1),buf3(2,2,2))
 
        if(2*(zthmax-zthmin).gt.th_period) then
 
           ! deal with branchcut
 
           do k=1,2
              do j=1,2
                 do ii=1,2
                    if(2*(zthmax-buf3(ii,j,k)).gt.th_period) then
                       iadj=iadj+1
                       buf3(ii,j,k)=buf3(ii,j,k)+th_period
                    endif
                 enddo
              enddo
           enddo
 
        endif
 
        zth(i) =( &
            (one-zfz(i)) * &
               ((one-zfr(i)) * &
                 ((one-zfphi(i))*buf3(1,1,1) + &
                        zfphi(i)*buf3(2,1,1)) + &
                      zfr(i) * &
                 ((one-zfphi(i))*buf3(1,2,1) + &
                        zfphi(i)*buf3(2,2,1)) ) + &
                  zfz(i) * &
               ((one-zfr(i)) * &
                 ((one-zfphi(i))*buf3(1,1,2) + &
                        zfphi(i)*buf3(2,1,2)) + &
                      zfr(i) * &
                 ((one-zfphi(i))*buf3(1,2,2) + &
                        zfphi(i)*buf3(2,2,2)) ) )
 
        if(iadj.gt.0) then
           if(zth(i).gt.thmax) zth(i)=zth(i)-th_period
        endif
 
        if(zrho(i).lt.one*rho_small) then
 
           ! revert to xplasma iterative map if too close to axis
 
           if(zrho(i).le.half*rho_small) then
              zxblend= one
              zblend = one
           else
              zxblend = min(one,max(zero, 2*(one - zrho(i)/rho_small) ))
              zblend = zxblend*zxblend*(3 - 2*zxblend)
              ! this a cubic satisfying f(0)=f'(0)=f'(1)=0 and f(1)=1
           endif
 
           ztol = invmap_tol
           call eq_inv(1,zr(i),zz(i),zphi(i),ztol,zrhoinv,zthinv, &
                iregion,ierr)
           if(ierr.ne.0) then
              zrhoinv=rhomin
              zthinv=(thmin+thmax)/2
           endif
 
           if(zblend.ge.one) zth(i)=zero
 
           ! if blending, bring zthinv into same branch as zth(i)
 
           do while ( zthinv < zth(i)-xpi )
              zthinv = zthinv + twopi
           enddo
           do while ( zthinv > zth(i)+xpi )
              zthinv = zthinv - twopi
           enddo
 
           if(zblend.ge.one) then
              zrho(i)=zrhoinv
              zth(i)=zthinv
           else
              zrho(i)=(one-zblend)*zrho(i) + zblend*zrhoinv
              zth(i) =(one-zblend)*zth(i)  + zblend*zthinv
 
              do while ( zth(i) < -xpi )
                 zth(i) = zth(i) + twopi
              enddo
              do while ( zth(i) > xpi )
                 zth(i) = zth(i) - twopi
              enddo
           endif
 
        endif   ! (zrho.ge.rho_small)
 
     enddo
 
  endif
 
  return
 
end subroutine xs_fxget_r8
