C  vacisld - original
C  vacisld2 - modified to get coil name from mgrid file name
C  vacisld3 - no stellarator symmetry version
C  vacisld4 - vacisld2 modified to interpolate vmec surfaces for iota
C             Also added option to bypass creating surfs file from wout file
C             and use existing surfs file - triggered by presence of 2nd arg
C             on command line.
C             Ntarget also read from input
C  vacisld5 - merge vacisld3 & vacisld4
C           NOTE: vacisld5 only does nonsymmetry calculations
C                 and assumes ntar data given for 360 deg, not 1 period
C  vacisld6 - added visldout ( without .ext ) file to collect all runs
C  vacisld7 - calculate B.grad(phi) explicitly
C  vacisld8 - use vmec field with coil field  ( coil field assumed perturbation field )
C  vacisld9 - CORRECTS Island Size Calculation (by factor of sqrt(2*pi*np))
c***********************************************************************
      module Vext
      character*80 ext, coilsfile
      integer, parameter :: iunit = 7
      end module Vext
c************************************************
      module testfield
      use kind_spec
      real(kind=rprec) :: rc, ac, xiota0,xiota1, bmnx
      integer ::  mx, nx
      end module testfield
c********************************************************
      module coil_data_mod
      use kind_spec
      implicit none
c
c coil filament data
      real(kind=rprec), dimension(:), allocatable :: xf,yf,zf,cf
      integer :: nfil, npts
c
c   working arrays for biot-savart
      real(kind=rprec), dimension(:),allocatable :: dx,dy,dz,vx,vy,vz,fa
c
c coil data
      real(kind=rprec), dimension(:), allocatable :: coilcur
      integer, dimension(:), allocatable :: i1, i2, igroup
      character*16, dimension(:), allocatable ::  group
      integer :: ncoil
