!----------------------------------------------------------------------
!  nblocb - given (r,y) get field and magnetic coordinates
!   inverse map used by depsub to ionize fast neutral and compute
!   guiding ctr parameters (dmc feb 1987)
!
subroutine nblocb(zr,zy,zxi,zth,zbmod,zbdir,icall,ifast,ier)
  !
  !  rga oct00 - remove levgb = 1,2
  !
  !  icall = 0 or 1 -- from depsub or cxline
  !   use inverse map to get zxi,zth
  !   if icall = 1 also
  !   return zbdir(1,2,3) = field in (toroidal,r,y) directions
  !  icall = 2 -- from nutrav/cxline
  !   jacobean known, inverse map not needed
  !   return zbdir(1,2) = bt/|b|, bp/|b|
  !   return zbdir(3,4) = flux surface tangent unit vector in poloidal plane
  !
  !  dmc may 2007: added control
  !    ifast=1, use bilinear maps
  !    ifast=0, use xplasma spline maps
  !
  !---------------------------------------
  !
  use xplasma_obj_instance
  use xstraln_calls
  use nbi_com
  use logmod
  !
  !============
  implicit none
  integer :: icall,ifast,ier,isign,id3(3)
  !============
  real(fp) :: zy,zxi,zth,zbmod,zr,zphiarg,zbr,zbz,zbphi,zbp2,zb2,zbp
  real(fp) :: zri,zzi,zdldth,zbnorm,zdth2(2)
  !============
  real(fp) :: zbdir(4)
  !
  real(fp) :: zjac(2,2)
  !
  real(fp) :: zbvec(3)
  !
  !  input (zr,zy)  location on flr of neutral being ionized
  !  output
  !        (zxi,zth)  magnetic coordinates corresponding to (zr,zy)
  !        zbmod      mod(b) at (zr,zy)
  !    if icall=1
  !        zbdir(...) normalized b components at (zr,zy)
  !               1   toroidal direction     (t)  **always positive**
  !
  !               2   major radial direction (r)  ** pll theta unit**
  !               3   vertical direction     (y)  ** vector        **
  !
  !    if icall=2
  !        zbdir as described above; similar to "btiltc" information
  !        computed in bparms
  !
  !------------------------------
  !
  !  inverse map... initial guess
  !   if near origin use full xplasma...
  !
  ier = 0
  zbmod=0.0_fp
  zbdir=0.0_fp
  !
  if(ifast.eq.1) then
     zphiarg = 0._fp
     call xs_fxget(0.01_fp,zphiarg,zr,zy,zxi,zth,ier)
     if(ier.gt.0) then
        call errorlog &
             (' nblocb(ifast=1): (r,z) to (xi,th) error, xi=',zxi)
        return
     end if

  else
     call xplasma_ctrans(s,ier, &
          r_in=zr*0.01_fp, z_in=zy*0.01_fp, &
          maptype=1, tol=1.0e-6_fp, &
          rho_out=zxi, theta_out=zth)
     if(zth.gt.xpi) zth=zth-twopi
     if(zth.lt.-xpi) zth=zth+twopi
     if(ier.gt.0) then
        call errorlog &
             (' nblocb(ifast=0): (r,z) to (xi,th) error, xi=',zxi)
        return
     end if
  end if

  if((ier.eq.0).and.(icall.gt.0)) then
     if(ifast.eq.1) then
        call xs_brz(0.01_fp,zphiarg,zr,zy,1.0e4_fp,zbr,zbz,zbphi, &
             ier)
        if(ier.gt.0) then
           call errorlog &
                (' nblocb(ifast=1): error to find (br,bz,bphi)')
           return
        end if

     else
        id3(1)=id_br
        id3(2)=id_bz
        id3(3)=id_bphi
        call xplasma_eval_prof(s,id3, &
             xplasma_r_coord,zr*0.01_fp, xplasma_z_coord,zy*0.01_fp, &
             zbvec,ier)
        zbr=zbvec(1)*1.0e4_fp   ! t->gauss
        zbz=zbvec(2)*1.0e4_fp
        zbphi=zbvec(3)*1.0e4_fp
        if(ier.gt.0) then !mg 04/22/2014
           call errorlog &
                (' nblocb(ifast=0): error to find (br,bz,bphi)')
           return
        end if
     end if

     if(ier.eq.0) then
        zbp2=zbr**2+zbz**2
        zb2=zbp2+zbphi**2
        zbmod=sqrt(zb2)
        if(zbmod.eq.0.0_fp) then
           call errorlog &
                (' nblocb: error modb = 0')
           return
        end if
        zbdir(1)=abs(zbphi)/zbmod ! |bt|/|b|
        if(icall.eq.2) then
           zbp=sqrt(zbp2)
           zbdir(2)=zbp/zbmod ! |bp|/|b|
        end if

        if(zxi.gt.0.01_fp) then
           isign=-nsnccw*nsjdotb
           if(icall.eq.1) then
              zbdir(2)=isign*zbr/zbmod ! tangent unit vector * |bp|/|b|
              zbdir(3)=isign*zbz/zbmod
           else
              zbdir(3)=isign*zbr/zbp ! tangent unit vector
              zbdir(4)=isign*zbz/zbp
           end if

        else
           id3(1)=id_r
           id3(2)=id_z
           zxi=max(1.0e-6_fp,zxi)
           call xplasma_eval_prof(s,id3(1:2), &
                xplasma_rho_coord,zxi, xplasma_theta_coord,zth, &
                zdth2,ier, ideriv2=1)
           if(ier.eq.0) then
              zdldth=sqrt(zdth2(1)*zdth2(1)+zdth2(2)*zdth2(2))
              if(icall.eq.1) then
                 zbp=sqrt(zbp2)
                 zbnorm=max(1.0e-6_fp,zbp/zbmod) ! tangent unit vector * |bp|/|b|
                 zbdir(2)=zdth2(1)*zbnorm/zdldth
                 zbdir(3)=zdth2(2)*zbnorm/zdldth
              else
                 zbdir(3)=zdth2(1)/zdldth ! tangent unit vector
                 zbdir(4)=zdth2(2)/zdldth
              end if
           else
              call errorlog &
                   (' nblocb: error to find (br,bz,bphi) at xi<=0.01')
              return
           end if
        end if
     end if
  end if
  !
  if(ier.eq.0) then
     if(abs(zth).gt.xpi+0.0001_fp) then
        write(nonlin,*) ' ?nblocb: zth not in range [-pi,pi].'
        ier=1
     else
        ier=0
     end if
  end if

  return
end subroutine nblocb
