
      SUBROUTINE oxneispec(ntemp, temp, nion, ionfrac, ionel, ionstage, 
     &                   ne, ear, z, linesh, nabun, abund, abunel, 
     &                   qcalc, callr, photar, photsum, workphot,
     &                   earold)

      INTEGER ne, nion, nabun, ntemp

      REAL temp(ntemp), ionfrac(nion,ntemp), ear(0:ne), z, linesh(2)
      REAL abund(nabun), photar(ne)
      REAL photsum(ne, nabun), workphot(ne, nabun), earold(0:ne)

      INTEGER ionel(nion), ionstage(nion), abunel(nabun)

      CHARACTER*(*) callr

      LOGICAL qcalc

c Routine to return the spectrum from a multi-temperature NEI plasma. 
c Initial version is simply a wrap-up of Kazik's code and doesn't use 
c the more sophisticated input arrays. NB The first entry in the abund 
c array should be for H.

c Arguments
c    temp         r        i: temperature (keV)
c    nion         i        i: size of ion arrays (may be > actual no. of ions)
c    ionfrac      r        i: fraction in ion
c    ionel        i        i: element for this ion
c    ionstage     i        i: level for this ion
c    ne           i        i: number of response energies
c    ear          r        i: response energies
c    z            r        i: redshift
c    linesh       r        i: lineshape parameters
c    nabun        i        i: number of elements with abundances
c    abund        r        i: abundances
c    abunel       i        i: element for this abundance
c    qcalc        l        i: if true force a recalculation of all arrays
c    callr        c        i: name of calling routine
c    photar       r        r: output spectrum
c    photsum      r        w: workspace array for accumulating spectrum
c    workphot     r        w: workspace array for spectrum
c    oldear       r        w: workspace array to save energies

      REAL solab(12)

      DOUBLE PRECISION pj
      REAL zf

      INTEGER itemp, i, j

      LOGICAL qear

c Solar abundances as defined by Anders and Grevesse.

      DATA solab /12.00, 10.99, 8.56, 8.05, 8.93, 8.09, 7.58, 7.55, 
     &            7.21, 6.36, 7.67, 6.25 /

c  Check whether the energy array has changed since the last invocation.
c  If it has we force a recalculation.

      qear = .FALSE.
      IF ( earold(0) .LE. 0 ) THEN
         qear = .TRUE.
      ELSE
         DO i = 0, ne
            IF ( ear(i) .NE. earold(i) ) qear = .TRUE.
         ENDDO
      ENDIF
      IF ( qear ) THEN
         DO i = 0, ne
            earold(i) = ear(i)
         ENDDO
      ENDIF
      qcalc = qcalc .OR. qear

c  Initialize summation array

      DO j = 1, nabun
         DO i = 1, ne
            photsum(i,j) = 0.d0
         ENDDO
      ENDDO

      DO itemp = 1, ntemp

         CALL nei1tspec(temp(itemp), nion, ionfrac(1,itemp), ionel, 
     &                  ionstage, ne, ear, z, linesh, qcalc, 
     &                  callr, nabun, workphot)

         DO j = 1, nabun
            DO i = 1, ne
               photsum(i,j) = photsum(i,j) + DBLE(workphot(i,j))
            ENDDO
         ENDDO
         
      ENDDO

c Sum all elements weighting by abundance

      DO i = 1, ne
         photar(i) = 0.0
      ENDDO

      DO j = 1, nabun
         pj = DBLE(abund(j))*10**(solab(j)-12.00)
         DO i = 1, ne
            photar(i)=photar(i)+SNGL(pj*photsum(i,j))
         ENDDO
      ENDDO

* Correct for redshift.

      zf = 1./(1+z)
      DO i = 1, ne
         photar(i) = 1.e14*photar(i)*zf
      ENDDO

      RETURN
      END

