subroutine magcf(ierr)
  !
  !
  !	last changed:
  !
  !       clf            04/11/08
  !              replace calls to bsub85 entry: gmhde,  qmhde,  imhde, phimhd
  !              with bsub85(...,option= 1, 2, 3, 4)
  !              to solve intel "heap - free" problem
  !
  !  june 1994:  commencing asymmetry upgrade
  !	dmc	** cleanup ** delete computations whose outputs are
  !		never used; corresponding entries in bsub85 also
  !		were deleted:  imhd pmhd omgmhd didtmh dpdtmh
  !			dphimh dodtmh
  !
  !	dmc            11/05/93    nmmjac, nmmbinv made const. vs. xi
  !       tbt            11/04/93    inmomj
  !	ted terpstra    5/30/90  - added ospler
  !	ted terpstra    5/28/90  - added + jerr to error messages.
  !                                - added code to determine the proper
  !                                  number of moments for j and 1/b.
  !                                - put in subroutine nmneed.
  !       bob mccann      6/28/84
  !
  !
  !       --------
  !       comments
  !       --------
  !
  !
  !	this uses m2spl and msplin to do 2-d fourier-spline
  !   transforms and 1-d spline transforms of various quantities
  !   which specify the magnetic field.  these coeffieicnts may
  !   then be used to interpolate these quantities at arbitrary
  !   locations in the plasma.  the primary intended use is to
  !   provide this capability to the beam particle orbit code in
  !   transp.
  !
  !	the functional values used here are supplied by the
  !   entry points in subroutine bsub85.  note that the entry points
  !   must also supply the (d/dxi) values at xi=0.0 and at xi=1.0.
  !
  !	the time derivative functions are supplied by the corresponding
  !    entry points in dbsub85.
  !
  !
  !       -------------------
  !       global declarations
  !       -------------------
  !
  !
  use nbi_com
  use nbi_com, only : ztwopi=>twopi
  !
  !
  !
  !       ---------------
  !       input variables
  !       ---------------
  !
  !
  !	parameters
  !	----------
  !
  !
  !	1-d (momdim)
  !	---
  !
  !	scalars
  !	-------
  !
  !
  !       ----------------
  !       output variables
  !       ----------------
  !
  !		1-d (4,nximax)
  !		---
  !
  !	real	splg		!spline cf's of g-mhd
  !	real	splphi		!spline cf's of phi(xi) (electric potential)
  !
  !  new dmc june 6, 1996 -- investigating replacement of
  !   q, and pol.flx. spline
  !
  !       real    spliot          !spline cf's of iota(bar)
  !       real    splimhd         !spline cf's of i(xi) -- related to
  !                                enclosed current
  !
  !       ------------------
  !       local declarations
  !       ------------------
  !
  !
  implicit none
  integer ::		inxi            !# xi's used
  integer ::		itheta          !# theta's used
  integer ::		iunerr          !unit number for error messages
  !
  real(fp), dimension(:), allocatable :: zxi,zfxi,ztheta
  !
  real(fp) :: zdxi,zximax,zdfxi0,zdfxi1,zdfxi0i,zdfxi1i
  real(fp) :: zkonst,zxi0,zdelxi,zi1,zi2,zi3,zi4,zdelp
  integer :: iximax
  integer :: isplg0,isplg1,ibreak,isplp0,isplp1,isplq0,isplq1
  integer :: i,j,ierr,jerr
  !
  character(len=80) ::	zmsg		!error messages
  !
  !	1-d
  !	---
  !
  !      external	gmhde,  qmhde,  imhde
  !      external	phimhd
  !
  !
  !       ---------------
  !       data statements
  !       ---------------
  !
  !
  !***
  !
  !
  !       -------------------
  !       statement functions
  !       -------------------
  !
  !
  !***
  !
  !
  !=========================================
  !
  !       ----------------------
  !       0.1     initialization
  !       ----------------------
  !
  ierr=0
  !
  allocate(zxi(mtbl),zfxi(mtbl),ztheta(mtheta))
  !
  !
  zdxi = dxizb(lcentr)		!use transp dxi for zone boundaries
  !
  !    compute the maximum xi used anywhere (1/b), etc. are needed
  !    **see dick wieland's code - bloat, trkrb, etc.
  !
  iximax = nznblo + lcentr
  zximax = xiblo(iximax)
  !
  !    now define some ***to become global*** variables
  !
  inxi = 1 + int(zximax/zdxi + 0.5_fp)
  itheta = mtheta
  iunerr = nout		!assume error messages goto this unit #
  !
  if(inxi.le.0 .or. inxi.gt.mtbl) then	
     iximax = mtbl
     write(iunerr,80011)	inxi, iximax
     ierr = 10
     deallocate(zxi,zfxi,ztheta)
     !
     return
  endif
  !
  do i = 1, inxi
     zxi(i) = real(i-1,fp)*zdxi
     
  end do
  !
  do j = 1, itheta
     ztheta(j) = ztwopi*real(j-1,fp)/real(itheta,fp)
     
  end do
  !
  !  dmc -- xplasma based spline setup
  !         will eventually replace magcf/bsub85 fourier-splines,
  !         though some 1d splines will contine
  !
  call nbspline(nonlin,zxi,inxi,itheta,ierr)
  if(ierr.ne.0) return
  !
  !       	-------------------------------
  !       	1.3     compute the spline of g
  !       	-------------------------------
  !
  !
  !
  !
  !   first compute the functional values
  !
  !     call gmhde(inxi,	mtbl,	zxi,
  call bsub85(inxi,	mtbl,	zxi, &
       zfxi,	zdfxi0,	zdfxi1, &
       iunerr,	jerr, 1)
  !
  if(jerr.ne.0)	then
     ierr = 13000+jerr
     zmsg = 'gmhde - compute g_mhd'
     write(iunerr,80999) ierr, zmsg
     deallocate(zxi,zfxi,ztheta)
     !
     return
  endif
  
  !
  !   store the functional values
  !
  do j = 1, inxi
     splg(1,j) = zfxi(j)
  end do
  !
  isplg0 = 1	!first derivative b.c.
  isplg1 = 1	!ditto
  ibreak = 1 + int(0.5_fp+ 1.0_fp/dxizb(lcentr))	!index of xi=1.0
  !
  !   now compute the spline transform
  !
  call msplin(ibreak,mtbl, inxi,zdxi,	splg, &
       isplg0,	zdfxi0, &
       isplg1,	zdfxi1, &
       iunerr,	jerr)
  !
  if(jerr.ne.0)	then
     ierr = 13000+jerr
     zmsg = 'gmhde - compute g_mhd'
     write(iunerr,80999) ierr, zmsg
     deallocate(zxi,zfxi,ztheta)
     !
     return
     
  endif
  !
  !
  !       	-------------------------------
  !       	1.3a     compute the spline of i
  !       	-------------------------------
  !
  !
  !
  !
  !   first compute the functional values
  !
  !      call imhde(inxi,	mtbl,	zxi,
  call bsub85(inxi,	mtbl,	zxi, &
       zfxi,	zdfxi0,	zdfxi1, &
       iunerr,	jerr, 2)
  !
  if(jerr.ne.0) then
     ierr = 13200+jerr
     zmsg = 'imhde - compute i_mhd'
     write(iunerr,80999) ierr, zmsg
     deallocate(zxi,zfxi,ztheta)
     !
     return
  endif
  !
  !   store the functional values
  !
  do j = 1, inxi
     splimhd(1,j) = zfxi(j)
  end do
  !
  isplg0 = 1	!first derivative b.c.
  isplg1 = 1	!ditto
  ibreak = 1 + int(0.5_fp+ 1.0_fp/dxizb(lcentr))	!index of xi=1.0
  !
  !   now compute the spline transform
  !
  call msplin(ibreak,mtbl, inxi,zdxi,	splimhd, &
       isplg0,	zdfxi0, &
       isplg1,	zdfxi1, &
       iunerr,	jerr)
  !
  if(jerr.ne.0) then
     ierr = 13300
     zmsg = 'msplin(i_mhd)       '
     write(iunerr,80999) ierr, zmsg
     deallocate(zxi,zfxi,ztheta)
     !
     return
     
  endif
  !
  !
  !       	-------------------------------
  !       	1.4     compute the spline of q
  !       	-------------------------------
  !
  !
  !
  !
  !   first compute the functional values
  !
  !      call qmhde(inxi,	mtbl,	zxi,
  call bsub85(inxi,	mtbl,	zxi, &
       zfxi,	zdfxi0,	zdfxi1, &
       iunerr,	jerr, 3)
  !
  if(jerr.ne.0) then
     ierr = 14000+jerr
     zmsg = 'qmhde - compute q_mhd'
     write(iunerr,80999) ierr, zmsg
     deallocate(zxi,zfxi,ztheta)
     !
     return
  endif
  !
  !   store the functional values
  !
  do j = 1, inxi
     spliot(1,j) = 1.0_fp/zfxi(j)
  end do
  zdfxi0i=zdfxi0
  zdfxi1i=-zdfxi1/(zfxi(inxi)*zfxi(inxi))
  !
  isplq0 = 1	!first derivative b.c.
  isplq1 = 1	!ditto
  ibreak = 1 + int(0.5_fp+ 1.0_fp/dxizb(lcentr))	!index of xi=1.0
  !
  !   compute spline transform:  iota(bar)
  !
  call msplin(ibreak, mtbl, inxi,zdxi,	spliot, &
       isplq0,	zdfxi0i, &
       isplq1,	zdfxi1i, &
       iunerr,	jerr)
  if(jerr.ne.0)	then
     ierr = 14200+jerr
     zmsg = 'msplin(iota(bar))   '
     write(iunerr,80999) ierr, zmsg
     deallocate(zxi,zfxi,ztheta)
     return
  endif
  !
  !  dmc -- integrate poloidal flux profile using
  !    d(psi)/d(xi) = iota(bar)*xi*tflux/pi
  !
  !    this is an analytic integration of xi* a spline in xi,
  !    i.e. a quintic expression
  !
  !     h=xi-xi(j); xi = h + xi(j)
  !     delxi=xi(j+1)-xi(j)
  !       iota(bar)=i1(j)+i2(j)*h+i3(j)*h**2+i4(j)*h**3
  !
  !     int[0 to delxi]{dh * [xi(j)*i1(j) + h*(i1(j)+xi(j)*i2(j))
  !                     +h**2*(i2(j)+xi(j)*i3(j)) + h**3*(i3(j)+xi(j)*i4(j))
  !                     +h**4*i4(j)] }
  !
  plflxi(1)=0.0_fp
  zkonst=1.e8_fp*tflux/xpi
  do j = 1,inxi-1
     zxi0=zxi(j)
     zdelxi=zxi(j+1)-zxi0
     zi1=spliot(1,j)
     zi2=spliot(2,j)
     zi3=spliot(3,j)
     zi4=spliot(4,j)
     !  quintic expression for integral:
     zdelp=zdelxi*(zxi0*zi1 + zdelxi*((zi1+zxi0*zi2)/2.0_fp &
          + zdelxi*((zi2+zxi0*zi3)/3.0_fp &
          + zdelxi*((zi3+zxi0*zi4)/4.0_fp &
          + zdelxi*zi4/5.0_fp))))
     plflxi(j+1)=plflxi(j)+zkonst*zdelp
  end do
  !
  !
  !
  !   first compute the functional values
  !
  !      call phimhd(	inxi,	mtbl,zxi,
  call bsub85(	inxi,	mtbl,zxi, &
       zfxi,	zdfxi0,	zdfxi1, &
       iunerr,	jerr, 4)
  !
  if(jerr.ne.0)	then
     ierr = 55000+jerr
     zmsg = 'phimhd - compute phi(xi)'
     write(iunerr,80999) ierr, zmsg
     deallocate(zxi,zfxi,ztheta)
     return
  endif
!
!   store the functional values
!
  do j = 1, inxi
     splphi(1,j) = zfxi(j)
  end do
!
  isplp0 = 1	!first derivative b.c.
  isplp1 = 1	!ditto
  ibreak = 1 + int(0.5_fp+ 1.0_fp/dxizb(lcentr))	!index of xi=1.0
  !
  !   now compute the spline transform
  !
  call msplin(ibreak,mtbl, inxi,zdxi,	splphi, &
       isplp0,	zdfxi0, &
       isplp1,	zdfxi1, &
       iunerr,	jerr)
  !
  if(jerr.ne.0) then
     ierr = 55100+jerr
     zmsg = 'msplin(phi(xi))       '
     write(iunerr,80999) ierr, zmsg
     deallocate(zxi,zfxi,ztheta)
     !
     return
  endif

  !
  !
  !
  deallocate(zxi,zfxi,ztheta)
  !
  return
80011 format(' magcf error: inxi=',i5,' out of bounds.',/, &
           ' iximax=',i5,/, &
           ' check trans.com parameters.',/, &
           ' recompile and link and lib/create magcf.')

!
80999 format(' magcf error: ierr=',i6,' on return from subroutine',/, &
           ' ',a)
end subroutine magcf
!
!=========================================
!
! 11sep2004 fgtok
