subroutine gintrp(icall)
! mg, 12/19/2013 , i restore (-nsnccw) switch as it is important 
! for promt losses due to [b x grad b] drift,
! i accomplished this by assigning a sign to q while using v_parallel .gt. 0 if it is in
! the "co" direction (v_phi*j_phi.gt.0). changes are made in gintrp.for, orbout.for, orbpre.for
  !
  ! this subroutine is similar to bintrp and jintrp in that it loads
  ! data into the orbital ion variables.  this functionality was pulled
  ! out of orbg.  the data loaded for icall=1 includes the function
  ! and first and second derivatives of
  !   g   = r*bt (gauss-cm)
  !   psi = poloidal flux with 2*pi normalization (gauss-cm**2/radian)
  !   phi = electric potential (volts)
  !
  ! when icall=2 only the data required by orbg after the time step is loaded
  !
  !	-----------------------
  !	1-d spline coefficients
  !	-----------------------
  !
  !     real    spliot(4,nximax) !iota(bar) spline
  !	real	splg(4,nximax)	 !ditto - g-mhd  pppl-2089, feb.1984
  !	real	splphi(4,nximax) !ditto - electric potential (phi)
  !
  !	-----------------------
  !	2-d spline coefficients
  !	-----------------------
  !
  !	real	bmodcs(4,momdim,nximax)	  !cos-spline of |b|
  !	real	xjaccs(4,4*momdim,nximax) !cos-spline of jacobian
  !
  !	real	bmodss(4,momdim,nximax)	  !sin-spline of |b|
  !	real	xjacss(4,4*momdim,nximax) !sin-spline of jacobian
  !
  use iso_c_binding, only: fp => c_double 
  use nbi_com
  !
 
  implicit none
  integer :: icall    ! = 1 -> before time step
                      ! = 2 -> after  time step
 
  integer, parameter :: ion = 1    ! only one ion
 
  integer :: jx       ! radial zone index
 
  real(fp)    :: zx0      ! xi at radial zone for spline interpolation
  real(fp)    :: zh       ! offset into radial zone
  real(fp)    :: zkonst   ! toroidal flux / pi with normalization
  real(fp)    :: ziotb    ! iota bar
  real(fp)    :: zdibdxi  ! d(iota_bar)/d(xi)
  real(fp)    :: zf1d,zpsint,zdf1d,zd2f1d
  real(fp), parameter :: zconv_g = 100._fp*1.e4_fp  ! convert m*tesla to cm*gauss
  real(fp), parameter :: zconv_psi = 1.e8_fp  ! convert webers to maxwells
  !
  !       -----
  !       begin
  !       -----
  !
  ! ------------------------------ spline -----------------------------------
  !
  !       --------------------------------------------------
  !       1.0     compute the particles' radial zone numbers
  !       --------------------------------------------------
  !
  !   first figure out which radial zone we're in, i.e. which radial
  !   subscript to use.  here i assume xi(i) = (i-1)*dxi.
  !   for spline interpolation...
  !
  jx  = int(bixi(ion)/dxi(lcentr,2)) + 1
  zx0 = (jx-1)*dxi(lcentr,2)
  zh  = bixi(ion) - zx0
 
  zkonst = zconv_psi*tflux/xpi
 
  if (icall == 1) then
        !
        !       -----------------------------------------
        !       2.0     interpolate the magnetics splines
        !       -----------------------------------------
        !
        !   the spline functions are evaluated as follows
        !
        !	f(x) = ((f(4,j)*h + f(3,j))*h + f(2,j))*h + f(1,j)
        !
        !	df(x) = (3*f(4,j)*h + 2*f(3,j))*h + f(2,j)
        !
        !	d2f(x) = 6*f(4,j)*h + 2*f(3,j)
        !
        !	d3f(x) = 6*f(4,j)
        !
        !   where f(1-4,j) are the spline coefficients, and where
        !   0 .le. h=x-(j-1)*dx .lt. dx.
        !
        !       -----------------------------------------
        !	2.1   1-d splines and their derivatives
        !       -----------------------------------------
        !
     ziotb    = zf1d(zh,spliot(1,jx),spliot(2,jx),spliot(3,jx),spliot(4,jx))

     ziotb    =ziotb*(-nsnccw) !mg, 12/19/2013
     biq(ion) = 1.0_fp/ziotb
     big(ion) = zf1d(zh,splg(1,jx),splg(2,jx),splg(3,jx),splg(4,jx))