c
c group data ( groupcur -> VMEC's EXTCUR )
      real(kind=rprec), dimension(:), allocatable :: groupcur
      integer, dimension(:), allocatable :: ncoilgrp
      integer :: ngroup
c
c header data
      character*32  :: top(3)
c
c Coils file extension
c
      character*16 coilext
c
      end module coil_data_mod
c***********************************************************************
      module Vmeshes
      integer, parameter :: nu = 64
      integer, parameter :: nv = 64*3
      integer, parameter :: md = 32
      integer, parameter :: nd = 32*3
      integer, parameter :: mnd = (md + 1)*(2*nd + 1)
      integer, parameter :: nuv = nu*nv
      integer, parameter :: nuvh = nuv/2 + nu
      end module Vmeshes
c***********************************************************************
      module Vsurfdata
      use Vmeshes
      use kind_spec
      integer ::  icm(mnd),icn(mnd),icm1(mnd),icn1(mnd),imn, ms, ns
      real(kind=rprec), dimension(0:md,-nd:nd) :: cr, cz, dcr, dcz
     +              ,cbu,cbv
      real(kind=rprec), dimension(nuv) :: x, y, z, r, snx, sny ,snz, bs
     +         , dsur, ex, ey, ez, bp, bsbp
      real(kind=rprec), dimension(nuv) :: ru, rv, xu, xv, yu, yv, zu, zv
      real(kind=rprec), dimension(nuv) :: rs, xs, ys, zs,bu,bv
      real(kind=rprec), dimension(100) :: wtar, star, dstar, ds2, bnavg
     +      ,bmns, bmnc, bmn,bnmax
      integer, dimension(100) :: mtar, ntar, idstar
      integer :: nts, ntarget, ntotal, nsurf, noreadwout
      real(kind=rprec) ::  js, s, phip, xiota, xiotap, curpol
      real(kind=rprec), dimension(100) :: xiotaps  ! save xiotap for analytic_island routine
      end module Vsurfdata
c***********************************************************************
      module Vprecal
      use Vmeshes
      use kind_spec
      real(kind=rprec) ::  pi,pi2,alp,alu,alv,alvp
      real(kind=rprec), dimension(nv) ::  coh, sih
      real(kind=rprec), dimension(nu,0:md) ::  comu, simu
      real(kind=rprec), dimension(nv,0:nd) ::  conv, sinv
      end module Vprecal
c***********************************************************************
c***********************************************************************
      program main
c
      use Vext
      use Vsurfdata, ONLY: noreadwout
      use testfield
      implicit none
c
      integer narg, id, iargc, getarg
c
c get file extension from command line
c
      narg=iargc()
      noreadwout=0
      if(narg.lt.2)then
        write(*,*) "Usage : xvacisld ext coilsfile [noreadwout]"
        stop
      elseif(narg.gt.2)then
          noreadwout=1  ! don't read wout file; read modified surfs file
      endif
      id=getarg(1,ext)
      id=getarg(2,coilsfile)
c
c read targeted modes and
c get targeted surfaces from wout file and convert to mag coords
c
      if(narg.eq.2)call vmec2mag ! use existing surfs file if not
c
c  read coils file and set currents based on vmec input file ( extcur values )
c
       call get_coils
c
C      read(*,*) rc, ac, xiota0, xiota1
C      read(*,*) mx, nx, bmnx
c
c calculate island size(s) on targeted surfacaes
c
      call island_size
c
      stop
      end
c  Routine to extract flux surfaces from vmec wout file
c   and convert to straight line coordinates
c  uses vmec routines kind_spec.f vparams.f read_wout_mod.f
c  which must be compiled prior to compiling this routine
c**************************************************************
c
      subroutine vmec2mag
      use kind_spec
      use read_wout_mod, ONLY: nfp, ns, mpol, ntor, mnmax,
     +      xm, xn, rmnx => rmnc, zmnx => zmns, lmnx => lmns,
     +      bsubvmn, iota => iotas, phip, readw_only
     + , bsupumn, bsupvmn
      use Vext
      use Vsurfdata, ONLY: nsurf
      implicit none
c
c      parameter (nu=32,nv=32, nmn=(nu/2)*(nv/2*2-1))
      integer, parameter :: nu=32, nv=32, nmn=(nu/2)*(nv/2*2-1)
      real(kind=rprec), dimension(0:nu/2,-nv/2+1:nv/2) ::  dummy,amus,
     #     rmnc_mag,zmns_mag,drmnc_mag,dzmns_mag,bumnc_mag,bvmnc_mag
      real(kind=rprec), dimension(nu,nv) ::  amu,amu_save,r,z,dr,dz

     + , bu, bv
      real(kind=rprec), allocatable, dimension(:) ::
     +   rmnc, zmns, rmnb, zmnb, lmns, drmnc, dzmns,bumnc,bvmnc
      integer, allocatable, dimension(:) :: marray, narray, mb, nb

      real(kind=rprec), allocatable, dimension(:) :: xt, iotap
      integer, allocatable, dimension(:) :: mt, nt
      character*8 char
      real(kind=rprec) :: pi,s,ds,sumr,sumz,udiff,u,v,um,sum1,sum2,f
     +      ,sumdr, sumdz, sumbu, sumbv, xiotamax
      integer :: ierr,imn,js,mnb,nmode,ntarget,i,isurf,ix,iu,iv,m,n,it,
     +        icount
c
      pi=4.*atan2(1.,1.)
c
c From VMEC output, read Fourier coefficients of boundary shape
c and fourier coefficients of contravariant (super) cpts of B
c
      open(unit=iunit,file="wout."//trim(ext),status='unknown')
c
      call readw_only(iunit, ierr)
      close(iunit)
c
c      write(102,*) "xn=",xn
c
c Allocate memory for local variables
c
      allocate(marray(mnmax), narray(mnmax), mb(mnmax), nb(mnmax),
     +         rmnc(mnmax), zmns(mnmax), rmnb(mnmax), zmnb(mnmax),
     +         lmns(mnmax), drmnc(mnmax), dzmns(mnmax),
     +         bumnc(mnmax),bvmnc(mnmax),
     +         xt(ns),iotap(ns),
     +         mt(ns),nt(ns), stat=ierr)
      if(ierr.ne.0)stop " vmec2mag : allocation error "
      do imn=1,mnmax
         marray(imn)=int(xm(imn))
         narray(imn)=int(xn(imn))/nfp
      enddo
c      write(102,"(2i5)") (marray(imn),narray(imn),imn=1,mnmax)
c
c
c calculate d(iota)/ds
c
      ds=1./(ns-1.)
      iotap(2)=(iota(3)-iota(2))/ds
      iotap(ns)=(iota(ns)-iota(ns-1))/ds
      do js=3,ns-1
      iotap(js)=(iota(js+1)-iota(js-1))/2./ds
      enddo
      xiotamax=maxval(iota(2:ns))
c
c weed out zero terms from boundary
c
      mnb=0
      do imn=1,mnmax
      if(abs(rmnx(imn,ns)).gt.1.e-10 .or.
     +   abs(zmnx(imn,ns)).gt.1.e-10 )then
        mnb=mnb+1
        mb(mnb)=marray(imn)
        nb(mnb)=-narray(imn)  ! REMEMBER TO CHANGE SIGN CONVENTION FOR NESCOIL
        rmnb(mnb)=rmnx(imn,ns)
        zmnb(mnb)=zmnx(imn,ns)
      endif
      enddo
c
c
c
c----------------------------------------------
c input resonant surfaces to extract
c
      open(unit=iunit,file="vacisld.param",status="old")
c      read(iunit,*) nmode  ! , ntarget
c      ntarget = 1 ! force 1 mode per surface
      read(iunit,*) nmode, ntarget
      read(iunit,*) (mt(i),nt(i),i=1,nmode)
      close(iunit)
c
      do i=1,nmode
      xt(i)=float(nt(i))/float(mt(i))
      if(xt(i).gt.xiotamax)xt(i)=xiotamax ! fix for targeting near resonances at edge
      nt(i)=-nt(i) ! convert vmec to internal nescoil notation
c      xt(i)=-nfp*float(nt(i))/float(mt(i))
      enddo
c
      write(102,*) " extracting surfaces for iota = "
      write(102,"(10f6.3)") (xt(i),i=1,nmode)
c
c----------------------------------------------
      open(unit=8, file="surfs."//trim(ext), status="replace")
c      write(8,"(2i5,'  /  number of interior surfaces')") nsurf, nfp
c-------------------------------------------------------
c  start big loop over nmode
c
c      write(8,"(2i5,' / nsurfs, nper')") nmode, nfp
      nsurf=0 ! count number of surfaces found
      do isurf=1,nmode
c
c look for surface with specified iota
c
      ix=0
      do js=2,ns-1
      if( (iota(js).le.xt(isurf) .and. iota(js+1).gt.xt(isurf)) .or.
     +    (iota(js).ge.xt(isurf) .and. iota(js+1).lt.xt(isurf)) .or.
     +    (iota(js+1).eq.xt(isurf).and. js+1.eq.ns) ) then ! big if block
          if(iota(js+1) .ne. iota(js) ) then
            f=( xt(isurf) - iota(js) ) / ( iota(js+1) - iota(js) )
          else
            f = 0.
          endif
        nsurf=nsurf+1
        ix=js
         write(102,*) " Using VMEC surface # ", ix,
     +      " with iota = ", iota(ix)
c
      do imn=1,mnmax
      rmnc(imn)=rmnx(imn,ix)*(1.-f) + rmnx(imn,ix+1)*f
      zmns(imn)=zmnx(imn,ix)*(1.-f) + zmnx(imn,ix+1)*f
      lmns(imn)=lmnx(imn,ix)*(1.-f) + lmnx(imn,ix+1)*f
      drmnc(imn) = ( rmnx(imn,ix+1)-rmnx(imn,ix) )*float(ns)
      dzmns(imn) = ( zmnx(imn,ix+1)-zmnx(imn,ix) )*float(ns)
      bumnc(imn)=bsupumn(imn,ix)*(1.-f) + bsupumn(imn,ix+1)*f
      bvmnc(imn)=bsupvmn(imn,ix)*(1.-f) + bsupvmn(imn,ix+1)*f
      enddo
c
c
c
c Vmec gives (as a Fourier series) the function lambda(u,v) s.t.
c u=um+lambda(u,v), where u,v are vmec poloidal and toroidal 
c angles, and um is the magnetic poloidal angle.
c Find (as a Fourier series) the function amu(um,v) s.t.
c um=u+amu(um,v).
c
      do iv=1,nv
         u=0.
         v=(iv-1)*2.*pi/nv
         do iu=1,nu
            um=(iu-1)*2.*pi/nu
            if(iv.eq.1.and.iu.eq.1) then
               amu(1,1)=0.
               goto 777
            endif
            icount=0
 77         continue
            icount=icount+1
            sum1=0.
            sum2=0.
            do imn=1,mnmax
               m=marray(imn)
               n=narray(imn)
               sum1=sum1+lmns(imn)*sin(m*u-n*v)
               sum2=sum2+m*lmns(imn)*cos(m*u-n*v)
            enddo
            udiff=(u-um+sum1)/(1.+sum2)
            u=u-udiff
            if(abs(udiff/u).gt.1.e-12) then
               if(icount.lt.500) goto 77
            endif
            amu(iu,iv)=u-um
 777        continue
         enddo
      enddo
c
      do iu=1,nu
         do iv=1,nv
            amu_save(iu,iv)=amu(iu,iv)
         enddo
      enddo
      call fft2(nu,nv,amu,dummy,amus)
c
cc      call check(dummy,amus,nu,nv,lmns,marray,narray,mnmax,ns)
c
c Transform rmnc, zmns to magnetic angles
      do iu=1,nu
         do iv=1,nv
            amu(iu,iv)=amu_save(iu,iv)
         enddo
      enddo
       do iv=1,nv
         v=(iv-1)*2.*pi/nv
         do iu=1,nu
            um=(iu-1)*2.*pi/nu
            u=um+amu(iu,iv) ! this is u corresponding to given um
            sumr=0.
            sumdr=0.
            sumbu=0.
            sumbv=0.
            do imn=1,mnmax
               m=marray(imn)
               n=narray(imn)
               sumr=sumr+rmnc(imn)*cos(m*u-n*v)
               sumdr=sumdr+drmnc(imn)*cos(m*u-n*v)
               sumbu=sumbu+bumnc(imn)*cos(m*u-n*v)
               sumbv=sumbv+bvmnc(imn)*cos(m*u-n*v)
            enddo
            r(iu,iv)=sumr
            dr(iu,iv)=sumdr
            bu(iu,iv)=sumbu
            bv(iu,iv)=sumbv
         enddo
      enddo     
      call fft2(nu,nv,r,rmnc_mag,dummy)
      call fft2(nu,nv,dr,drmnc_mag,dummy)
      call fft2(nu,nv,bu,bumnc_mag,dummy)
      call fft2(nu,nv,bv,bvmnc_mag,dummy)
c

c
      do iv=1,nv
         v=(iv-1)*2.*pi/nv
         do iu=1,nu
            um=(iu-1)*2.*pi/nu
            u=um+amu(iu,iv) ! this is u corresponding to given um
            sumz=0.
            sumdz=0.
            do imn=1,mnmax
               m=marray(imn)
               n=narray(imn)
               sumz=sumz+zmns(imn)*sin(m*u-n*v)
               sumdz=sumdz+dzmns(imn)*sin(m*u-n*v)
            enddo
            z(iu,iv)=sumz
            dz(iu,iv)=sumdz
         enddo
      enddo     
      call fft2(nu,nv,z,dummy,zmns_mag)
      call fft2(nu,nv,dz,dummy,dzmns_mag)
c
C      write(6,*) '   m    n      rmnc_vmec      zmns_vmec'
      do imn=1,mnmax
         write(6,617) marray(imn),narray(imn),rmnc(imn),zmns(imn)
      enddo
C      write(6,*) '   m    n       rmnc_mag       zmns_mag'
c
      s=float(ix-1)/float(ns-1)
      write(8,"(2i5,1p4e15.7,
     +'  / nmn, ix, s, phip, iota, iotap -VMEC Surface')")
     +  nmn, ix, s, phip(ix), xt(isurf), (iotap(ix)*(1-f)+iotap(ix+1)*f)
c
      do m=0,nu/2-1
         do n=-nv/2+1,nv/2-1
            write(8,617) m,n*nfp,rmnc_mag(m,n),zmns_mag(m,n),
     +                     drmnc_mag(m,n),dzmns_mag(m,n),
     +                     bumnc_mag(m,n),bvmnc_mag(m,n)
 617        format(2i5,1p6e22.12)
         enddo
      enddo
c
      write(8,"(i5,'  / number of target modes ')") ntarget
      do it=1,ntarget
      write(8,"( 2i5,' 1.e4 0. 0  /  m, n, weight, dstar, idstar' )")
     +         mt(isurf)*it, nt(isurf)*it
      enddo
c
      endif   !     big if block
      enddo
c
      if(ix.eq.0)then
        write(102,*) "No surface matches iota = ", xt(isurf)
c        stop
      endif
c
      enddo  !  big loop over nmode
c
      close(8)
c
c-------------------------------------------------------
c
c deAllocate memory for local variables
c
      deallocate(marray,narray,mb,nb,rmnc,zmns,
     +           lmns,xt,iotap,mt,nt)
      return
      end
cc
      subroutine fft2(m,n,x,xmnc,xmns)
c
c 2-D FFT of real function x(theta,phi) tabulated at m theta points and
c n phi points to give Fourier series
c x(theta,phi)=Sum_{k=0,...mmax,l=nmin,nmax} [xmnc*cos(k*theta+l*phi)
c              +xmns*sin(k*theta+l*phi)] where mmax=m/2,
c nmax=n/2 and nmin=-n/2+1.
c
c n and m must be even.
c
c Input real data x(theta, phi) in array x(i,j) 
c with i=1,2,...m; j=1,2,...n such that
c theta_i=(i-1)*2.*pi/m, phi_j=(j-1)*2.*pi/n
c The y(i,j) array elements should all be set sero on input
c
c On output, the Fourier cosine and sine series are stored
c in xmnc(i,j) and xmns(i,j), where i=0,...,m/2 and
c j=-n/2+1,...,n/2
c 
      use kind_spec
      implicit none
      integer :: i, j, k, l, m, n, ifail
      real(kind=rprec) ::
     +    x(m,n),y(m,n),trigm(2*m),trign(2*n),work(2*m*n)
      real(kind=rprec), dimension(0:m/2,-n/2+1:n/2) ::  xmnc,xmns
      character*1,init
c
      do i=1,m
         do j=1,n
            y(i,j)=0.
         enddo
      enddo
c
      init='i'
      ifail=1
      call c06fuF(m,n,x,y,init,trigm,trign,work,ifail)
c      write(102,*) 'ifail = ',ifail
c     
      i=1
      k=0
      do j=1,n
         if(j.le.n/2) l=-(j-1)
         if(j.gt.n/2) l=n-j+1
         xmnc(k,l)=x(i,j)/sqrt(float(m*n))
         xmns(k,l)=y(i,j)/sqrt(float(m*n))
      enddo     
c
      do i=m/2+1,m
         k=m-i+1
         do j=1,n
            if(j.le.n/2) l=-(j-1)
            if(j.gt.n/2) l=n-j+1
            xmnc(k,l)=2.*x(i,j)/sqrt(float(m*n))
            xmns(k,l)=2.*y(i,j)/sqrt(float(m*n))
         enddo
      enddo
      return
      end

c***************************************************************
c
c************************************************
      subroutine get_coils
c
c Routine to generate coils files
c
      use coil_data_mod
!      use vmec_input
      use Vext
      implicit none
      integer :: istat, i
!c
!c
!c get extcur data from vmec input file
!c
!      open(unit=iunit,file="input."//trim(ext),status="old",
!     +      iostat=istat)
!      if(istat.ne.0) stop "docoils : error opening input.ext"
!c
!      read(iunit,indata,iostat=istat)
!      if(istat.ne.0) stop "docoils : error reading input.ext"
!c
!      close(iunit)
!c
!c read in coil set consistent wit extcur data
!c
!C      open(unit=iunit,file="coils."//trim(ext),status="old",
!C     +      iostat=istat)
!c
!c get coil file name from MGRID_FILE inp vmec input
!c
!      i=index(mgrid_file,"mgrid",.TRUE.)
!      if(i.eq.0)i=index(mgrid_file,"MGRID",.TRUE.)
!      mgrid_file(i:i+4)="coils"
!      open(unit=iunit,file=mgrid_file,status="old",
!     +      iostat=istat)
      open(unit=iunit,file=coilsfile,status="old",iostat=istat)
      if(istat.ne.0) stop "docoils : error opening coils.ext"
      call read_coil_data(iunit)
      close(iunit)
c
c save coil file extension
c
      coilext = coilsfile(7:) ! assumes first six characters are "coils."
c
c
      write(*,*)
      write(*,"('Number of Coil Filament:', i8)") npts
      write(*,"('Number of Coils Read in:', i8)") ncoil
      write(*,"('Number of Groups Found :', i8)") ngroup
      write(*,"('Initial Coil Group Currents:')")
      write(*,606) (groupcur(i),i=1,ngroup)
 606  format(1p5e15.5)
c
c
!      do i=1,ngroup
!C  don't use input file extcur currents
!C          groupcur(i)=EXTCUR(i)*groupcur(i)  ! old makegrid routine
!          groupcur(i)=EXTCUR(i)  ! new makegrid routine
!      enddo
c
      call write_coil_data(iunit)
c
      write(*,"('Final Coil Group Currents:')")
      write(*,606) (groupcur(i),i=1,ngroup)
c
      return
      end
c************************************************
      subroutine read_coil_data(iunit)
      use coil_data_mod
      implicit none
c
      real(kind=rprec) :: x,y,z,c
      integer :: istat, i, iunit
c
      ncoil=0          ! coil counter
      nfil=0           ! filament counter
      npts=0           ! point counter ( = nfil+1 )
      istat=0
c
c
c
c count lines and coils ( c=0. ) for array allocation
c
      read(iunit,"(a32)") (top(i),i=1,3)
c
      do while(istat .eq. 0 )
        read(iunit,*,iostat=istat) x,y,z,c
        if(c.eq.0.)ncoil=ncoil+1
        npts=npts+1
      enddo
C      ngroup=ncoil
      ngroup=max(ncoil,100) ! fix for coil files group# exceeds #coils present
      npts=npts-1
c
      nfil=npts-1
      write(102,*) npts," coil filaments read in "
c
      allocate(xf(npts),yf(npts),zf(npts),cf(npts),
     +  dx(nfil),dy(nfil),dz(nfil),
     +  vx(nfil),vy(nfil),vz(nfil),fa(nfil) ,
     +  coilcur(ncoil),i1(ncoil),i2(ncoil),igroup(ncoil),
     +  group(ncoil), groupcur(ngroup), ncoilgrp(ngroup) , stat=istat)
      if(istat.ne.0)stop "read_coil_data : allocatation error"
c
      groupcur(:)=0.   ! current in group
      igroup(:)=0      ! group id by coil
      ncoilgrp(:)=0    ! number of coils in group
c
      rewind(iunit)
c
      read(iunit,"(a32)") (top(i),i=1,3)
c
      ncoil=0
      npts=0
c
      do while(.TRUE.) ! do always
        ncoil=ncoil+1
        i1(ncoil)=npts+1  ! starting filament of coil ncoil
        i=0
        c=1. ! some non zero number
        do while (c.ne.0.)
          i=i+1
          npts=npts+1
         read(iunit,*,err=98,end=98) xf(npts),yf(npts),zf(npts),cf(npts)
          c=cf(npts)
          if(i.eq.1)coilcur(ncoil)=c
        enddo
c
c           read group info
c
        backspace(iunit)
        read(iunit,*) x,y,z,c,igroup(ncoil),group(ncoil)
        ncoilgrp(igroup(ncoil))=ncoilgrp(igroup(ncoil))+1
c
c         Normalize coil currents to first coil in group
c
        if(ncoilgrp(igroup(ncoil)).eq.1) then
            groupcur(igroup(ncoil))=coilcur(ncoil)
            coilcur(ncoil)=1.
        else
            coilcur(ncoil)=coilcur(ncoil)/groupcur(igroup(ncoil))
        endif
c
c
      i2(ncoil)=npts  ! endding filament of coil ncoil
      enddo
 98   ncoil=ncoil-1  ! should end coil with 0. current ( if not, this throw
      npts=npts-1
      ngroup=maxval(igroup) ! number of coils groups
c
      return
      end
c************************************************
      subroutine write_coil_data(iunit)
      use coil_data_mod
      implicit none
c
      integer :: ic, i, ig, iunit
c
c write out new coil_data file
c
      do ic=1,ncoil
        ig=igroup(ic)
        do i=i1(ic),i2(ic)-1
          cf(i)=coilcur(ic)*groupcur(ig)
c          write(iunit,"(4(1pe13.6,2x))") xf(i),yf(i),zf(i),cf(i)
        enddo
        cf(i2(ic))=0.
CC       write(iunit,"(4(1pe13.6,2x),i5,2x,a16)") xf(i),yf(i),zf(i),cf(i),
c     +             igroup(ic),group(ic)
      enddo
c
      return
      end
c***********************************************************************
c***********************************************************************
c***********************************************************************
      subroutine precal1(np)
c
c                                                             01/01/89
c     excepted from NESCOIL subroutine "precal"
c
c
c ----------------------------------------------------------------------
      use Vmeshes
      use Vprecal
      implicit none
      integer :: m, n, ku, kv, np
c
      pi     = 2.*asin(1.)
      pi2    = 2. * pi
      alp    = pi2 / float(np)
      alu    = pi2 / float(nu)
      alv    = pi2 / float(nv)
      alvp   = pi2 / float(nv*np)
c
      do  m=0,md
      do  ku=1,nu
         comu(ku,m)   = cos(m*alu*(ku-1))
         simu(ku,m)   = sin(m*alu*(ku-1))
      enddo
      enddo
c
      do  n=0,nd
      do  kv=1,nv
         conv(kv,n)   = cos(n*alv*(kv-1))
         sinv(kv,n)   = sin(n*alv*(kv-1))
      enddo
      enddo
c
      do kv=1,nv
         coh(kv)= cos(alvp*(kv-1))
         sih(kv)= sin(alvp*(kv-1))
      enddo
c
      return
      end
c***********************************************************************
      subroutine coilfld4 (xp,yp,zp, bx,by,bz)
c
      use kind_spec
      use Vmeshes
      use coil_data_mod
      implicit none
c
      real(kind=rprec) :: xp,yp,zp, bx,by,bz
       real(kind=rprec) :: ax, ay, az, r1, r2, r12, sdot
      real(kind=rprec) :: xp1, yp1, zp1, xp2, yp2, zp2
      integer :: i
c
      do i=1,nfil
      if(cf(i).ne.0.)then
         dx(i) = (xf(i+1) - xf(i))*cf(i)
         dy(i) = (yf(i+1) - yf(i))*cf(i)
         dz(i) = (zf(i+1) - zf(i))*cf(i)

         vx(i) = yf(i)*dz(i) - zf(i)*dy(i)
         vy(i) = zf(i)*dx(i) - xf(i)*dz(i)
         vz(i) = xf(i)*dy(i) - yf(i)*dx(i)

         xp1 = xp - xf(i)
         yp1 = yp - yf(i)
         zp1 = zp - zf(i)

         r1 = sqrt(xp1*xp1 + yp1*yp1 + zp1*zp1)

         xp2 = xp - xf(i+1)
         yp2 = yp - yf(i+1)
         zp2 = zp - zf(i+1)

         r2  = sqrt(xp2*xp2 + yp2*yp2 + zp2*zp2)

         r12    = r2*r1
         fa(i)  = (r2 + r1)/
     $           (r12*(r12 + xp2*xp1 + yp2*yp1 + zp2*zp1))
      else
          fa(i)=0.
      endif

      enddo

      ax     = dot_product(fa,dx)
      ay     = dot_product(fa,dy)
      az     = dot_product(fa,dz)

      bx     = dot_product(fa,vx) - yp*az+zp*ay
      by     = dot_product(fa,vy) - zp*ax+xp*az
      bz     = dot_product(fa,vz) - xp*ay+yp*ax

      bx=bx*1.e-7
      by=by*1.e-7
      bz=bz*1.e-7
c
      return
      end
c
c************************************************************
      subroutine island_size
c
c routine reads in additional interior surfaces
c and gerates entries in G matrix for each resonance targeted
c
      use Vmeshes
      use Vsurfdata, ONLY: nts,mtar,ntar,ds2,nsurf,bnavg,bs,ntarget,
     + bmns, bmnc, bnmax, noreadwout, dsur, rprec,star, xiotap, xiotaps
      use read_wout_mod, ONLY: np => nfp
      use Vext
      use coil_data_mod, ONLY : coilext
      implicit none
      integer ::  i, j, k, isurf,istat,nmode, ii, kk
      real(kind=rprec) :: ds, drho
c
      open(unit=iunit, file="surfs."//trim(ext), status="old")
c      if(noreadwout.eq.1)then ! need data normally found in wout file
c        read(iunit,*) nmode, np  ! number of additional interior surfaces
c        if(nsurf.eq.0)return
c      endif
c
      nts=0   ! total number of targets for all surfaces
C      call precal1(np)
      call precal1(1)
c
      istat=0
      nsurf=0
      do while(istat.eq.0)
        call readsurf(istat)  ! read surface and targeted resonances ( Modified NESCOIL Routine )
       if(istat.eq.0)then
        call surface1   ! evaluate mesh points and normals on surface (  NESCOIL Routine )
        call bdotn      !  calculate b.n on interior surfaces
        call bdotn2bmn   !  convert b.n to bmn for targeted resona
        bs=abs(bs) ! lets assume were finished with bs at this point
C        bnavg(i)=sum(bs(1:nu)) + 2*sum(bs(nu+1:nuvh-nu)) +
C     +           sum(bs(nuvh-nu+1:nuvh))
        nsurf=nsurf+1
        bnavg(nsurf)=sum(bs(1:nuv))
        bnavg(nsurf) = bnavg(nsurf) / sum(dsur(1:nuv))
        bnmax(nsurf)=maxval(bs(1:nuv)/dsur(1:nuv))
        xiotaps(nsurf)=xiotap
       endif
      enddo
c
      close(iunit)
c
c   output island sizes
c
c      open(unit=8,file="visldout."//trim(ext), status="replace")
      open(unit=8,file="visldout."//trim(coilext), status="replace")
c
      write(8,"(i5,'  / Number of surfaces')") nsurf
      write(8,"('Vacuum Island Size and Avg B.n'
     +           /' m    n    ds          drho            b.n_avg, T
     +     b.n_max, T      bmns           bmnc')")
      do isurf=1,nsurf
        do j=1,ntarget
          i=(isurf-1)*ntarget + j  ! print out all modes  per surface
          ds=sqrt(abs(ds2(i))) ! island size based on single mode
          drho=ds/2/sqrt(star(i))
          write(8,"(2i5,1p6e15.6)") mtar(i),-ntar(i), ds,drho,
     +                bnavg(isurf), bnmax(isurf), bmns(i),bmnc(i)
        enddo
c
c  calculate net island size from all modes
c
        k=(isurf-1)*ntarget + 1  !  first mode per surface
         xiotap=xiotaps(isurf)
        call analytic_island(mtar(k),bmns(k),bmnc(k),ntarget,xiotap, ds)
        drho=ds/2/sqrt(star(i))
          write(8,"(2i5,1p5e15.6)") mtar(k),-ntar(k), ds,drho,
     +                bnavg(isurf)
      enddo
c
      close(8)
c
      open(unit=8,file="visldout",status="unknown",position="append")
      write(8,"(a16,1p16e15.6)")  coilext,
     +     ( sqrt(abs(ds2(i))), i=1,nsurf*ntarget)
      close(8)
c
      open(unit=8,file="bmns",status="unknown",position="append")
      write(8,"(a16,1p6e15.6,10(/16x,1p6e15.6))")  coilext,
     +   ( bmns(i),bmnc(i), i=1,nsurf*ntarget)
      close(8)
c
      return
      end
c***********************************************************************
      subroutine readsurf(istat)
c
c                                                             09/15/99
c     modified NESCOIL subroutine
c
c
c ----------------------------------------------------------------------
      use Vmeshes
      use Vsurfdata
      use Vext
      use read_wout_mod, ONLY: np => nfp
      implicit none
c
      integer :: m,n,mr,nr,i, k,istat
      real(kind=rprec) :: crin, czin, dcrin, dczin,buin,bvin
c ----------------------------------------------------------------------
c read plasma interior surface
c
      do 30 m = 0,md
      do 30 n = -nd,nd
         cr(m,n) = 0.
         cz(m,n) = 0.
   30 continue
c
         ms    = 0
         ns     = 0
      read(iunit,*,iostat=istat) ntotal, js, s, phip, xiota, xiotap
      if(istat.ne.0)return
c      phip=phip*np ! TMP corrects Jacobian in B dot Grad phi == 1/J * dPhi/ds.
c NOTE: if s = 0., Island Width not calculated. Bmn only calculated.
         imn   = ntotal
      do  40 k=1,ntotal
      read(iunit,*) mr,nr,crin,czin, dcrin, dczin, buin,bvin
         icm(k)     = mr
         icn(k)     = nr
         cr(mr, nr) = crin
         cz(mr, nr) = czin
         dcr(mr, nr) = dcrin
         dcz(mr, nr) = dczin
         cbu(mr, nr) = buin
         cbv(mr, nr) = bvin
         ms    = max0(ms,iabs(mr))
         ns    = max0(ns,iabs(nr))
   40 continue
c
      do 90 i=0,ms
         cr (i,0) = .5*cr (i,0)
         cz (i,0) = .5*cz (i,0)
         dcr (i,0) = .5*dcr (i,0)
         dcz (i,0) = .5*dcz (i,0)
         cbu (i,0) = .5*cbu (i,0)
         cbv (i,0) = .5*cbv (i,0)
   90 continue
c
c
c read resonances to target on this surface
C NOTE: for vacisld5, no symmetry assumed.
C        Must input ntar for 360 deg, not single period
c
      read(iunit,*) ntarget  ! number of modes to target on surface
      do i=1+nts, ntarget+nts
      read(iunit, *) mtar(i), ntar(i), wtar(i), dstar(i),idstar(i)
      enddo
c
c save s as a function of nts
c
      do i=1+nts, ntarget+nts
      star(i)=s
      enddo
c
      return
      end
c***********************************************************************
      subroutine surface1
c
c                                                             01/01/89
c     modified NESCOIL subroutine
c
c
c ----------------------------------------------------------------------
      use Vmeshes
      use Vprecal
      use Vsurfdata
      implicit none
c
      integer :: i, m, n, ku, kv,isid
      real(kind=rprec) :: cofp, cofm, sifp, sifm, cm, cn,f,fx,fy,fz,pow
c ----------------------------------------------------------------------
c   outer surface   1
      do 10 i=1,nuv
         r(i)   = 0.
         z(i)   = 0.
         rs(i)  = 0.
         zs(i)  = 0.
         ru(i)  = 0.
         zu(i)  = 0.
         rv(i)  = 0.
         zv(i)  = 0.
         bu(i)  = 0.
         bv(i)  = 0.
   10 continue
      do 20 m = 0,ms
      do 20 n = 0,ns
         cm     = m*pi2
         cn     = n*pi2
         i      = 0
      do 20 kv = 1 , nv
      do 20 ku = 1 , nu
         i      = i + 1
         cofp   = comu(ku,m)*conv(kv,n)-simu(ku,m)*sinv(kv,n)
         cofm   = comu(ku,m)*conv(kv,n)+simu(ku,m)*sinv(kv,n)
         sifp   = simu(ku,m)*conv(kv,n)+comu(ku,m)*sinv(kv,n)
         sifm   = simu(ku,m)*conv(kv,n)-comu(ku,m)*sinv(kv,n)
         r(i)   = r(i)    +           cr(m,n)*cofp + cr(m,-n)*cofm
         z(i)   = z(i)    +           cz(m,n)*sifp + cz(m,-n)*sifm
         rs(i)  = rs(i)   +          dcr(m,n)*cofp + dcr(m,-n)*cofm
         zs(i)  = zs(i)   +          dcz(m,n)*sifp + dcz(m,-n)*sifm
         ru(i)  = ru(i)   - cm *    ( cr(m,n)*sifp + cr(m,-n)*sifm )
         rv(i)  = rv(i)   - cn *    ( cr(m,n)*sifp - cr(m,-n)*sifm )
         zu(i)  = zu(i)   + cm *    ( cz(m,n)*cofp + cz(m,-n)*cofm )
         zv(i)  = zv(i)   + cn *    ( cz(m,n)*cofp - cz(m,-n)*cofm )
         bu(i)  = bu(i)   +          cbu(m,n)*cofp + cbu(m,-n)*cofm
         bv(i)  = bv(i)   +          cbv(m,n)*cofp + cbv(m,-n)*cofm
   20 continue
      do 50 kv = 1 , nv
      do 50 ku = 1,nu
         i      = nu*(kv-1)+ku
         x(i)   =   coh(kv) * r(i)
         y(i)   =   sih(kv) * r(i)
         xu(i)  =   coh(kv) * ru(i)
         yu(i)  =   sih(kv) * ru(i)
         xv(i)  =   coh(kv) * rv(i) - alp*y(i)
         yv(i)  =   sih(kv) * rv(i) + alp*x(i)
         xs(i)  =   coh(kv) * rs(i)
         ys(i)  =   sih(kv) * rs(i)
   50 continue
      do 60  i = 1 , nuv
         snx(i) = yu(i)*zv(i)-zu(i)*yv(i)
         sny(i) = zu(i)*xv(i)-xu(i)*zv(i)
         snz(i) = xu(i)*yv(i)-yu(i)*xv(i)
         dsur(i)= sqrt(snx(i)*snx(i)+sny(i)*sny(i)+snz(i)*snz(i))
         ex(i) = ys(i)*zu(i)-zs(i)*yu(i)
         ey(i) = zs(i)*xu(i)-xs(i)*zu(i)
         ez(i) = xs(i)*yu(i)-ys(i)*xu(i)
CCAWB make sn a unit normal ( if solving for B dot n )
C         snx(i) = snx(i) / dsur(i)
C         sny(i) = sny(i) / dsur(i)
C         snz(i) = snz(i) / dsur(i)
c
c convert from d/du to d/dth...
c
      snx(i)=snx(i)/pi2/pi2
      sny(i)=sny(i)/pi2/pi2
      snz(i)=snz(i)/pi2/pi2
      dsur(i)=dsur(i)/pi2/pi2
      ex(i)=ex(i)/pi2
      ey(i)=ey(i)/pi2
      ez(i)=ez(i)/pi2
   60 continue
      return
      end
c***********************************************************************
      subroutine bdotn
c
c  routine evaluates B.n on interior surfaces
c
      use Vmeshes
      use coil_data_mod
      use Vsurfdata
      implicit none

      integer :: i
      real(kind=rprec) :: xp, yp, zp, pi2, bx, by, bz
      real(kind=rprec) :: f,fx,fy,fz,pow, rp, sn, ap
      integer :: isid
c
      pi2=asin(1.)*4.
      bs(1:nuv)=0.
c
      do i=1,nuv
        xp   = x(i)
        yp   = y(i)
        zp   = z(i)
        call coilfld4(xp,yp,zp, bx,by,bz)
        write(13,*) xp,yp,zp, bx,by,bz
C        call testfld(xp,yp,zp, bx,by,bz)
        bs(i)=snx(i)*bx + sny(i)*by + snz(i)*bz
CC        bp(i)=ex(i)*bx + ey(i)*by + ez(i)*bz
      bp(i)=phip ! *pi2   Not needed according to Steve - see chkwout code ( suv.x )
c      bp=sqrt(xp*xp+yp*yp)*0.125 ! b dot Grad Phi ( times J ) for test field
c      write(77,"(1p2e12.4)") bp(i),bv(i)
c      bp(i)=bp(i)+bv(i) ! add vmec field
      bsbp(i)=bs(i)/bp(i)
C      write(21,"('GRID    ',I8,8X,3f8.5)")i,x(i),y(i),z(i)
C      isid=1
C        f=sqrt(bx**2+by**2+bz**2)
C        if(f.ne.0.)then
C          pow=float(int(log10(f)+20.)-20)
C          f=10.**pow
C          fx=bx/f
C          fy=by/f
C          fz=bz/f
C        endif
C      write(21,"('FORCE   ',2I8,8X,1pe8.1,0p3f8.5)")isid,i,f,fx,fy,fz
C      isid=2
C        f=sqrt(snx(i)**2+sny(i)**2+snz(i)**2)
C        if(f.ne.0.)then
C          pow=float(int(log10(f)+20.)-20)
C          f=10.**pow
C          fx=bx/f
C          fy=by/f
C          fz=bz/f
C        endif
C      write(21,"('FORCE   ',2I8,8X,1pe8.1,0p3f8.5)")isid,i,f,fx,fy,fz
C      sn=sqrt(snx(i)**2+sny(i)**2+snz(i)**2)
C      rp=sqrt(xp*xp+yp*yp)
C      ap=sqrt( (rp-1.5)**2+zp**2)
C      write(21,"(i5,1p4e12.4)") i, rp, ap, sn, sn/rp/ap
      enddo
c
      return
      end
c***********************************************************************
      subroutine bdotn2bmn
c
      use Vmeshes
      use Vsurfdata
      use Vprecal
      implicit none
c
      integer :: i
c
      do i=1,ntarget
        nts=nts+1
       call fourtr(bsbp,nu,nv, mtar(nts),ntar(nts), bmns(nts),bmnc(nts))
        bmn(nts)=sqrt(bmns(nts)**2+bmnc(nts)**2)
        if(phip.ne.0.)then
c          ds2(nts)=bmn(nts)/phip*1.e-7/(pi2**3)/(mtar(nts)*xiotap )*16.
          ds2(nts)=bmn(nts)/(mtar(nts)*xiotap )*16.
        else
          ds2(nts)=bmn(nts)
        endif
c
      enddo
c
      return
      end
c***********************************************************************
      subroutine fourtr(b,nu,nv, m,n, bmns,bmnc)
c calculate fourier sin harmonics
c assumes m >= 0 modes retained
      use kind_spec
      implicit none
      integer :: nu, nv, m, n, i, iu, iv
      real(kind=rprec), dimension(nu*nv) :: b
      real(kind=rprec) :: bmns, bmnc, pi2, f, du, dv, u, v
c
      pi2=asin(1.)*4.
      du=1./nu
      dv=1./nv
      bmns=0.
      bmnc=0.
      i=0
      do iv=1,nv
      do iu=1,nu
      i=i+1
      u=du*(iu-1)
      v=dv*(iv-1)
      bmns=bmns+b(i)*sin((m*u+n*v)*pi2)
      bmnc=bmnc+b(i)*cos((m*u+n*v)*pi2)
      enddo
      enddo
      bmns=bmns*du*dv*2.
      bmnc=bmnc*du*dv*2.
      if(m.eq.0 .and. n.eq.0)bmns=bmns/2.
      if(m.eq.0 .and. n.eq.0)bmnc=bmnc/2.
c note: assumes routine will calculate harmonics for 0,n but not 0,-n
c       these are not independent and are combined
c       by eliminating 0,-n and doubling bmns(0,-n)
c       by eliminating 0,-n and doubling bmnc(0,-n)
c     It also assumes m < 0 modes not present.
c      ( these have been combined with m > 0 modes )
      return
      end
!c***********************************************************************
!      subroutine berr2d2(fac)
!c convert from berr to Island size ( or just [( B dot Grad s ) / ( B dot
!c  Grad phi )] )
!c
!      use Vmeshes
!      use Vfilaments
!      use Vsvddata
!      use Vsurfdata
!      use Vflags
!c
!c-----------------------
!c Up to here , we have B dot ( dR/du cross dR/dv )
!c which differs from B dot n by |dR/du cross dR/dv|
!c
!c what we are trying to target is really  [( B dot Grad s ) / ( B dot Gr
!c ad phi )]
!c  where s is the normalized flux ( Psi ) coordinate
!c  and phi is the toroidal coordinate in straight line coordinates
!c
!c
!c  Grad s  is just 1/J * ( dR/dtheta cross dR/dphi )
!c             where J is the Jacobian
!c                   R is the position vector of a point on the surface
!c       and theta and phi are the straight line coordinates
!c
!c  B dot Grad phi is taken from VMEC as just 1/J * ( dPhi/ds )
!c      ( reminder: Psi is flux, but psi is the toroidal coord )
!c   Note we are assuming B dot Grad phi does not change (much) relative
!c to B dot Grad s
!c   We neeed to do this to be able to linearize the problem.
!c
!c  Therefore [( B dot Grad s ) / ( B dot Grad phi )] is just
!c     B dot ( dR/dtheta cross dR/dphi ) / ( dPsi/ds )
!c   ( independent of J !! )
!c
!c We need to keep track of all the "constants" floating about.
!c   To convert to dR/du to dR/dtheta must divide by 2*pi
!c   To convert to dR/dv to dR/dphi   must divide again by 2*pi
!c   dPsi/ds is taken from VMEC ( phip ) which must be multiplyed by 2*pi
!c B is calculated from normalized currents and must be multiplied by Ipo
!c l
!c There is also a factor of 1e-7 for B.
!c Therefore,
!c
!c            [( B dot Grad s ) / ( B dot Grad phi )] =
!c
!c            B dot ( dR/du cross dR/dv ) *Ipol/phip*1.e-7/(2*pi)**3
!c
!c-----------------------
!c to target Island size instead of ( B dot Grad s ) / ( B dot Grad phi )
!c  we target d**2 = 16 * [( B dot Grad s ) / ( B dot Grad phi )] / ( m *
!c  xiotap )
!c
!c    where d == ds  is in normalized flux units ( alla s )
!c
!c To convert to rho=sqrt(s), drho=ds/2/sqrt(s)
!c-----------------------
!c
!      pi2=asin(1.)*4.
!c
!c
!      fac=curpol/phip*1.e-7/(pi2**3)/(mtar(nts)*xiotap )*16.
!c      fac=curpol/phip*1.e-7/(pi2**3)
!c
!      return
!      end
      subroutine testfld (xp,yp,zp, bx,by,bz)
      use testfield
      implicit none
      real(kind=rprec) :: xp, yp, zp, bx, by, bz
      real(kind=rprec) :: rp, ap, xiota, bphi, bth, bn, th, phi, br,btot
c
c routine generates a simple field pattern which
c  produces nested circular toroidal flux surfaces
c  with prescibed linear iota profile
c   also includes ideal perturbation field
c
c psuedo field
c
      rp=sqrt(xp*xp+yp*yp)
      ap=sqrt( (rp-rc)**2 + zp**2 )
c
      xiota=xiota0 + (xiota1-xiota0)*ap/ac
c
      phi=atan2(yp,xp)
      if(rp-rc .eq. 0. .and. zp .eq. 0.)then
            th=0.
      else
            th=atan2(zp,rp-rc)
      endif
c
C      bphi=rp
C      bth=ap*xiota
C      bn=0.
C      if(ap.gt.ac*.1)bn=bmnx*sin(mx*th+nx*phi)
      bphi=rc/rp
      bth=bphi*ap/rp*xiota
      bn=0.
      if(ap.gt.ac*.1)bn=bmnx*sin(mx*th+nx*phi)*bphi/rp
c
      br=-bth*sin(th)+bn*cos(th)
      bz=bth*cos(th)+bn*sin(th)
c
      bx=br*cos(phi) - bphi*sin(phi)
      by=br*sin(phi) + bphi*cos(phi)
      btot=sqrt(bx*bx+by*by+bz*bz)
c
      bx=bx*1.e7
      by=by*1.e7
      bz=bz*1.e7


c
      return
      end
c*************************************************
      subroutine analytic_island(m,bms,bmc,nm, xiotap, ds)
c routine constructs ( 'traces' ) islands based on shear and resonant perturbations
      use kind_spec
      implicit none
c
      integer i, j, k, nm, nsurf, nth
      real(kind=rprec) :: pi2,xiotap,ds,smax,s,s2,th,th0

      real(kind=rprec), dimension(nm) ::  bms, bmc
      integer, dimension(nm) ::  m

      nth=100
      nsurf=100
      pi2=asin(1.)*4.
c
      smax=0.
      do j=0,nsurf
        th0=pi2/m(1)*j/nsurf
        do i=0,nth
          th=pi2/m(1)*i/nth
          s2=0.
          do k=1,nm
            s2=s2+2.*bms(k)/m(k)/xiotap*(cos(m(k)*th0)-cos(m(k)*th))
            s2=s2-2.*bmc(k)/m(k)/xiotap*(sin(m(k)*th0)-sin(m(k)*th))
          enddo
          s=0.
          if(s2.gt.0.)then
            s=sqrt(s2)
            smax=max(smax,s)
          endif
        enddo
      enddo
c
      ds=smax*2.
      return
      end