subroutine xs_frz1_r8(id,zscalrz,zphi,zr,zz,zscalf,zf,iwarn)
 
  use xstraln_mod
  use xstraln_params
 
  implicit NONE
 
  !  bilinear lookup on indicated function (id)
 
  integer, intent(in) :: id                ! lookup function id
  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,zr,zz           ! (cyl.) coordinates
 
  REAL*8, intent(in) :: zscalf               ! factor to scale units of output
 
  REAL*8, intent(out) :: zf                  ! output item
  integer, intent(out) :: iwarn            ! out-of-range warning flag
                                           ! iwarn=-99 means "id" invalid
 
  !-----------------------------------------------------
  !  local items:
 
  REAL*8 :: buf2(2,2)       ! R,Z lookup
  REAL*8 :: buf3(2,2,2)     ! phi,R,Z lookup
 
  integer ir1,iz1,iphi1
  REAL*8 zfr,zfz,zfphi
 
  integer ir,iz,iphi
  REAL*8, parameter :: one = 1
  !-----------------------------------------------------
 
  iwarn=0
 
  call rzgrid1_r8(zphi,zscalrz*zr,zscalrz*zz, &
       iphi1,zfphi, ir1,zfr, iz1,zfz, iwarn)
 
  if(nphi.eq.1) then
 
     if(id.eq.idxs_dlim) then
        buf2(1:2,1:2)=frz_r8%dlim(1,ir1:ir1+1,iz1:iz1+1)
     else if(id.eq.idxs_dbdy) then
        buf2(1:2,1:2)=frz_r8%dbdy(1,ir1:ir1+1,iz1:iz1+1)
     else if(id.eq.idxs_modb) then
        buf2(1:2,1:2)=frz_r8%modb(1,ir1:ir1+1,iz1:iz1+1)
     else if(id.eq.idxs_binv) then
        buf2(1:2,1:2)=frz_r8%binv(1,ir1:ir1+1,iz1:iz1+1)
     else
        iwarn=-99
        zf=0
        return
     endif
 
     zf = zscalf*( &
            (one-zfz) * &
              ((one-zfr)*buf2(1,1) + &
                     zfr*buf2(2,1)) + &
             zfz * &
              ((one-zfr)*buf2(1,2) + &
                     zfr*buf2(2,2)) )
 
  else
 
     if(id.eq.idxs_dlim) then
        buf3(1:2,1:2,1:2)=frz_r8%dlim(iphi1:iphi1+1,ir1:ir1+1,iz1:iz1+1)
     else if(id.eq.idxs_dbdy) then
        buf3(1:2,1:2,1:2)=frz_r8%dbdy(iphi1:iphi1+1,ir1:ir1+1,iz1:iz1+1)
     else if(id.eq.idxs_modb) then
        buf3(1:2,1:2,1:2)=frz_r8%modb(iphi1:iphi1+1,ir1:ir1+1,iz1:iz1+1)
     else if(id.eq.idxs_binv) then
        buf3(1:2,1:2,1:2)=frz_r8%binv(iphi1:iphi1+1,ir1:ir1+1,iz1:iz1+1)
     else
        iwarn=-99
        zf=0
        return
     endif
 
     zf = zscalf*( &
         (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)) ) )
 
  endif
 
  return
 
end subroutine xs_frz1_r8
 
subroutine xs_frz_r8(id,ivec,zscalrz,zphi,zri,zzi,zscalf,zf,iwarn)
 
  use xstraln_mod
  use xstraln_params
 
  implicit NONE
 
  !  bilinear lookup on indicated function (id) -- VECTORIZED
 
  integer, intent(in) :: id                ! lookup function id
 
  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(in) :: zscalf               ! factor to scale units of output
 
  REAL*8, intent(out) :: zf(ivec)            ! output item
  integer, intent(out) :: iwarn(ivec)      ! out-of-range warning flag
                                           ! iwarn=-99 means "id" invalid
 
  !-----------------------------------------------------
  !  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
  REAL*8, parameter :: one = 1
  integer i
  !-----------------------------------------------------
 
  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)
        if(id.eq.idxs_dlim) then
           buf2(1:2,1:2)=frz_r8%dlim(1,ir:ir+1,iz:iz+1)
        else if(id.eq.idxs_dbdy) then
           buf2(1:2,1:2)=frz_r8%dbdy(1,ir:ir+1,iz:iz+1)
        else if(id.eq.idxs_modb) then
           buf2(1:2,1:2)=frz_r8%modb(1,ir:ir+1,iz:iz+1)
        else if(id.eq.idxs_binv) then
           buf2(1:2,1:2)=frz_r8%binv(1,ir:ir+1,iz:iz+1)
        else
           iwarn=-99
           zf=0
           return
        endif
 
        zf(i) = zscalf*( &
               (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)) )
     enddo
 
  else
 
     do i=1,ivec
 
        ir=ir1(i)
        iz=iz1(i)
        iphi=iphi1(i)
 
        if(id.eq.idxs_dlim) then
           buf3(1:2,1:2,1:2)= frz_r8%dlim(iphi:iphi+1,ir:ir+1,iz:iz+1)
        else if(id.eq.idxs_dbdy) then
           buf3(1:2,1:2,1:2)= frz_r8%dbdy(iphi:iphi+1,ir:ir+1,iz:iz+1)
        else if(id.eq.idxs_modb) then
           buf3(1:2,1:2,1:2)= frz_r8%modb(iphi:iphi+1,ir:ir+1,iz:iz+1)
        else if(id.eq.idxs_binv) then
           buf3(1:2,1:2,1:2)= frz_r8%binv(iphi:iphi+1,ir:ir+1,iz:iz+1)
        else
           iwarn=-99
           zf=0
           return
        endif
 
        zf = zscalf*( &
            (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)) ) )
 
     enddo
 
  endif
 
  return
 
end subroutine xs_frz_r8
