subroutine bsub85(ixi,	iximax, &
     zxi, &
     zfxi,	zdfxi0,	zdfxi1, &
     lunerr,	ierr, option)
  !
  !       clf 11 apr 2008 -- replace entry: gmhde,  qmhde,  imhde, phimhd
  !                          with "case option" (1,2,3,4)
  !                          intel had problems with dynamic arrays
  !       dmc 29 dec 2001 -- much of this code has been replaced with
  !       bicubic spline interpolation software.  only 1d spline fits
  !       are still done here.
  !
  !       bob mccann, may 1985
  !
  !		    only the entries are actually used!
  !		    the order they are called is set in magcf
  !
  !	all quantities are computed for t=tgeo in order to use the
  !	most current metric which is known.
  !	in transp there is a relevant naming convention: 2-d quantities
  !	which have the ",iz)" ending in port.for will have two 1-d arrays
  !	equivalenced to it---one will end in "zc" for zone centers and
  !	the other will end in "zb" for zone boundary.  some of the 1-d
  !	variables end in "g1", "g2", and "geo" corresponding to the times
  !	tg1, tg2, and tgeo which are the beginning, end, and current times
  !	within the current magnetics time interval.
  !
  !
  !       -------------------
  !       global declarations
  !       -------------------
  !
  !
  use nbi_com
  !
  !
  !       ---------------
  !       input variables
  !       ---------------
  !
  !
  !============
  implicit none
  integer :: lunerr,ithmaxp,ipbdy,i,ibdy,ierr
  !============
  real(fp) :: zsvtov,zlope0,delta,f0,f0pj,f0p2j,f0p3j,zlope1,f1,f1mj
  real(fp) :: f1m2j,f1m3j,zfe,zfem1,zfem2,zfem3,zdxtnd,zf0,zf1
  real(fp) :: zf2,zf3,zkonst,zqx1,zqx2,zxmix,zxmix3,zqmix,zmix
  real(fp) :: ztrght2
  !============
  integer ::	iximax		!max # xi's available
  !
  real(fp) ::	zxi(iximax)		!ditto - xi
  !
  integer ::	ixi		!# xi's to use
  !
  !
  !       ----------------
  !       output variables
  !       ----------------
  !
  !
  real(fp) ::	zfxi(iximax)		!dependent variable f(xi)
  real(fp) ::	zdfxi0			!df(xi)/d(xi) at xi=0
  real(fp) ::	zdfxi1			!df(xi)/d(xi) at xi=zdx*(ixi-1)
  !
  !
  integer ::   option
  !       ------------------
  !       local declarations
  !       ------------------
  !
  !	parameter	ithmaxp = 32	!.ge.2*(mimmsta-1) *** see magcf.for
  parameter	(ithmaxp = 64)  	!.tbt      5/8/90
  !
  real(fp) ::	zdx		!dxizb(lcentr)
  !
  !-----------------------------------
  !  dmc 29 dec 2001 -- bicubic splines replace all this...
  !
  !
  !-----------------------------------
  !
  !    the following are strictly local to one entry point:
  !
  real(fp) ::	zigeo(mtbl)	!toroidal current within flux surface at t=tgeo
  !
  !
  !    debug stuf
  !
  !       ---------------
  !       data statements
  !       ---------------
  !
  !
  data	zsvtov/300.0_fp/		!convert statvolts to volts
  !
  !
  !
  !	---------------------------------
  !	2.0	compute the 1-d functions
  !	---------------------------------
  !
  !
  !
  !       	---------------------
  !       	2.1     compute g-mhd
  !			in gauss-cm
  !       	---------------------
  !
  !
  !	according to roscoe, 2*pi*g is the poloidal current outside a given
  !    flux surface.
  !
  !
  select case (option)
     !      entry gmhde
  case (1)
     !
     !
     ierr = 0
     zdx = dxizb(lcentr)			!independent of xi
     ipbdy = 1 + int(0.5_fp+ 1.0_fp/zdx)		!index for xi=1.0
     !
     do i = lcentr, lep1
        ibdy = i - lcentr + 1
        zfxi(ibdy) = bzxr*fbz(i,1)*1.0e4_fp	!convert tesla-cm to gauss-cm

     end do
     !
     !    extrapolate into the bloated region
     !
     zfe =   fbz(lep1,1)*bzxr*1.0e4_fp	!convert tesla-cm to gauss-cm
     zfem1 = fbz(ledge,1)*bzxr*1.0e4_fp	!convert tesla-cm to gauss-cm
     zfem2 = fbz(ledge-1,1)*bzxr*1.0e4_fp	!convert tesla-cm to gauss-cm
     zfem3 = fbz(ledge-2,1)*bzxr*1.0e4_fp	!convert tesla-cm to gauss-cm
     !
     do i = lep1+1, ixi-1+lcentr
        ibdy = i - lcentr + 1
        zdxtnd = real(i-lep1)*zdx
        !
        !	  zfxi(ibdy) = ztrght2(zdxtnd,zfe,zdx,
        !     1				zfem1,zfem2,zfem3)
        !
        !    use the fact that therer is no poloidal current between the plasma
        !    boundary and the tf coils
        !
        zfxi(ibdy) = zfxi(lep1-lcentr+1)

     end do
     !
     zf0 = zfxi(1)
     zf1 = zfxi(2)
     zf2 = zfxi(3)
     zf3 = zfxi(4)
     !
     !	zdfxi0 = zlope0(zdx,zf0,zf1,zf2,zf3)
     !
     !    the current inside a flux surface is proportional to area
     !    so the slope goes to zero
     !
     zdfxi0 = 0.0_fp		!physically correct
     !
     !	zdfxi1 = zlope1(zdx,zfe,zfem1,zfem2,zfem3)
     !
     !    use the fact that there is no current between the plasma edge
     !    and the tf coils.
     !
     zdfxi1 = 0.0_fp		!physically correct
     !
     !    entry imhde
  case (2)
     !
     !
     !       	---------------------
     !       	2.1     compute i-mhd
     !			in gauss-cm
     !       	---------------------
     !
     !
     !	according to roscoe, 2*pi*g is the poloidal current outside a given
     !    flux surface.
     !
     !
     !
     !
     ierr = 0
     zdx = dxizb(lcentr)			!independent of xi
     ipbdy = 1 + int(0.5_fp+ 1.0_fp/zdx)		!index for xi=1.0
     !
     zfxi(1)=0.0_fp
     !
     zigeo = 0.0_fp
     zigeo(1:lep1)=curt(1:lep1)
     !
     zkonst=1.0e6_fp*xmu/twopi            ! amps -> gauss*cm
     !
     do i = lcentr+1, lep1
        ibdy = i - lcentr + 1
        zfxi(ibdy) = zkonst*zigeo(i)
     end do
     !
     !    extrapolate into the bloated region
     !
     do i = lep1+1, ixi-1+lcentr
        ibdy = i - lcentr + 1
        !
        !    use the fact that therer is no toroidal current between the plasma
        !    boundary and the tf coils
        !
        zfxi(ibdy) = zfxi(lep1-lcentr+1)
     end do
     !
     zf0 = zfxi(1)
     zf1 = zfxi(2)
     zf2 = zfxi(3)
     zf3 = zfxi(4)
     !
     zdfxi0 = zlope0(zdx,zf0,zf1,zf2,zf3)
     !
     !    use the fact that there is no current between the plasma edge
     !    and the tf coils.
     !
     zdfxi1 = 0.0_fp		!physically correct
     !
     !
     !     entry qmhde
  case (3)
     !
     !
     !       	---------------------
     !       	2.2     compute q-mhd
     !			no units
     !       	---------------------
     !
     !
     !
     !      entry qmhde
     !
     !
     ierr = 0
     zdx = dxizb(lcentr)			!independent of xi
     ipbdy = 1 + int(0.5_fp+ 1.0_fp/zdx)		!index for xi=1.0
     !
     do i = lcentr, lep1
        ibdy = i - lcentr + 1
        zfxi(ibdy) = qgeo(i)		!we do not care about up-down
        !					!assymetry
        !	  zfxi(ibdy) = qgeo(i)*(-nsnccw)       !adjust sign to conform
        !					       !to rob's ansatz for orbc
        !					       !get cw or ccw orbits ok!
        !				!q positive => j.b > 0  => nsnccw < 0
        !				!q negative => j.b < 0  => nsnccw > 0

     end do
     !
     !    extrapolate into the bloated region
     !
     do i = lep1+1, ixi-1+lcentr
        ibdy = i - lcentr + 1
        zdxtnd = real(i-lep1)*zdx
        !  cont. slope extrapolation
        zqx1 = zfxi(ipbdy)+real(i-lep1)*(zfxi(ipbdy)-zfxi(ipbdy-1))
        !  q~x**2 extrapolation
        zqx2 = zfxi(ipbdy)*(1.0_fp+zdxtnd)**2
        !  mixing function
        zxmix=min(1.0_fp,zdxtnd/0.2_fp)
        zxmix3=zxmix*zxmix*zxmix
        zqmix=(6.0_fp*zxmix-15.0_fp)*zxmix+10.0_fp
        zmix=zqmix*zxmix3
        !
        zfxi(ibdy) = (1.0_fp-zmix)*zqx1 + zmix*zqx2
        !

     end do
     !
     !	zdfxi0 = zlope0(zdx,zfxi(1),zfxi(2),zfxi(3),zfxi(4))
     zdfxi0 = 0.0_fp  !physically correct rtm april 86
     zdfxi1 = zlope1(zdx,zfxi(ixi),zfxi(ixi-1),zfxi(ixi-2), &
          zfxi(ixi-3))
     !
     !
     !     entry phimhd
  case (4)
     !
     !
     !
     !		--------------------------------------------------
     !		2.7	compute phi(xi): radial electric potential
     !					  due to poloidal rotation
     !			    **also convert volts to statvolts**
     !		--------------------------------------------------
     !
     !
     !
     !
     !	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.
     !
     ierr = 0
     zdx  = dxizb(lcentr)			!independent of xi
     ipbdy = 1 + int(0.5_fp+ 1.0_fp/zdx)		!index for xi=1.0
     !
     do i = lcentr, lep1
        ibdy = i - lcentr + 1
        zfxi(ibdy) = phiprg(i)/zsvtov	!convert volts to statvolts

     end do
     !
     !    extrapolate into bloated region
     !
     do i = ipbdy+1, ixi
        zdxtnd = real(i-ipbdy)*zdx
        zfxi(i) = ztrght2(zdxtnd,zfxi(ipbdy),zdx, &
             zfxi(ipbdy-1),zfxi(ipbdy-2), &
             zfxi(ipbdy-3))

     end do
     !
     zdfxi0 = zlope0(zdx,zfxi(1),zfxi(2),zfxi(3),zfxi(4))
     zdfxi1 = zlope1(zdx,zfxi(ixi),zfxi(ixi-1),zfxi(ixi-2), &
          zfxi(ixi-3))
     !
  end select
  return
  !
