subroutine nbi_ncift(jzone,jzb,zft,zftminus,zftplus)
  use iso_c_binding, only: fp => c_double
  !
  !  from nbi_com compute ft (trapping fraction) for the upcoming
  !  neoclassical calculation.
  !
  !  called from jbshld
  !
  !  see houlberg's report, eqns 31--33.
  !  use approximate <b>/bmax and <b**2>/bmax**2 for <h> and <h**2>.
  !
  !  input:  transp zone index and zone/bdy shift
  !  output:  trapping fraction zft
  !
  !--------------------------------------------------------
  !
  use nbi_com, only : gb1zc, gb2zc, gncfhzc, bmaxm ,ledge

  implicit none
  !
  !--------------------------------------------------------
  !
  !============
  integer :: jzb,jzone,jzp,jzm
  !============
  real(fp) :: zft,zftminus,zftplus,zbmax,zb1,zb2,zfh,zh,zh2
  !============
  if(jzb.eq.1) then
     !  bdy; avg btw 2 neighboring zones
     jzp=min(ledge,jzone)
     jzm=jzone-1
     zbmax=bmaxm(jzone)  ! bdy ctr'd
     zb1=0.5_fp*(gb1zc(jzp)+gb1zc(jzm))
     zb2=0.5_fp*(gb2zc(jzp)+gb2zc(jzm))
     zfh=0.5_fp*(gncfhzc(jzp)+gncfhzc(jzm))  ! eqn 32 <fcn of h>
  else if(jzb.eq.2) then  ! zone centered
     zbmax=0.5_fp*(bmaxm(jzone)+bmaxm(jzone+1))
     zb1=gb1zc(jzone)
     zb2=gb2zc(jzone)
     zfh=gncfhzc(jzone)
  end if
  !
  zh=min(1.0_fp,zb1/zbmax)
  zh2=min(1.0_fp,zb2/(zbmax*zbmax))
  !
  zftplus  = 1.0_fp- zh2/(zh*zh) * &
       (1.0_fp- sqrt(1.0_fp-zh)*(1.0_fp+0.5_fp*zh))
  zftminus = 1.0_fp- zh2*zfh
  !
  zft = 0.75_fp*zftplus + 0.25_fp*zftminus
  !
  return
end  subroutine nbi_ncift
