!  gthets -- generate straight line coordinate system theta spline
!            generate field line phi length spline
!
subroutine gthets(ierr)
  !
  use nbi_com
  !
  !
  !============
  implicit none
  real(fp) :: stat
  !============
  integer, intent(out) :: ierr   ! status code (0=ok)
  !
  !  dmc nov 2001 -- moved into nubeam; converting to xplasma integrator
  !
  !            dmc 13 june 1996
  !
  !  based on geometric integration formula, generate bicubic splines
  !  representing the straight field line (boozer) theta as a function
  !  of transp theta and normalized radial coordinate (xi).
  !
  !  straight field line coordinates (boozer coordinates) are described
  !  in the published literature in various places...
  !
  !    r. b. white, phys. fluids b, vol. 2 no. 4 p. 845 (1990).
  !    r. b. white & m. s. chance, phys. fluids vol. 27, p. 2455 (1984).
  !
  !  in straight line coords the b field takes the forms:
  !
  !    b = grad(psi.tor)^grad(theta*) + grad(psi.pol)^grad(zeta)
  !         and
  !    b = g(x)*grad(zeta) + i(x)*grad(theta*) +
  !            correction to zero out radial component.
  !            (here x is a flux surface label)
  !
  !  and the jacobian j**-1 = grad(psi.pol).(grad(theta*)^grad(zeta))
  !                   j = (i + g*q)/b**2 = f(psi.pol)/b**2.
  !
  !  where psi.pol is the poloidal flux enclosed within each flux surface,
  !  and theta* is chosen to satisfy these relations and zeta is
  !  chosen to make the field lines straight, i.e.
  !    b.grad(zeta)/b.grad(theta*) = q everywhere
  !
  !  and zeta = phi + nu(theta*,x) shows the relationship to the
  !  geometrical toroidal angle.
  !
  !  after some algebra based on the above (notes by d. mccune dated may
  !  1996) integral formulae for theta* and nu in terms of quantities known
  !  to transp ... can be derived:
  !
  !  tflux = enclosed toroidal flux in plasma (scalar)
  !  xi = normalized toroidal flux coordinate
  !       inside given surface, enclused toroidal flux is xi**2*tflux,
  !       this relation defines xi.
  !
  !  and the flux surface functions:
  !
  !  (dv/dxi) = derivative of volume elment
  !  gx2r2i = flux surface average < grad(xi)**2 / r**2 >, already computed.
  !  i = iota(bar) = 1/q
  !  g = r*btoroidal
  !
  !  and the derived flux surface function:
  !  i(xi) = 1/(4*pi**2) * tflux * (dv/dxi) * i * xi * gx2r2i
  !    which is proportional to the enclosed toroidal current
  !
  !  ------------------------------------------------------------------
  !  d(theta*)/d(theta) = (dl(pol.)/d(theta))*[b**2/btheta]/[i + g*q]
  !  ------------------------------------------------------------------
  !
  !  with dl(pol.)/d(theta) = sqrt((dr/dtheta)**2+(dy/dtheta)**2)
  !    from the flux surface representation,
  !
  !  btheta = grad(psi.pol)/r = grad(xi)/r * d(psi.pol)/d(xi)
  !    d(psi.pol)/d(xi) = tflux * i * xi / pi
  !    grad(xi) = sqrt((dr/dtheta)**2+(dy/dtheta)**2)/d
  !      d = ((dr/dxi)*(dy/dtheta) - (dy/dxi)*(dr/dtheta))
  !
  !  bphi = g/r
  !
  !  and b**2 = bphi**2 + btheta**2, ... and ...
  !
  !  ------------------------------------------------------------------
  !  d(nu)/d(theta) = q*d(theta*)/d(theta) - (dl(pol.)/d(theta))*bphi/(r*btheta)
  !  ------------------------------------------------------------------
  !
  !  now the component dl(pol.)*bphi/(r*btheta)
  !   gives the change d(philb) in geometric toroidal angle, along the field
  !   line, for a step dl(pol.) along the surface in the poloidal direction.
  !
  !  thus nu contains the relationship between the poloidal length and the
  !  toroidal length of a field line (this gives its field line straightening
  !  properties).
  !
  !  ** modified ** dmc 20 june -- but:  nu and theta* are stored & splined.
  !  it can be shown nu=0 on axis for all theta, and nu(xi,-pi)=nu(xi,pi)=0;
  !  better for splining than the phi length itself).
  !
  !  it can be shown that theta* integrates up to 2*pi and nu integrates
  !  up to zero, i.e. nu is periodic.
  !
  !----------------
  !
  !  the integration formulae are singular on axis, so, the integration
  !  only done on non-singular surfaces, while the axial variation with
  !  theta is approximated based on the results on the innermost two
  !  non-singular surfaces (this is done for theta*; for nu, nu=0 on
  !  axial b.c. can be shown and is used).
  !
  !----------------
  !
  !  the starting point of the integration on each surfaces is at
  !  theta=-pi; the initial value for theta* at theta=-pi is taken as -pi,
  !
  !  nu is taken as zero at theta = +/-pi
  !
  !----------------
  !----------------
  !  old local stuff
  !
  integer :: ith,iadr,ifree,isize,iwksize,iwork,iwork2,j
  integer :: ix,ier
  integer :: iii,ilinx,ilinth,ia,ia2,istat,idum
  real(fp) :: zxi,zth,zq
  !
  !----------------
  !  new local stuff
  !
  integer :: inchi,inrho,id_rho=0,id_chi=0
  real(fp), dimension(:), allocatable :: zxchi,zxrho  ! xplasma grids
  real(fp), dimension(:,:), allocatable :: zdgthsta,zdgphilb
  !
  integer :: inumchi,inumrho
  real(fp), dimension(:), allocatable :: zichi,zirho  ! integration grids
  real(fp), dimension(:,:), allocatable :: zigthsta,zigphilb
  real(fp) :: zpi,z2pi,zximin
  real(fp), dimension(:), allocatable :: zwk,zwkth
  real(fp) :: zthr8,zthr8p,zphr8,zphr8p,znucr8,znorm8,ztmp
  !
  !----------------
  !  external integrand functions
  !
  external dgthsta                  ! thets (theta-star) integrand
  !
  external dgphilb                  ! field-line length in phi integrand
  !
  !----------------
  !
  integer :: iaddr,iibase,iix,iinx,iith
  !
  !  statement fcn
  !
  iaddr(iibase,iix,iinx,iith) = iibase + 16*((iix-1)+iinx*(iith-1))
  !
  !----------------
  !
  !  start of executable code
  !
  ierr = 0
  !
  zpi=abs(acos(-1.0_fp))
  z2pi=zpi+zpi
  !
  !  0.  xplasma grids
  !
  call eq_ganum('__rho',id_rho)
  call eq_ganum('__chi',id_chi)
  call eq_ngrid(id_rho,inrho)
  call eq_ngrid(id_chi,inchi)
  allocate(zxrho(inrho),zxchi(inchi))
  call eq_grid(id_rho,zxrho,inrho,idum,ierr)
  call eq_grid(id_chi,zxchi,inchi,idum,ierr)
  !
  allocate(zdgthsta(inchi,inrho))   ! for integrand data
  allocate(zdgphilb(inchi,inrho))   ! for integrand data
  !
  !  1.  the theta grid and rho integration grids:
  !
  inumchi=nthspl
  allocate(zichi(inumchi))
  inumrho=nzp1
  allocate(zirho(inumrho))
  zirho=xi(lcentr:lep1,1)
  !
  allocate(zigthsta(inumchi,inumrho))   ! for integration result
  allocate(zigphilb(inumchi,inumrho))   ! for integration result
  !
  lsptheta=1
  !
  ifree=lsptheta+nthspl
  !
  !  size of spline structures
  !
  isize=16*nthspl*nzp1
  !
  !  allocate space
  !
  lspthets=ifree
  ifree=ifree+isize
  !
  lspphilb=ifree
  ifree=ifree+isize
  !
  !  workspace
  !
  iwksize=max(nzp1*(nthspl+1),5*max(nzp1,nthspl))
  iwork=ifree
  iwork2=ifree+nzp1*nthspl
  ifree=ifree+iwksize
  !
  !  ok now make the datspl array...
  !
  midspl=2*(ifree-1)
  if(allocated(datspl)) deallocate(datspl)
  allocate(datspl(midspl))
  datspl=0.0_fp
  !
  !  simple cobble to avoid axial singularity
  !
  zirho(1)=zirho(2)/4
  zximin=0.5_fp*zirho(1)
  !
  call eq_flxint_init(0,inumrho,zirho,zximin,ierr)
  if(ierr.ne.0) then
     write(nonlin,*) ' ?gthets: eq_flxint_init error.'
     goto 1000
  end if
  !
  !  theta grid
  !
  do ith=1,nthspl
     zthr8 = -zpi + real(ith-1)*(z2pi)/real(nthspl-1)
     zichi(ith)=zthr8
     zth=zthr8
     iadr=lsptheta+ith-1
     datspl(iadr)=zth
  end do
  !
  call eq_flxint_chinit(0,inumchi,zichi,ierr)
  if(ierr.ne.0) then
     write(nonlin,*) ' ?gthets: eq_flxint_cinit error.'
     goto 1000
  end if
  !
  !  prepare xplasma integrand
  !
  allocate(zwk(inchi),zwkth(inchi))
  zwkth=zxchi
  do j=1,inrho
     zxi=zxrho(j)
     call dgthsta(zxi,0,zwkth,zwk,inchi,nonlin,ierr)
     zdgthsta(1:inchi,j)=zwk
     !
     call dgphilb(zxi,0,zwkth,zwk,inchi,nonlin,ierr)
     zdgphilb(1:inchi,j)=zwk
  end do
  deallocate(zwk,zwkth)
  !
  !  integrate
  !
  call eq_flxint_arr2(zdgthsta,1,zigthsta,inumchi,inumrho,ierr)
  if(ierr.ne.0) then
     write(nonlin,*) ' ?gthets: eq_flxint_arr2 error (dgthsta)'
     goto 1000
  end if
  call eq_flxint_arr2(zdgphilb,1,zigphilb,inumchi,inumrho,ierr)
  if(ierr.ne.0) then
     write(nonlin,*) ' ?gthets: eq_flxint_arr2 error (dgphilb)'
     goto 1000
  end if
  !
  !  convert:  remove twopi factor from dphi integration; turn integral
  !     pieces into running sum, convert nu.
  !
  do j=1,inumrho
     zq=1.0_fp/spliot(1,j)
     zthr8=0
     zphr8=0
     do ith=1,inumchi
        zthr8p=zthr8
        zphr8p=zphr8
        if(ith.lt.inumchi) then
           zthr8=zthr8p+zigthsta(ith,j)/(z2pi)
           zphr8=zphr8p+(zq*zigthsta(ith,j)-zigphilb(ith,j))/(z2pi)
        end if
        zigthsta(ith,j)=zthr8p
        zigphilb(ith,j)=zphr8p
     end do
  end do
  !
  !
  !  renormalization & phase shift
  !
  !  thets should integrate to 2*pi; apply shift of -pi to thets.
  !  nu should integrate to zero.
  !
  do j=lcentr,lep1
     ix=j-lcentr+1
     znucr8=zigphilb(nthspl,ix)
     znorm8=(z2pi)/zigthsta(nthspl,ix)
     do ith=1,nthspl
        zigthsta(ith,ix)=znorm8*zigthsta(ith,ix)-zpi
        if(j.eq.lcentr) then
           zigphilb(ith,ix)=0
        else
           !  linear shift to nu to get exact periodicity
           zigphilb(ith,ix)=zigphilb(ith,ix)- &
                (ith-1)*znucr8/(nthspl-1)
        end if
     end do
  end do
  !
  !  copy data into spline buffer...
  !
  do ith=1,nthspl
     do ix=1,nzp1
        ia=iaddr(lspthets,ix,nzp1,ith)
        datspl(ia)=zigthsta(ith,ix)
        ia=iaddr(lspphilb,ix,nzp1,ith)
        datspl(ia)=zigphilb(ith,ix)
     end do
  end do
  !
  !  debug plots
  !
  !dbg      do ith=1,nthspl
  !dbg         do ix=1,nzp1
  !dbg            ia=iaddr(lspthets,ix,nzp1,ith)
  !dbg            ia2=iwork+(ix-1)+nzp1*(ith-1)
  !dbg            datspl(ia2)=datspl(ia)
  !dbg         end do
  !dbg      end do
  !
  !dbg      call grf3f1(xizb(lcentr),datspl(iwork),datspl(lsptheta),
  !dbg     >            datspl(iwork2),nzp1,nthspl,nzp1,
  !dbg     >            'xi','theta(boozer)','theta(transp)',
  !dbg     >            'xi','radians','radians',
  !dbg     >            'debug plot (gthets.for)',runid,0)
  !
  !dbg      do ith=1,nthspl
  !dbg         do ix=1,nzp1
  !dbg            ia=iaddr(lspphilb,ix,nzp1,ith)
  !dbg            ia2=iwork+(ix-1)+nzp1*(ith-1)
  !dbg            datspl(ia2)=datspl(ia)
  !dbg         end do
  !dbg      end do
  !
  !dbg      call grf3f1(xizb(lcentr),datspl(iwork),datspl(lsptheta),
  !dbg     >            datspl(iwork2),nzp1,nthspl,nzp1,
  !dbg     >            'xi','nu(s.l. phi corr.)','theta(transp)',
  !dbg     >            'xi','radians','radians',
  !dbg     >            'debug plot (gthets.for)',runid,0)
  !
  !  now compute the spline fits
  !
  call bpspline(xizb(lcentr),nzp1,datspl(lsptheta),nthspl, &
       datspl(lspthets),nzp1,datspl(iwork),iwksize, &
       ilinx,ilinth,ier)
  if(ier.ne.0) then
     write(nonlin,*) 'gthets:  bpspline error.'
     ierr=1; return
  end if
  !
  if(ilinx.ne.1) then
     write(nonlin,*) &
          'gthets:  bpsline xizb even spacing test failure.'
     ierr=1; return
  end if
  if(ilinth.ne.1) then
     write(nonlin,*) &
          'gthets:  bpsline theta array even spacing test failure.'
     ierr=1; return
  end if
  !
  call bpspline(xizb(lcentr),nzp1,datspl(lsptheta),nthspl, &
       datspl(lspphilb),nzp1,datspl(iwork),iwksize, &
       ilinx,ilinth,ier)
  if(ier.ne.0) then
     write(nonlin,*) 'gthets:  bpspline error.'
     ierr=1; return
  end if
  !
  if(ilinx.ne.1) then
     write(nonlin,*) &
          'gthets:  bpsline xizb even spacing test failure.'
     ierr=1; return
  end if
  if(ilinth.ne.1) then
     write(nonlin,*) &
          'gthets:  bpsline theta array even spacing test failure.'
     ierr=1; return
  end if
  !
  !  normal exit
  !
  lspuse=iwork-1    ! reserve datspl(1:lspuse)
  !
  !mg       nphipzc=1
  !
  !dbg      call chkspl
  !
  goto 1000
  !
  !  error
  !
999 continue
  !
  iii=midspl
  write(nonlin,9999) iii

9999 format(' ?gthets -- datspl (spline) buffer size exceeded'/ &
       '  currently idspl = ',i7/ &
       '  (this error is unexpected -- likely bug in startup code.')
  ierr=1
  !
1000 continue
  deallocate(zirho,zichi,stat=istat)
  deallocate(zxrho,zxchi,stat=istat)
  deallocate(zdgthsta,zdgphilb,stat=istat)
  deallocate(zigthsta,zigphilb,stat=istat)
  return
  !
end subroutine gthets
! 11sep2004 fgtok