c **********************************************************************

      SUBROUTINE nei1tspec(temp, nion, ionfrac, ionel, ionstage, ne, 
     &                     ear, z, linesh, qcalc, callr, nabun, 
     &                     workphot)

      INTEGER ne, nion, nabun

      REAL temp, ionfrac(nion), ear(0:ne), z, linesh(2)
      REAL workphot(ne, nabun)

      INTEGER ionel(nion), ionstage(nion)

      CHARACTER*(*) callr

      LOGICAL qcalc

c Routine to return the spectrum from a single-temperature NEI plasma. 
c Initial version is simply a wrap-up of Kazik's code and doesn't use 
c the more sophisticated input arrays.

c Arguments
c    temp         r        i: temperature (keV)
c    nion         i        i: size of ion arrays (may be > actual no. of ions)
c    ionfrac      r        i: fraction in ion
c    ionel        i        i: element for this ion
c    ionstage     i        i: level for this ion
c    ne           i        i: number of response energies
c    ear          r        i: response energies
c    z            r        i: redshift
c    linesh       r        i: lineshape parameters
c    qcalc        l        i: if true force a recalculation of all arrays
c    callr        c        i: name of calling routine
c    nabun        i        i: number of elements with abundances
c    workphot     r        r: output spectrum for each element

c
      real t, tprev
      real zf, zprev
      character*20 cprev
*        Conversion factors.
      real keV, boltz, tconv
      parameter ( keV = 1.60219e-9, boltz = 1.38062e-16 )
      parameter ( tconv = keV/boltz )

*
*          nelem         -  number of heavy elements
*          ielem         -  curent index of elements (2 for He^++,
*                           3-12 for heavy elements
*          istart, iend  -  dummy integer indexes
      integer nelem, nelemp
      parameter (nelem=10,nelemp=nelem+2)
      integer nlines(3:nelemp), ielem, istart, iend
*        nline - number of lines, el(iline) - energy of line iline in keV
      integer nline, iline
      parameter (nline=566)
      real el(nline)
c
c			xline(iline) = line luminosity of line iline
c	xline(iline) = lines
c
      real xline(nline)
      integer nmaxfe
      parameter ( nmaxfe = 146000 )
      real elfe(nmaxfe), xlinefe(nmaxfe)
      integer nfe
      integer ifreq, ibine, l, j
      common /elin/ el

      save nlines, xline, elfe, xlinefe, nfe
      save tprev, zprev, cprev

      data nlines / 15, 15, 15, 20, 22, 48, 70, 71, 157, 133 /
      data tprev, zprev / 2*0. /
      data cprev / ' ' /


      t = tconv*temp
      zf = 1./(1. + z)

      if( t .ne. tprev .or. z .ne. zprev .or. qcalc .or.
     &    callr .ne. cprev )then

* Line emission.
         call emisl(t,ionfrac,xline)
* Ne- to Li-like Fe ions
         call felines(t,ionfrac(104),elfe,xlinefe,nfe)

*           Continua.
         DO ielem = 1, nelemp
            DO ifreq = 1, ne
               workphot(ifreq,ielem) = 0.0
            ENDDO
         ENDDO
         call continua(t,ionfrac,ear,ne,z,workphot)

*           Add lines into energy bins.
         istart = 1
         iend = 0
         do 2000 ielem=3,nelemp
           iend = iend + nlines(ielem)
	   do 1200 iline=istart,iend           
              call hunt(zf*el(iline),ear,ne,ifreq)
	      if( ifreq .ne. 0) workphot(ifreq,ielem) =
     .	      workphot(ifreq,ielem) + xline(iline)/(keV*el(iline))
1200       continue
           istart = iend + 1
 2000    continue

*           Add Fe L-shell lines.
        do 3200 iline=1,nfe
           call hunt(zf*elfe(iline),ear,ne,ifreq)
	   if( ifreq .ne. 0) workphot(ifreq,11) = 
     .	      workphot(ifreq,11) + xlinefe(iline)/(keV*elfe(iline))
 3200    continue

         tprev = t
         zprev = z
         cprev = callr
      endif

      RETURN
      END
