!-----------------------------------------------------------------------
!  dgphilb:  boozer "nu" integrand
!    (related to field phi-lengt integrand: dlp*bt/(r*bp))
!
!      dnu/dtheta = q*d(theta-star)/dtheta - dlp*bt/(r*bp))
!         where d(theta-star)/dtheta is as evaluated in subroutine
!         dgthsta
!
!  (mod dmc 8 sep 1998):  d(theta-star)/dtheta has already been evaluated;
!         just evaluate the dlp*bt/(r*bp) part here, which is pos. definite.
!
!  derivation of transp -> boozer coords transformation integrals is available
!  from dmc as a set of notes...
!
!
subroutine dgphilb(zxiarg,id,zth,zf,ipts,iunerr,ierr)
  use iso_c_binding, only: fp => c_double
  !
  !
  !     this integral over a range of theta gives the transformation
  !     to a change in toroidal angle phi gotten by following a
  !     magnetic field line.
  !
  !     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.
  !
  !
  !       ---------------
  !       input variables
  !       ---------------
  !
  !
  implicit none
  real(fp) ::	zxiarg      ! normalized radial position argument
  integer :: id                        ! zone index
  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
  real(fp) :: zxi,zdl,zrpos,zypos,zb,zbp,zbt,zjac(2,2)
  !
  !  currently unused:  iunerr, id
  !
  !   form the integrand, at each point on the list of theta values
  !
  ierr=0
  !
  zxi=max(1.0e-4_fp,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
     !  here is dnu/dtheta (less dtheta*/dtheta part);
     !  dtheta element will be mult. in by integrator
     zf(i)=zdl*zbt/(zrpos*zbp)
  end do
  !
  return
end subroutine dgphilb
