*
      subroutine noneqr(tempr,tau,n,nzmax,nzpmax,nionp,vr,vl,eig,feqs,
     &                  fs,fspec,work,lt,fout,ionel,ionstage)

      integer n, nzmax, nzpmax, nionp
      real tau(n),tempr(n)
      real fout(nionp)
      integer ionel(nionp), ionstage(nionp)
      double precision vr(nzmax,nzmax),vl(nzmax,nzmax),
     .       eig(nzmax), feqs(0:nzmax),
     .       fs(0:nzmax),fspec(nzmax),work(nzmax)
      integer lt(n)

*
*        Calculates ionization fractions fout at electron temperature
*        tempr(n) and ionization parameter tau(n), for electron 
*        temperatures tempr given in a tabular form as a function of 
*        ionization parameter tau. 
*
*        INPUT:  tempr  -  electron temperature (K)
*                tau    -  ionization parameter (integral of electron
*                          density over time) (s cm^-3) 
*                n      -  number of elements in TEMPR and TAU (must be
*                          larger than 1)
*                nzmax    - largest atomic number of elements
*                nzpmax   - nzmax + 1
*                nionp    - total number of ions
*        WORK:   vr       - the matrix of right eigenvectors
*                vl       - the matrix of left eigenvectors
*                eig      - the set of eigenvalues
*                feqb     - equilibrium ionization fractions
*                feqs     - ionization fractions
*                fs       - equilibrium ionization fractions
*                fspec    - ionization fractions with equilibrium subtracted
*                work     - work array
*                lt       - indices for the temperatures to read
*        OUTPUT: fout   -  ionization fractions
*                ionel    - the element for each fout entry
*                ionstage - the ion stage for each fout entry
*                           (ionel+1=fully stripped)
*
*        Note that first element of tau should be zero.
*
        integer i, j, k
        integer max, min
        real alog10
c
c	these are parameters for the ionization states of the elements
c		nz = z of ion
c		nzp = z + 1
c		nzm = z - 1
c		nelem = total number of elements
c		ielem = current element
        integer nz,nzm,nzp,nelem,ielem
        parameter(nelem=11)
        integer nzz(nelem),iunit
c
c	the lower and upper limiting values of log temperature,
c		the increment in log temp., and the number of 
c		temperature values in the ionization tables.
c
        real tll,tlu,dtl
        integer nt
        parameter (tll=4.0e0,tlu=9.0e0,dtl=4.0e-3,nt=1251)
c
        double precision sum, dble, dexp 

c	index for temperature value in ionization files
        integer it
        integer ion, ionout, nret

      logical qanyf, q0502
      integer status, ioff, block, hdutyp, lenn

      character pathd*255, filenm*255, cz(nelem)*2
      character pname*128, pvalue*128, neivers*128

      character fgmodf*255, fgmstr*128
      integer lenact
      external lenact, fgmodf

      save nzz, cz, neivers

c	the z of each element
      data nzz/6,7,8,10,12,14,16,18,20,26,28/
      data cz/'c ','n ','o ','ne','mg','si','s ','ar','ca','fe','ni'/
      data neivers / '1.1' /

      pathd = fgmodf()
      lenn = lenact(pathd)

      pname = 'NEIVERS'
      pvalue = fgmstr(pname)
      IF ( pvalue .NE. ' ' ) neivers = pvalue
      q0502 = .TRUE.
      IF ( neivers .EQ. '1.0' ) q0502 = .FALSE.

      status = 0


        fout(1) = 1.0
        fout(2) = 1.0
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c	calculate temperature index
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      do 300 k=2,n
         it=(alog10(0.5*(tempr(k)+tempr(k-1)))-tll)/dtl+1.5
         it=max(1,it)
         lt(k-1)=min(nt,it)
300   continue
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c	loop over elements
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
           ionout=2
           do 1500 ielem=1,nelem
              nz=nzz(ielem)
              nzm=nz-1
              nzp=nz+1
              CALL getlun(iunit)
              if ( q0502 ) then
                 filenm = pathd(:lenn) // 'eigen' // 
     &              cz(ielem)(:lenact(cz(ielem))) // '0502.fits'
              else
                 filenm = pathd(:lenn) // 'eigen' // 
     &              cz(ielem)(:lenact(cz(ielem))) // '.fits'
              endif
              CALL ftopen(iunit, filenm, 0, block, status)
              CALL ftmrhd(iunit, 1, hdutyp, status)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c	start neutral
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
              fs(0)=1.0d0
              do 400 i=1,nz
400           fs(i)=0.0d0
              do 1000 k=1,n-1
                it=lt(k)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c	read in data for correct temperature
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
                CALL ftgcvd(iunit, 1, it, 1, nzp, 0.d0, feqs(0), 
     &                      qanyf, status)
                CALL ftgcvd(iunit, 2, it, 1, nz, 0.d0, eig, 
     &                      qanyf, status)
                ioff = 1
                DO j = 1, nz
                   CALL ftgcvd(iunit, 3, it, ioff, nz, 0.d0, vr(1,j), 
     &                         qanyf, status)
                   CALL ftgcvd(iunit, 4, it, ioff, nz, 0.d0, vl(1,j), 
     &                         qanyf, status)
                   ioff = ioff + nz
                ENDDO
cc                read(iunit,rec=it) (feq(i),i=1,nzp),(eig(i),i=1,nz),
cc     1            ((vr(i,j),i=1,nz),j=1,nz),((vl(i,j),i=1,nz),j=1,nz)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c	decompose into eigenvectors
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
                do 500 i=1,nz
500             work(i)=fs(i)-feqs(i)
                do 700 i=1,nz
                   fspec(i)=0.0d0
                   do 600 j=1,nz
600                   fspec(i)=fspec(i)+vl(j,i)*work(j)
700             continue
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c	propagate forward
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
                do 800 i=1,nz
800                work(i)=fspec(i)*dexp(eig(i)*dble((tau(k+1)-tau(k))))
                sum=0.0d0
                do 900 i=1,nz
                   fs(i)=0.0d0
                   do 850 j=1,nz
                      fs(i)=fs(i)+work(j)*vr(i,j)
850                continue
                   fs(i)=fs(i)+feqs(i)
                   if (fs(i).lt.0.0d0) fs(i)=0.0d0
900                sum=sum+fs(i)
                fs(0)=1.0d0-sum
                if (fs(0).lt.0.0d0) fs(0)=0.0d0
1000          continue
              CALL ftclos(iunit, status)
              CALL frelun(iunit)

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c	store ion fractions in output array
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
           do 1100 i=1,nzp
              ionel(ionout+i) = nz
              ionstage(ionout+i) = i
              fout(ionout+i)=fs(i-1)
 1100      continue
           ionout=ionout+nzp
1500       continue


c set ionel and ionstage for H and He, which are assumed completely ionized

      ionel(1) = 1
      ionstage(1) = 2
      ionel(2) = 2
      ionstage(2) = 3

      return
      end
*

