subroutine xs_brz1_r8(zscalrz,zphi,zr,zz,zscalb,br,bz,bphi,iwarn)
 
  use xstraln_mod
  use xstraln_params
 
  implicit NONE
 
  !  bilinear lookup on indicated 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 (m)*zscalrz
 
  REAL*8, intent(in) :: zscalb               ! factor to scale units of output
 
  REAL*8, intent(out) :: br,bz,bphi          ! output [BR,BZ,Bphi] (T)/zscalb
 
  integer, intent(out) :: iwarn            ! out-of-range warning flag
                                           ! iwarn=-99 means "id" invalid
                                           ! iwarn=-100 means NO_BRZ = .TRUE.
 
  !-----------------------------------------------------
  !  local items:
 
  integer ir1,iz1,iphi1
  REAL*8 zfr,zfz,zfphi
 
  integer ir,iz,iphi
  REAL*8, parameter :: one = 1
  !-----------------------------------------------------
 
  iwarn=0
  if(NO_BRZ) then
     br=0
     bz=0
     bphi=0
     iwarn=-100
     return
  endif
 
  call rzgrid1_r8(zphi,zscalrz*zr,zscalrz*zz, &
       iphi1,zfphi, ir1,zfr, iz1,zfz, iwarn)
 
  if(nphi.eq.1) then
 
     br = zscalb*( &
            (one-zfz) * &
              ((one-zfr)*frz_r8%br(1,ir1,iz1) + &
                     zfr*frz_r8%br(1,ir1+1,iz1)) + &
             zfz * &
              ((one-zfr)*frz_r8%br(1,ir1,iz1+1) + &
                     zfr*frz_r8%br(1,ir1+1,iz1+1)) )
 
     bz = zscalb*( &
            (one-zfz) * &
              ((one-zfr)*frz_r8%bz(1,ir1,iz1) + &
                     zfr*frz_r8%bz(1,ir1+1,iz1)) + &
             zfz * &
              ((one-zfr)*frz_r8%bz(1,ir1,iz1+1) + &
                     zfr*frz_r8%bz(1,ir1+1,iz1+1)) )
 
     bphi = zscalb*( &
            (one-zfz) * &
              ((one-zfr)*frz_r8%bphi(1,ir1,iz1) + &
                     zfr*frz_r8%bphi(1,ir1+1,iz1)) + &
             zfz * &
              ((one-zfr)*frz_r8%bphi(1,ir1,iz1+1) + &
                     zfr*frz_r8%bphi(1,ir1+1,iz1+1)) )
 
  else
 
     br = zscalb*( &
         (one-zfz) * &
            ((one-zfr) * &
              ((one-zfphi)*frz_r8%br(iphi1,ir1,iz1) + &
                     zfphi*frz_r8%br(iphi1+1,ir1,iz1)) + &
             zfr * &
              ((one-zfphi)*frz_r8%br(iphi1,ir1+1,iz1) + &
                     zfphi*frz_r8%br(iphi1+1,ir1+1,iz1)) ) + &
          zfz * &
            ((one-zfr) * &
              ((one-zfphi)*frz_r8%br(iphi1,ir1,iz1+1) + &
                     zfphi*frz_r8%br(iphi1+1,ir1,iz1+1)) + &
             zfr * &
              ((one-zfphi)*frz_r8%br(iphi1,ir1+1,iz1+1) + &
                     zfphi*frz_r8%br(iphi1+1,ir1+1,iz1+1)) ) )
 
     bz = zscalb*( &
         (one-zfz) * &
            ((one-zfr) * &
              ((one-zfphi)*frz_r8%bz(iphi1,ir1,iz1) + &
                     zfphi*frz_r8%bz(iphi1+1,ir1,iz1)) + &
             zfr * &
              ((one-zfphi)*frz_r8%bz(iphi1,ir1+1,iz1) + &
                     zfphi*frz_r8%bz(iphi1+1,ir1+1,iz1)) ) + &
          zfz * &
            ((one-zfr) * &
              ((one-zfphi)*frz_r8%bz(iphi1,ir1,iz1+1) + &
                     zfphi*frz_r8%bz(iphi1+1,ir1,iz1+1)) + &
             zfr * &
              ((one-zfphi)*frz_r8%bz(iphi1,ir1+1,iz1+1) + &
                     zfphi*frz_r8%bz(iphi1+1,ir1+1,iz1+1)) ) )
 
     bphi = zscalb*( &
         (one-zfz) * &
            ((one-zfr) * &
              ((one-zfphi)*frz_r8%bphi(iphi1,ir1,iz1) + &
                     zfphi*frz_r8%bphi(iphi1+1,ir1,iz1)) + &
             zfr * &
              ((one-zfphi)*frz_r8%bphi(iphi1,ir1+1,iz1) + &
                     zfphi*frz_r8%bphi(iphi1+1,ir1+1,iz1)) ) + &
          zfz * &
            ((one-zfr) * &
              ((one-zfphi)*frz_r8%bphi(iphi1,ir1,iz1+1) + &
                     zfphi*frz_r8%bphi(iphi1+1,ir1,iz1+1)) + &
             zfr * &
              ((one-zfphi)*frz_r8%bphi(iphi1,ir1+1,iz1+1) + &
                     zfphi*frz_r8%bphi(iphi1+1,ir1+1,iz1+1)) ) )
 
  endif
 
  return
 
