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

      IMPLICIT NONE

      INTEGER ne, nion, nabun, ntemp

      REAL temp(ntemp), ionfrac(nion,ntemp), ear(0:ne), z, linesh(2)
      REAL abund(nabun), photar(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 This is a wrap-up routine to call either nneispec (the default) or
c oneispec if NEIVERS is set to 1.x.
c 

c Arguments
c    temp         r        i: temperature (keV)
c    nion         i        i: number of different ions
c    ionfrac      r        i: fraction in ion
c    ionel        i        i: element for this ion
c    ionstage     i        i: level for this ion (ionel+1 = fully stripped)
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

      INTEGER TOTEL
      PARAMETER(TOTEL=30)

      REAL andgrv(TOTEL)

      INTEGER iold, newnion, newnabun, j

      CHARACTER pname*128, pvalue*128, ovalue*128, tmpstr*255
      CHARACTER*2 celt(TOTEL)

      LOGICAL qdone

      INTEGER lenact
      REAL fgabnd
      CHARACTER fgmstr*128
      EXTERNAL lenact, fgmstr, fgabnd

      SAVE newnion, newnabun, ovalue

      DATA celt/'H ', 'He', 'Li', 'Be', 'B ', 'C ', 'N ', 'O ', 'F ', 
     &          'Ne', 'Na', 'Mg', 'Al', 'Si', 'P ', 'S ', 'Cl', 'Ar', 
     &          'K ', 'Ca', 'Sc', 'Ti', 'V ', 'Cr', 'Mn', 'Fe', 'Co', 
     &          'Ni', 'Cu', 'Zn'/
      DATA andgrv/1.00e+00, 9.77e-02, 1.45e-11, 1.41e-11, 3.98e-10, 
     &            3.63e-04, 1.12e-04, 8.51e-04, 3.63e-08, 1.23e-04, 
     &            2.14e-06, 3.80e-05, 2.95e-06, 3.55e-05, 2.82e-07, 
     &            1.62e-05, 3.16e-07, 3.63e-06, 1.32e-07, 2.29e-06, 
     &            1.26e-09, 9.77e-08, 1.00e-08, 4.68e-07, 2.45e-07, 
     &            4.68e-05, 8.32e-08, 1.78e-06, 1.62e-08, 3.98e-08/
      DATA ovalue / '1.1' /


c Check for the NEIVERS string parameter. 

      pname = 'NEIVERS'
      pvalue = fgmstr(pname)
      IF ( pvalue .NE. ' ' .AND. pvalue .NE. ovalue ) THEN
         ovalue = pvalue
         tmpstr = 'Changing to NEI models version '//
     &            pvalue(:lenact(pvalue))
         CALL xwrite(tmpstr,10)
      ENDIF

      tmpstr = 'Running NEI models version '//ovalue(:lenact(ovalue))
      CALL xwrite(tmpstr, 20)

c Convert from whatever the current abundance ratios are to the Anders &
c Grevesse ratios which are used internally

      DO j = 1, nabun
         abund(j) = abund(j) * fgabnd(celt(abunel(j))) 
     &                       / andgrv(abunel(j))
      ENDDO

      IF ( ovalue(1:1) .EQ. '2' ) THEN

         CALL nneispec(ntemp, temp, nion, ionfrac, ionel, ionstage, 
     &                 ne, ear, z, linesh, nabun, abund, abunel, 
     &                 qcalc, photar)

      ELSEIF ( ovalue(1:1) .EQ. '1' ) THEN

c The old version which had no Ar so we have to strip this out of the ion*
c and ab* arrays (they may already have been removed from the ion* arrays
c if they have not been recalculated since the last invocation)

         qdone = .TRUE.
         DO iold = 1, nion
            IF ( ionel(iold) .EQ. 18 ) qdone = .FALSE.
         ENDDO
         IF ( .NOT.qdone ) THEN
            newnion = 0
            DO iold = 1, nion
               IF ( ionel(iold) .NE. 18 ) THEN
                  newnion = newnion + 1
                  DO j = 1, ntemp
                     ionfrac(newnion, j) = ionfrac(iold, j)
                  ENDDO
                  ionel(newnion) = ionel(iold)
                  ionstage(newnion) = ionstage(iold)
               ENDIF
            ENDDO
         ENDIF

         newnabun = 0
         DO iold = 1, nabun
            IF ( abunel(iold) .NE. 18 ) THEN
               newnabun = newnabun + 1
               abund(newnabun) = abund(iold)
               abunel(newnabun) = abunel(iold)
            ENDIF
         ENDDO

         CALL oneispec(ntemp, temp, nion, ionfrac, ionel, ionstage, 
     &                 ne, ear, z, linesh, newnabun, abund, abunel, 
     &                 qcalc, callr, photar)

      ENDIF

      RETURN
      END

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

      SUBROUTINE xneispec(ntemp, temp, nion, ionfrac, ionel, ionstage, 
     &                    ne, ear, z, linesh, nabun, abund, abunel, 
     &                    qcalc, photar, nmtval, tval, nmions, 
     &                    cont, cenerg, conEl, conIon, conpt, 
     &                    consiz, pseudo, penerg, psupt, psusiz, lineE,
     &                    lineEm, lineEl, lineIon, linpt, linsiz, 
     &                    qtherm, velocity, status)

      IMPLICIT NONE

      INTEGER ne, nion, nabun, ntemp, nmtval, nmions, status

      REAL temp(ntemp), ionfrac(nion,ntemp), ear(0:ne), z, linesh(2)
      REAL abund(nabun), photar(ne)
      REAL tval(nmtval), velocity
      REAL cont(*), cenerg(*), pseudo(*), penerg(*)
      REAL lineE(*), lineEm(*)

      INTEGER lineEl(*), lineIon(*)
      INTEGER conpt(nmions, nmtval), consiz(nmions, nmtval)
      INTEGER psupt(nmions, nmtval), psusiz(nmions, nmtval)
      INTEGER conEl(nmions, nmtval), conIon(nmions, nmtval)
      INTEGER linpt(nmtval), linsiz(nmtval)
      INTEGER ionel(nion), ionstage(nion), abunel(nabun)

      LOGICAL qcalc, qtherm

c Routine to return the spectrum from a multi-temperature NEI plasma. 

c Arguments
c    ntemp        i        i: number of temperatures
c    temp         r        i: temperature (keV)
c    nion         i        i: number of different ions
c    ionfrac      r        i: fraction in ion
c    ionel        i        i: element for this ion
c    ionstage     i        i: level for this ion (ionel+1 = fully stripped)
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    photar       r        r: output spectrum
c    nmtval       I        i: number of tabulated temperatures
C    tval         R        i: tabulated temperatures
C    nmions       I        i: number of tabulated ions
C    cont         R        i: Tabulated continua
C    cenerg       R        i: Energies for tabulated continua
C    conEl        I        i: Atomic numbers for tabulated continua
c    conIon       I        i: Ion stages for tabulated continua
C    conpt        I        i: Offsets in the cont array for each ion and kT
C    consiz       I        i: Sizes in the cont array for each ion and kT
C    pseudo       R        i: Tabulated pseudo-continua
C    penerg       R        i: Energies for tabulated pseudo-continua
C    psupt        I        i: Offsets in the pseudo array for each ion and kT
C    psusiz       I        i: Sizes in the pseudo array for each ion and kT
C    lineE        R        i: Tabulated line energies
C    lineEm       R        i: Tabulated line emissivities
C    lineEl       I        i: Tabulated line elements
c    lineIon      I        i: Tabulated line ion
C    linpt        I        i: Offsets in the line arrays for each kT
C    linsiz       I        i: Sizes in the line arrays for each kT
C    qtherm       L        i: If true then apply thermal broadening to lines
C    velocity     R        i: Gaussian velocity broadening to apply
C    status       I        r: 0==OK

c This code is adapted from SUMAPE but loops over ions rather than elements.

      REAL c_cm_s, keV_erg, amu_g, km_cm
c Speed of light in cm s^-1
      PARAMETER (c_cm_s=2.9979246e10)
c kT in erg corresponding to temperature of 1 keV
      PARAMETER (keV_erg=1.60217646e-9)
c Unified atomic mass constant in g
      PARAMETER (amu_g=1.660538e-24)
c km in cm
      PARAMETER (km_cm=1e5)

      INTEGER NUMZ
      PARAMETER (NUMZ=28)

      REAL abunZ(NUMZ), atomM(NUMZ)
      REAL tdelta, tcoeff(2), energy, edelta, zf, coeff
      REAL width, lineflux

      INTEGER atomZ(13)
      INTEGER itemp(2), lim(2)
      INTEGER ien, it, ielt, itabt, il, ihigh, ilow, ion
      INTEGER ionout, i, itpt

      CHARACTER wrtstr*72

      SAVE atomZ, atomM

      DATA atomZ /2, 6, 7, 8, 10, 12, 13, 14, 16, 18, 20, 26, 28/

c Atomic masses in AMU
      DATA atomM /1.00794, 4.002602, 6.941, 9.012182, 10.811,
     &            12.0107, 14.0067, 15.9994, 18.9984032, 20.1797,
     &            22.989770, 24.3050, 26.981538, 28.0855, 30.973761,
     &            32.065, 35.453, 39.948, 39.0983, 40.078,
     &            44.955910, 47.867, 50.9415, 51.9961, 54.938049,
     &            55.845, 58.933200, 58.6934/

      status = 0

c Initialize output array

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

c set up the abundances in terms of atomic number so we can match them
c to the lineEl array

      DO ielt = 1, NUMZ
         abunZ(ielt) = 0.0
      ENDDO
      DO ielt = 1, nabun
         abunZ(abunel(ielt)) = abund(ielt)
      ENDDO

c Loop over required temperatures

      DO it = 1, ntemp

c Find tabulated temperatures to be used in the interpolation and
c calculate interpolation coefficients

         IF ( temp(it) .LT. tval(1) ) THEN
            WRITE(wrtstr, '(a,g12.5)')
     &        ' Warning: temperature < minimum in model file: ', tval(1)
            CALL xwrite(wrtstr, 10)
            itemp(1) = 1
            itemp(2) = 1
         ELSEIF ( temp(it) .GT. tval(nmtval) ) THEN
            WRITE (wrtstr, '(a,g12.5)')
     &        ' Warning: temperature > maximum in model file: ', 
     &        tval(nmtval)
            CALL xwrite(wrtstr, 10)
            itemp(1) = nmtval
            itemp(2) = nmtval
         ELSE
            itemp(2) = 1
            DO WHILE ( temp(it) .GE. tval(itemp(2)) )
               itemp(2) = itemp(2) + 1
            ENDDO
            itemp(1) = itemp(2) - 1
            IF ( temp(it) .EQ. tval(itemp(1)) ) itemp(2) = itemp(1)
         ENDIF

         tdelta = tval(itemp(2)) - tval(itemp(1))
         IF ( tdelta .GT. 0. ) THEN
            tcoeff(1) = (tval(itemp(2)) - temp(it))/tdelta
            tcoeff(2) = (temp(it) - tval(itemp(1)))/tdelta
         ELSE
            tcoeff(1) = 1.0
            tcoeff(2) = 0.0
         ENDIF

         WRITE(wrtstr, '(a, f8.4, a, f8.4, a, f8.4)' ) 
     &       'Interpolating between ', tval(itemp(1)), 
     &       ' and ', tval(itemp(2)), ' for ', temp(it)
         CALL xwrite(wrtstr, 15)
         WRITE(wrtstr, '(a, f8.4, a, f8.4)' )
     &       'with coefficients ', tcoeff(1), ' and ', tcoeff(2)
         CALL xwrite(wrtstr, 15)

c Loop over the two temperatures to use in interpolation

         DO itabt = 1, 2

            itpt = itemp(itabt)

c Loop over tabulated ions. There are at most nmions - if there are less
c then the arrays are set to -1.

            DO ion = 1, nmions

               IF ( conEl(ion,itpt) .EQ. -1 ) GOTO 100

               coeff = tcoeff(itabt) * abunZ(conEl(ion,itpt))

c If we are not using this element (ie abunZ is zero then jump to end
c of loop

               IF ( coeff .EQ. 0. ) GOTO 100

c Find which entry in the ionfrac array this ion matches

               ionout = 0
               DO i = 1, nion
                 IF ( (ionel(i) .EQ. conEl(ion,itpt)) .AND.
     &                (ionstage(i) .EQ. conIon(ion,itpt)) ) ionout = i
               ENDDO

               IF ( ionout .EQ. 0 ) GOTO 100

               coeff = coeff * ionfrac(ionout,it)

c Interpolate the continuum onto the output energy array. Note that the
c input is flux/keV values at the energies tabulated so need to multiply
c by the bin width to convert to the units that XSPEC requires. Also note
c that this algorithm assumes that the response energies are monotonically
c increasing.

               lim(1) = conpt(ion, itpt)
               lim(2) = lim(1) + consiz(ion, itpt) - 1

               ihigh = lim(1)

               DO ien = 1, ne

                  energy = ( (ear(ien-1) + ear(ien))/2. ) * (1 + z)

                  IF ( energy .GE. cenerg(lim(1)) .AND.
     &                 energy .LE. cenerg(lim(2)) ) THEN

                     DO WHILE ( energy .GT. cenerg(ihigh) .AND. 
     &                          ihigh .LE. lim(2) )
                        ihigh = ihigh + 1
                     ENDDO
                     ilow = ihigh - 1

                     photar(ien) = photar(ien) + coeff *
     &                       ( (cenerg(ihigh)-energy)*cont(ilow)
     &                        +(energy-cenerg(ilow))*cont(ihigh) )
     &                       * ( ear(ien) - ear(ien-1) ) / (1 + z)
     &                       / (cenerg(ihigh)-cenerg(ilow))


                  ENDIF

               ENDDO

c Repeat for the pseudo-continuum

               lim(1) = psupt(ion, itpt)
               lim(2) = lim(1) + psusiz(ion, itpt) - 1

               ihigh = lim(1)

               DO ien = 1, ne

                  energy = ( (ear(ien-1) + ear(ien))/2. ) * (1 + z)

                  IF ( energy .GE. penerg(lim(1)) .AND.
     &                 energy .LE. penerg(lim(2)) ) THEN

                     DO WHILE ( energy .GT. penerg(ihigh) .AND. 
     &                          ihigh .LE. lim(2) )
                        ihigh = ihigh + 1
                     ENDDO
                     ilow = ihigh - 1

                     photar(ien) = photar(ien) + coeff *
     &                       ( (penerg(ihigh)-energy)*pseudo(ilow)
     &                        +(energy-penerg(ilow))*pseudo(ihigh) )
     &                       * ( ear(ien) - ear(ien-1) ) / (1 + z)
     &                       / (penerg(ihigh)-penerg(ilow))

                  ENDIF

               ENDDO

c End loop over elements.

 100           CONTINUE

            ENDDO

c Now do the lines. Loop over the lines for this temperature

            DO il = linpt(itpt), linpt(itpt)+linsiz(itpt)-1

               IF ( abunZ(lineEl(il)) .NE. 0. ) THEN

                  energy = lineE(il) / (1 + z)

c Find the width of the line using the equation
c  W = E/c  sqrt [ kT / (A mp ) + velocity^2 ] if thermal broadening
c is on else
c  W = (E/c) velocity

                  IF ( qtherm ) THEN
                     width = energy / c_cm_s *
     &                 sqrt( temp(it) * keV_erg /
     &                 atomM( lineEl(il) ) / amu_g + 
     &                 velocity*velocity*km_cm*km_cm)
                  ELSE
                     width = (energy / c_cm_s) * velocity * km_cm
                  ENDIF

                  IF ( energy-3*width .LE. ear(ne) .AND.
     &                 energy+3*width .GE. ear(0) ) THEN

c Find which entry in the ionfrac array this line is from

                     ionout = 0
                     DO i = 1, nion
                        IF ( (ionel(i) .EQ. lineEl(il)) .AND.
     &                      (ionstage(i) .EQ. lineIon(il)) ) ionout = i
                     ENDDO

c Add in multiplying by abundance and ion fraction

                     IF ( ionout .GT. 0 ) THEN

                        lineflux = tcoeff(itabt) * abunZ(lineEl(il)) 
     &                            * lineEm(il) * ionfrac(ionout,it)

                        CALL gaussline(ear, ne, energy, width, 
     &                                 lineflux, 3.0, .TRUE., photar)

                     ENDIF

                  ENDIF

               ENDIF

c End loop over lines.

            ENDDO

c End loop over two temperatures used in the interpolation

         ENDDO

c End loop over temperatures.

      ENDDO

c Fix the normalization by multiplying by 1e14

      DO ien = 1, ne
         photar(ien) = photar(ien) * 1.e14
      ENDDO

c Correct for time dilation

      IF ( z .NE. 0 ) THEN
         zf = 1./(1.+z)
         DO ien = 1, ne
            photar(ien) = photar(ien) * zf
         ENDDO
      ENDIF

      RETURN
      END