!mg, 12/19/2013     bipsi(ion) = zpsint(plflxi(jx)*,zkonst,zx0,zh,spliot(1,jx),  &
!mg, 12/19/2013          spliot(2,jx),spliot(3,jx),spliot(4,jx)) 
     bipsi(ion) = zpsint(plflxi(jx)*(-nsnccw),zkonst,zx0,zh,spliot(1,jx)*(-nsnccw),  &
          spliot(2,jx)*(-nsnccw),spliot(3,jx)*(-nsnccw),spliot(4,jx)*(-nsnccw)) !mg, 12        !
        ! please note that the correct sign according to mike zarnstorff
        ! is co injection results in positive central potential and counter
        ! injection results in negative central potential with respect to the
        ! plasma edge.
        !
 
     if(nlvphi) then
        biphi(ion)    = zf1d(zh,splphi(1,jx),splphi(2,jx),splphi(3,jx),splphi(4,jx))
        bi1phi1x(ion) = zdf1d(zh,splphi(2,jx),splphi(3,jx),splphi(4,jx))  !d(phi)/d(xi)
        bi2phi2x(ion) = zd2f1d(zh,splphi(3,jx),splphi(4,jx))              !d2(phi)/d2(xi)
     else
        biphi(ion)    = 0.0_fp
        bi1phi1x(ion) = 0.0_fp
        bi2phi2x(ion) = 0.0_fp
     endif
 
     zdibdxi = zdf1d(zh,spliot(2,jx),spliot(3,jx),spliot(4,jx)) ! d(iotab)/d(xi)
     zdibdxi =zdibdxi*(-nsnccw) !mg, 12/19/2013
     bidqdx(ion) = -biq(ion)*biq(ion)*zdibdxi                   ! d(q)/d(xi)
 
     bi1g1x(ion) = zdf1d(zh,splg(2,jx),splg(3,jx),splg(4,jx))  !d(g)/d(xi)
     bi2g2x(ion) = zd2f1d(zh,splg(3,jx),splg(4,jx))            !d2(g)/d2(xi)
 
     bi1psi1x(ion) = zkonst*ziotb*bixi(ion)              ! d(psi)/d(xi)
     bi2psi2x(ion) = zkonst*(ziotb + bixi(ion)*zdibdxi)  ! d2(psi)/d2(xi)
 
  else
        ! icall=2
 
     big(ion)   = zf1d(zh,splg(1,jx),splg(2,jx),splg(3,jx),splg(4,jx))
!mg, 12/19/2013     bipsi(ion) = zpsint(plflxi(jx),zkonst,zx0,zh,spliot(1,jx), &
!mg, 12/19/2013          spliot(2,jx),spliot(3,jx),spliot(4,jx))
     bipsi(ion) = zpsint(plflxi(jx)*(-nsnccw),zkonst,zx0,zh,spliot(1,jx)*(-nsnccw), &
          spliot(2,jx)*(-nsnccw),spliot(3,jx)*(-nsnccw),spliot(4,jx)*(-nsnccw))
     if(nlvphi) then
        biphi(ion) = zf1d(zh,splphi(1,jx),splphi(2,jx),splphi(3,jx),splphi(4,jx))
     else
        biphi(ion) = 0.0_fp
     endif
 
  end if   ! end of icall test
 
end subroutine gintrp

subroutine psi_spline(zx,jx,zpsi)

  !  using the same method as gintrp, return psi(x) & also the x zone index
  !  but: psi(x) returned in wb = t*m**2!
  !    intent: replace fnbint calls on plflxg with this...

  use iso_c_binding, only: fp => c_double
  use nbi_com
  implicit none

  real(fp), intent(in) :: zx    ! flux coordinate in (sqrt(phitor/phitor(a)))
  integer, intent(out) :: jx  ! integer zone index out (can be beyond zx=1)
  real(fp), intent(out) :: zpsi ! poloidal flux psi(x), wb

  !------------------------------
  !  local:

  real(fp), parameter :: zconv_psi = 1.e8_fp  ! convert webers to maxwells
  real(fp), parameter :: zunconv = 1.e-8_fp
  real(fp) :: zxi,zkonst,zx0,zh,zpsint
  integer :: jxi
  !------------------------------

  zxi = max(0.0_fp,zx)
  jxi  = int(zxi/dxi(lcentr,2)) + 1
  zx0 = (jxi-1)*dxi(lcentr,2)
  zh  = zxi - zx0
  
  zkonst = zconv_psi*tflux/xpi
 
  zpsi = zpsint(plflxi(jxi),zkonst,zx0,zh,spliot(1,jxi),  &
       spliot(2,jxi),spliot(3,jxi),spliot(4,jxi))

  zpsi = zunconv*zpsi
  jx = lcm1 + jxi

end subroutine psi_spline
