!-----------------------------------------------------------------------
!  dgthsta:  flux surface integrand:  dlp*(b**2/bp)/[i + gq]
!    this is dtheta/dtheta, where theta is the straight field line
!    theta used by roscoe white et al (aka boozer coordinates; see
!    gthets.for comments).
!
subroutine dgthsta(zxiarg,id,zth,zf,ipts,iunerr,ierr)
  use iso_c_binding, only: fp => c_double
  !
  !  dmc june 1996 -- based loosely on bob mccann's dvtor routine
  !
  !     it is assumed that momry can be called to get the jacobian for
  !     the volume element and for evaluation of the fields.  this
  !     involves spline interpolation of the moments and evaluation
  !     of 1st derivatives in the xi and theta direction when the
  !     jacobian is extracted.
  !
  !  **note** this routine is intended to integrate along the numeric
  !     flux surfaces, not btw them; the id argument should point to the
  !     zone whose *outer* surface is to be traced, and zxiarg should
  !     give the correct xi of that surface.  example:  if nzones=20,
  !     id=lcentr and zxiarg=0.05 would be correct.
  !
  !------------------------------
  use nbi_com
  !
  !       ---------------
  !       input variables
  !       ---------------
  !
  !
  implicit none
  real(fp) ::	zxiarg      ! normalized radial position argument
  integer :: id                        ! zone index (ignore)
  integer ::	ipts        !  length of list
  real(fp) ::	zth(ipts)   ! list of theta positions
  integer :: iunerr
  !
  !
  !       ----------------
  !       output variables
  !       ----------------
  !
  !
  real(fp) ::	zf(ipts)
  integer :: ierr
  !
  !
  !       ------------------
  !       local declarations
  !       ------------------

  integer :: i,jx,iaxis
  real(fp) :: zh,zx0,zg,zq,zi,zxi,zxismall,zfluxu,zrpos,zypos,zb,zbp,zbt
  real(fp) :: zd2,zlim,zdl,zbf,zjac(2,2)
  real(fp) :: zf1d
  !
  !  currently unused:  id, iunerr
  !
  ierr=0
  !
  !  check arguments:  zxiarg=0.0 or zxi=outer bdy of zone required.
  !
  jx=int(zxiarg/dxi(lcentr,2)) + 1
  zx0=(jx-1)*dxi(lcentr,2)
  zh=zxiarg-zx0
  !
  zg=1.0e-4_fp*zf1d(zh,splg(1,jx),splg(2,jx),splg(3,jx),splg(4,jx))
  zq=zf1d(zh,spliot(1,jx),spliot(2,jx),spliot(3,jx),spliot(4,jx))
  zq=1.0_fp/zq
  zi=1.0e-4_fp*zf1d(zh, &
       splimhd(1,jx),splimhd(2,jx),splimhd(3,jx),splimhd(4,jx))
  !
  if(zxiarg.eq.0.0_fp) then
     iaxis=1
  else
     iaxis=0
  end if
  !
  if(iaxis.eq.1) then
     !
     !  limiting expression, numerically evaluated, for axis
     !
     zxismall=1.0e-4_fp
     !
     zfluxu=1.e4_fp*tflux
     !
     do i = 1,ipts
        if(nfast_map.eq.1) then
           call nbmomry(zxismall,zth(i),1,zrpos,zypos,zjac) ! use bilinear maps
        else
           call nbmomry_xs(zxismall,zth(i),1,zrpos,zypos,zjac) !use xplasma spline maps
        end if
        zd2=(zjac(1,1)*zjac(2,2)-zjac(1,2)*zjac(2,1))
        zlim=(zd2*xpi*zg)/(zrpos*zfluxu*zxismall)
        zf(i)=zlim
     end do
  else
     !
     zxi=zxiarg
     do i = 1,ipts
        call mombs_nb(zxi,zth(i),zrpos,zypos,zjac,zb,zbp,zbt)
        zdl=sqrt(zjac(1,2)**2+zjac(2,2)**2) ! dlp/dtheta
        zbf=zb*zb/zbp               ! b**2/bp
        zf(i)=zdl*zbf/(zi+zg*zq)    ! dtheta mult in by integrator
     end do
  end if
  !
  return
end subroutine dgthsta
! 11sep2004 fgtok
