
      SUBROUTINE sumape(ear, ne, abun, z, ninputt, inputt, dem, 
     &     nmtval, tval, nmelt, cont, cenerg, conEl, conpt, 
     &     consiz, pseudo, penerg, psupt, psusiz, lineE,
     &     lineEm, lineEl, linpt, linsiz, qtherm, velocity,
     &     photar, status)

      IMPLICIT NONE

      INTEGER ne, ninputt, nmtval, nmelt, status
      REAL z, velocity
      REAL ear(0:*), abun(*), inputt(*), dem(*), photar(*)
      REAL tval(nmtval)
      REAL cont(*), cenerg(*), pseudo(*), penerg(*)
      REAL lineE(*), lineEm(*)
      INTEGER lineEl(*)
      INTEGER conpt(nmelt, nmtval), consiz(nmelt, nmtval)
      INTEGER psupt(nmelt, nmtval), psusiz(nmelt, nmtval)
      INTEGER conEl(nmelt)
      INTEGER linpt(nmtval), linsiz(nmtval)
      LOGICAL qtherm

C Subroutine to interpolate and sum the emission from a plasma with a
C DEM given by the input array. This routine actually does the work for
C the case when the tabulated data is in the APED format of separated
C continuum, pseudo-continuum, and lines.
C Argument :
c       ear       R       i: model energy ranges 
c       ne        I       i: number of model energies 
c       abun      R       i: abundances 
c       z         R       i: redshift 
c       ninputt   I       i: number of temperatures 
c       inputt    R       i: temperatures 
c       dem       R       i: emission measures for input temperatures 
c       nmtval    I       i: number of tabulated temperatures
C       tval      R       i: tabulated temperatures
C       nmelt     I       i: number of tabulated elements
C       cont      R       i: Tabulated continua
C       cenerg    R       i: Energies for tabulated continua
C       conEl     I       i: Atomic numbers for tabulated elements
C       conpt     I       i: Offsets in the cont array for each Z and kT
C       consiz    I       i: Sizes in the cont array for each Z 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 Z and kT
C       psusiz    I       i: Sizes in the pseudo array for each Z and kT
C       lineE     R       i: Tabulated line energies
C       lineEm    R       i: Tabulated line emissivities
C       lineEl    I       i: Tabulated line elements
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       photar    R       r: spectrum 
C       status    I       r: 0==OK

      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, zf, coeff, width
      REAL lineflux, z1, z15

      INTEGER limup, limdown, itemp(2)
      INTEGER ien, it, ielt, itabt, il, ihigh, ilow
      INTEGER atomZ(13)

      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

      abunZ(1) = 1.
      DO ielt = 2, NUMZ
         abunZ(ielt) = 0.
      ENDDO
      DO ielt = 1, 13
         abunZ(atomZ(ielt)) = abun(ielt)
      ENDDO

c Precalculate for speed
      z1 = 1. + z
      z15 = 0.5*z1

c Loop over DEM temperatures
      DO it = 1, ninputt

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

         IF ( inputt(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 ( inputt(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 ( inputt(it) .GE. tval(itemp(2)) )
               itemp(2) = itemp(2) + 1
            ENDDO
            itemp(1) = itemp(2) - 1
            IF ( inputt(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)) - inputt(it))/tdelta
            tcoeff(2) = (inputt(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 ', inputt(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

c Loop over elements

            DO ielt = 1, nmelt

               coeff = tcoeff(itabt) * dem(it) * abunZ(conEl(ielt)) * z1

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.

               limdown = conpt(ielt, itemp(itabt))
               limup = limdown + consiz(ielt, itemp(itabt)) - 1

               ihigh = limdown

               DO ien = 1, ne

                  energy = (ear(ien-1) + ear(ien))*z15

                  IF ( energy .GE. cenerg(limdown) .AND.
     &                 energy .LE. cenerg(limup) ) THEN

                     DO WHILE ( energy .GT. cenerg(ihigh) .AND. 
     &                          ihigh .LE. limup )
                        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) )
     &                       / (cenerg(ihigh)-cenerg(ilow))

                  ENDIF

               ENDDO

c Repeat for the pseudo-continuum

               IF ( psusiz(ielt, itemp(itabt)) .GT. 0 ) THEN

                  limdown = psupt(ielt, itemp(itabt))
                  limup = limdown + psusiz(ielt, itemp(itabt)) - 1

                  ihigh = limdown

                  DO ien = 1, ne

                     energy = (ear(ien-1) + ear(ien))*z15

                     IF ( energy .GE. penerg(limdown) .AND.
     &                    energy .LE. penerg(limup) ) THEN

                        DO WHILE ( energy .GT. penerg(ihigh) .AND. 
     &                       ihigh .LE. limup )
                           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) )
     &                       / (penerg(ihigh)-penerg(ilow))

                     ENDIF

                  ENDDO

               ENDIF

c End loop over elements.

            ENDDO

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

            DO il = linpt(itemp(itabt)), 
     &              linpt(itemp(itabt))+linsiz(itemp(itabt))-1

               energy = lineE(il) / z1

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( inputt(it) * keV_erg /
     &              atomM( lineEl(il) ) / amu_g + 
     &              velocity*velocity*km_cm*km_cm)
               ELSE
                  width = (energy / c_cm_s) * velocity * km_cm
               ENDIF

               lineflux = tcoeff(itabt) * dem(it) *
     &                    abunZ(lineEl(il)) * lineEm(il)

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

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./z1
         DO ien = 1, ne
            photar(ien) = photar(ien) * zf
         ENDDO
      ENDIF

      RETURN
      END