end subroutine xs_brz1_r8
 
subroutine xs_brz_r8(ivec,zscalrz,zphi,zri,zzi,zscalb,br,bz,bphi,iwarn)
 
  use xstraln_mod
  use xstraln_params
 
  implicit NONE
 
  !  bilinear lookup on indicated function (id) -- 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(in) :: zscalb               ! factor to scale units of output
 
  REAL*8, intent(out) :: br(ivec),bz(ivec),bphi(ivec) ! output [BR,BZ,Bphi]
                                           ! (T)/zscalb
 
  integer, intent(out) :: iwarn(ivec)      ! out-of-range warning flag
                                           ! iwarn=-99 means "id" invalid
 
  !-----------------------------------------------------
  !  local items:
 
  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)
 
        br(i) = zscalb*( &
            (one-zfz(i)) * &
              ((one-zfr(i))*frz_r8%br(1,ir,iz) + &
                     zfr(i)*frz_r8%br(1,ir+1,iz)) + &
             zfz(i) * &
              ((one-zfr(i))*frz_r8%br(1,ir,iz+1) + &
                     zfr(i)*frz_r8%br(1,ir+1,iz+1)) )
 
        bz(i) = zscalb*( &
            (one-zfz(i)) * &
              ((one-zfr(i))*frz_r8%bz(1,ir,iz) + &
                     zfr(i)*frz_r8%bz(1,ir+1,iz)) + &
             zfz(i) * &
              ((one-zfr(i))*frz_r8%bz(1,ir,iz+1) + &
                     zfr(i)*frz_r8%bz(1,ir+1,iz+1)) )
 
        bphi(i) = zscalb*( &
            (one-zfz(i)) * &
              ((one-zfr(i))*frz_r8%bphi(1,ir,iz) + &
                     zfr(i)*frz_r8%bphi(1,ir+1,iz)) + &
             zfz(i) * &
              ((one-zfr(i))*frz_r8%bphi(1,ir,iz+1) + &
                     zfr(i)*frz_r8%bphi(1,ir+1,iz+1)) )
 
     enddo
 
  else
 
     do i=1,ivec
 
        ir=ir1(i)
        iz=iz1(i)
        iphi=iphi1(i)
 
        br(i) = zscalb*( &
         (one-zfz(i)) * &
            ((one-zfr(i)) * &
              ((one-zfphi(i))*frz_r8%br(iphi,ir,iz) + &
                     zfphi(i)*frz_r8%br(iphi+1,ir,iz)) + &
             zfr(i) * &
              ((one-zfphi(i))*frz_r8%br(iphi,ir+1,iz) + &
                     zfphi(i)*frz_r8%br(iphi+1,ir+1,iz)) ) + &
          zfz(i) * &
            ((one-zfr(i)) * &
              ((one-zfphi(i))*frz_r8%br(iphi,ir,iz+1) + &
                     zfphi(i)*frz_r8%br(iphi+1,ir,iz+1)) + &
             zfr(i) * &
              ((one-zfphi(i))*frz_r8%br(iphi,ir+1,iz+1) + &
                     zfphi(i)*frz_r8%br(iphi+1,ir+1,iz+1)) ) )
 
        bz(i) = zscalb*( &
         (one-zfz(i)) * &
            ((one-zfr(i)) * &
              ((one-zfphi(i))*frz_r8%bz(iphi,ir,iz) + &
                     zfphi(i)*frz_r8%bz(iphi+1,ir,iz)) + &
             zfr(i) * &
              ((one-zfphi(i))*frz_r8%bz(iphi,ir+1,iz) + &
                     zfphi(i)*frz_r8%bz(iphi+1,ir+1,iz)) ) + &
          zfz(i) * &
            ((one-zfr(i)) * &
              ((one-zfphi(i))*frz_r8%bz(iphi,ir,iz+1) + &
                     zfphi(i)*frz_r8%bz(iphi+1,ir,iz+1)) + &
             zfr(i) * &
              ((one-zfphi(i))*frz_r8%bz(iphi,ir+1,iz+1) + &
                     zfphi(i)*frz_r8%bz(iphi+1,ir+1,iz+1)) ) )
 
        bphi(i) = zscalb*( &
         (one-zfz(i)) * &
            ((one-zfr(i)) * &
              ((one-zfphi(i))*frz_r8%bphi(iphi,ir,iz) + &
                     zfphi(i)*frz_r8%bphi(iphi+1,ir,iz)) + &
             zfr(i) * &
              ((one-zfphi(i))*frz_r8%bphi(iphi,ir+1,iz) + &
                     zfphi(i)*frz_r8%bphi(iphi+1,ir+1,iz)) ) + &
          zfz(i) * &
            ((one-zfr(i)) * &
              ((one-zfphi(i))*frz_r8%bphi(iphi,ir,iz+1) + &
                     zfphi(i)*frz_r8%bphi(iphi+1,ir,iz+1)) + &
             zfr(i) * &
              ((one-zfphi(i))*frz_r8%bphi(iphi,ir+1,iz+1) + &
                     zfphi(i)*frz_r8%bphi(iphi+1,ir+1,iz+1)) ) )
 
     enddo
 
  endif
 
  return
 
end subroutine xs_brz_r8