end subroutine bsub85
!
!=========================================
!
real(fp)  function ztrght2(dx,f1,delta,f1mj,f1m2j,f1m3j)
  !
  !  replacement for statement function -- dmc june 1990:
  !
  !x	ztrght(dx,f1,delta,f1mj,f1m2j,f1m3j) = f1 + dx*
  !x     1			zlope1(delta,f1,f1mj,f1m2j,f1m3j)
  !
  !  this extrapolation function does not permit the quantity extrapolated
  !  to change sign.
  !
  !
  !  statement function:
  !    quadratic extrapolation of slope at xi=zdx*(ixi-1)
  !
  !============
  use iso_c_binding, only: fp => c_double
  implicit none
  real(fp) :: f1,delta,f1mj,f1m2j,f1m3j,dx,zlope1,zslope,zslopa,zans
  real(fp) :: zl,zalf,zfa
  !--------------------------------------------------------
  !
  !
  !  executable code
  !
  zslope=zlope1(delta,f1,f1mj,f1m2j,f1m3j)
  !
  zslopa=abs(zslope)
  zfa=abs(f1)
  !
  if((zslopa.lt.1.e-7_fp*zfa).or.(zslopa.lt.1.e-30_fp)) then
     !
     !  no slope -- default to flat extrapolation
     !
     zans = f1
     !
  else
     !
     zl = f1 / zslope
     !
     zalf = f1*zl
     !
     if (zl .le. 0.0_fp) then
        !
        !  assymptote towards 0.5*f1
        !
        zans = 0.5_fp*f1 - 0.25_fp*zalf / ( -0.5_fp*zl + dx )
        !
     else
        !
        !  assymptote towards 2*f1
        !
        zans = 2.0_fp*f1 - zalf / ( zl +dx )
        !
     end if
  end if
  !
  ztrght2=zans
  !
  return
end function ztrght2
! 11sep2004 fgtok
