
      SUBROUTINE mekal

c  Runs the MEKAL code to produce tabulated data files for use by 
C  XSPEC.
c
C    v1.00  1/10/96   kaa
c     1.01  9/22/01   kaa   Fixed uninitialized variables
c
 

C Pointers for dynamic arrays

      INTEGER iefilear, iestep, iemin, insteps, iparval
      INTEGER ienorm, iwork1, iwork2


      REAL tstart, tend, density

      INTEGER nset, ntsteps, nenrgy, status, nret

      CHARACTER filename*72, efilenm*72, contxt*72

      LOGICAL qclob
 

* Dynamic memory allocation stuff
C  the following MEM common block definition is in the system iraf77.inc file
C
C Get IRAF MEM common into main program.
C
      LOGICAL          MEMB(1)
      INTEGER*2        MEMS(1)
      INTEGER*4        MEMI(1)
      INTEGER*4        MEML(1)
      REAL             MEMR(1)
      DOUBLE PRECISION MEMD(1)
      COMPLEX          MEMX(1)
      EQUIVALENCE (MEMB, MEMS, MEMI, MEML, MEMR, MEMD, MEMX)
      COMMON /MEM/ MEMD

C     datatype gives a symbolic code for the data type, e.g.,  4 = Integer*4
C     1 is boolean
C     3 is short integer
C     4 is integer
C     5 is long integer
C     6 is single precision real
C     7 is double precision real
C     8 complex

      CHARACTER*40 taskname
      COMMON /task/ taskname

      status = 0
 
      CALL fcecho('MEKAL vers 1.01    22 Sep 2001.')

      taskname = 'mekal v1.01'

C This ftool needs the XANADU variable set in order to find the input
C files so check for it

      CALL trlog('XANADU', 6, contxt, nret)
      IF ( nret .EQ. 0 ) THEN
         CALL fcecho('Please set up the XANADU environment variable')
         RETURN
      ENDIF

C Get the parameters

      CALL gpmek(filename, tstart, tend, ntsteps, efilenm, density, 
     &           qclob, status)
      contxt = 'Error getting parameters'
      IF ( status .NE. 0 ) GOTO 999

C Get the number of energy bins required

      CALL gtnebn(efilenm, nset, nenrgy, status)
      contxt = 'Error getting number of energy bins required'
      IF ( status .NE. 0 ) GOTO 999

      WRITE(contxt,'(a,i6,a)') 'Calculating MEKAL model for ', nenrgy, 
     &                         ' energy bins'
      CALL fcecho(contxt)

C Get the memory that will be required

      iefilear = 0
      ienorm = 0
      iwork1 = 0
      iwork2 = 0
      iestep = 0
      iemin = 0
      insteps = 0
      iparval = 0

      CALL udmget(nenrgy+1, 6, iefilear, status)
      contxt = 'Insufficient memory for efilear array'
      IF ( status .NE. 0 ) GOTO 999

      CALL udmget(nenrgy, 6, ienorm, status)
      contxt = 'Insufficient memory for enorm array'
      IF ( status .NE. 0 ) GOTO 999

      CALL udmget(nenrgy, 6, iwork1, status)
      contxt = 'Insufficient memory for work1 array'
      IF ( status .NE. 0 ) GOTO 999

      CALL udmget(nenrgy, 6, iwork2, status)
      contxt = 'Insufficient memory for work2 array'
      IF ( status .NE. 0 ) GOTO 999

      CALL udmget(nset, 6, iestep, status)
      contxt = 'Insufficient memory for estep array'
      IF ( status .NE. 0 ) GOTO 999

      CALL udmget(nset, 6, iemin, status)
      contxt = 'Insufficient memory for emin array'
      IF ( status .NE. 0 ) GOTO 999

      CALL udmget(nset, 4, insteps, status)
      contxt = 'Insufficient memory for nsteps array'
      IF ( status .NE. 0 ) GOTO 999

      CALL udmget(ntsteps, 6, iparval, status)
      contxt = 'Insufficient memory for parval array'
      IF ( status .NE. 0 ) GOTO 999

C Read the file with the energy bin information and set up
C the arrays

      CALL stebin(efilenm, nset, nenrgy, density, MEMR(iestep), 
     &            MEMR(iemin), MEMI(insteps), MEMR(iefilear), 
     &            MEMR(ienorm))

C Call the routine that actually does the work

      CALL domekl(filename, nset, ntsteps, nenrgy, tstart, tend, 
     &            density, MEMR(iestep), MEMR(iemin), MEMI(insteps), 
     &            MEMR(iefilear), MEMR(ienorm), MEMR(iparval), 
     &            qclob, MEMR(iwork1), MEMR(iwork2))


 999  CONTINUE
      IF ( status .NE. 0 ) THEN
         CALL fcerr(contxt)
         CALL fcerrm(status)
      ENDIF
 
      RETURN
      END

C *********************************************************************

      SUBROUTINE gpmek(filename, tstart, tend, ntsteps, efilenm, 
     &                 density, qclob, status)

      REAL tstart, tend, density

      INTEGER ntsteps, status

      CHARACTER filename*(*), efilenm*(*)

      LOGICAL qclob

C Gets the parameters from the raysmith parameter file
C      filename       c       r: Output filename
C      tstart         r       r: Low temperature
C      tend           r       r: High temperature
C      ntsteps        i       r: The number of (logarithmic) temperature steps
C      efilenm        c       r: Filename specifying energy bins
C      density        r       r: Density (cm^-3)
C      qclob          l       r: Whether to overwrite existing output file
C      status         i       r: i/o status

      CHARACTER contxt*72

      status = 0

C Get the output filename

      CALL uclgst('filename', filename, status)
      contxt = 'Cannot get the filename parameter'
      IF ( status .NE. 0 ) GOTO 999

C Get the lower temperature

      CALL uclgsr('tstart', tstart, status)
      contxt = 'Cannot get the tstart parameter'
      IF ( status .NE. 0 ) GOTO 999
      
C Get the upper temperature

      CALL uclgsr('tend', tend, status)
      contxt = 'Cannot get the tend parameter'
      IF ( status .NE. 0 ) GOTO 999
      
C Get the number of temperatures

      CALL uclgsi('ntsteps', ntsteps, status)
      contxt = 'Cannot get the ntsteps parameter'
      IF ( status .NE. 0 ) GOTO 999
      
C Get the filename with energy bins

      CALL uclgst('efilenm', efilenm, status)
      contxt = 'Cannot get the efilnm parameter'
      IF ( status .NE. 0 ) GOTO 999

C Get the density

      CALL uclgsr('density', density, status)
      contxt = 'Cannot get the density parameter'
      IF ( status .NE. 0 ) GOTO 999
      
C Get the clobber parameter

      CALL uclgsb('clobber', qclob, status)
      contxt = 'Cannot get the clobber parameter'
      IF ( status .NE. 0 ) GOTO 999

 999  CONTINUE
      IF ( status .NE. 0 ) THEN
         CALL fcerr(contxt)
         CALL fcerrm(status)
      ENDIF

      RETURN
      END

C *********************************************************************

      SUBROUTINE gtnebn(efilenm, nset, nenrgy, status)

      INTEGER nset, nenrgy, status

      CHARACTER efilenm*(*)

C Subroutine to read the EFILENM file and find out how many energy bins
C and sets of energy bins are specified.

      REAL edelta

      INTEGER ilun, enumber

      CHARACTER contxt*72

      LOGICAL qdone

      INTEGER lenact
      EXTERNAL lenact

      status = 0

      CALL getlun(ilun)
      OPEN(unit=ilun, file=efilenm, status='old', iostat=status)
      contxt = 'Failed to open '//efilenm(:lenact(efilenm))
      IF ( status .NE. 0 ) GOTO 999

      READ(ilun, *, iostat=status) edelta
      contxt = 'Failed to read first line of '
     &            //efilenm(:lenact(efilenm))
      IF ( status .NE. 0 ) GOTO 999

      READ(ilun, *, iostat=status) edelta, enumber
      contxt = 'Failed to read second line of '
     &            //efilenm(:lenact(efilenm))
      IF ( status .NE. 0 ) GOTO 999

      nset = 0
      nenrgy = 0
      qdone = .FALSE.

      DO WHILE ( .NOT.qdone )

         nset = nset + 1
         nenrgy = nenrgy + enumber

         READ(ilun, *, iostat=status) edelta, enumber
         IF ( status .NE. 0 ) qdone = .TRUE.

      ENDDO
      status = 0

 999  CONTINUE
      IF ( status .NE. 0 ) THEN
         CALL fcerr(contxt)
         CALL fcerrm(status)
      ENDIF

      CLOSE(ilun)
      CALL frelun(ilun)

      RETURN
      END

C *********************************************************************

      SUBROUTINE stebin(efilenm, nset, nenrgy, density, estep, emin, 
     &                  nsteps, efilear, enorm)

      INTEGER nset, nenrgy

      REAL efilear(0:nenrgy), enorm(nenrgy)
      REAL estep(nset), emin(nset)
      REAL density

      INTEGER nsteps(nset)

      CHARACTER*(*) efilenm

C Subroutine to read the EFILENM file and set up the energy bin arrays
C Also calculates enorm which converts the 1e-4 ph/m**2/s/keV output from
C fmekal to ph/cm**2/s/bin for XSPEC. also correct for the effect that 
C input normalisation is n_en_H V / d**2 (output of fline is 
C ~ n_H n_H V / d**2). 119.64 factor is to give meka model is the same 
C normalisation as the R-S model (factor = 4 * pi * (3.096)**2)

      REAL dfact

      INTEGER iset, ien, i, ilun, status

      INTEGER lenact
      EXTERNAL lenact

      status = 0

      CALL getlun(ilun)
      OPEN(unit=ilun, file=efilenm, status='old', iostat=status)

      READ(ilun, *) efilear(0)
      efilear(0) = efilear(0)/1000.


      ien = 0
      dfact = 119.64 / density**2
      DO iset = 1, nset

         READ(ilun, *) estep(iset), nsteps(iset)
         estep(iset) = estep(iset)/1000.

         emin(iset) = efilear(ien)
         DO i = 1, nsteps(iset)

            ien = ien + 1
            efilear(ien) = efilear(ien-1) + estep(iset)

            enorm(ien) = dfact * (efilear(ien)-efilear(ien-1))
 
         ENDDO

      ENDDO

      CLOSE(ilun)
      CALL frelun(ilun)

      RETURN
      END
 
C *********************************************************************

      SUBROUTINE domekl(filename, nset, ntsteps, nenrgy, tstart, tend, 
     &                  density, estep, emin, nsteps, efilear, enorm, 
     &                  parval, qclob, work1, work2)

      INTEGER nset, ntsteps, nenrgy

      REAL efilear(0:nenrgy), parval(ntsteps)
      REAL estep(nset) , emin(nset)
      REAL work1(nenrgy) , work2(nenrgy) , enorm(nenrgy)
      REAL density

      INTEGER nsteps(nset)

      REAL tstart, tend

      CHARACTER filename*64

      LOGICAL qclob

C Subroutine to calculate the MEKAL models for temperatures from
C tstart to tend in ntsteps logarithmic steps.


      INTEGER NELEMENT
      PARAMETER (NELEMENT=14)

      REAL frac(NELEMENT+1)
      REAL addivl(6, NELEMENT), intivl(6)
      REAL tstep, t, ed
 
      INTEGER ierr, iunit
      INTEGER i, j, l, status
 
      CHARACTER contxt*72
      CHARACTER*8 elements(NELEMENT)

      DATA elements/'He' , 'C ' , 'N ' , 'O ' , 'Ne' , 'Na', 'Mg' ,
     &              'Al' , 'Si' , 'S ' , 'Ar' , 'Ca' , 'Fe' , 'Ni'/

      status = 0

c Set the initial parameter values

      DO i = 1, NELEMENT
         addivl(1,i) = 1.
         addivl(2,i) = -0.01
         addivl(3,i) = 0.
         addivl(4,i) = 0.
         addivl(5,i) = 1000.
         addivl(6,i) = 1000.
      ENDDO
      intivl(1) = 1.
      intivl(2) = 0.01
      intivl(3) = 0.008
      intivl(4) = 0.008
      intivl(5) = 64.0
      intivl(6) = 64.0

c Store the temperatures to be used in PARVAL.

      tstep = (ALOG10(tend)-ALOG10(tstart))/(ntsteps-1)

      DO i = 1, ntsteps
 
         t = ALOG10(tstart) + (i-1)*tstep
         parval(i) = 10**t

      ENDDO

C If clobber is set then remove any existing file

      IF ( qclob ) THEN
         CALL clobber(filename, status)
         contxt = 'Failed to delete existing FITS file'
         IF ( status .NE. 0 ) GOTO 999
      ENDIF

c Write the header information, parameter extension, and energy
c extension of the table file.

      CALL wftbmd(filename, ' ', 'tab_mekal', 'emis meas', 1,
     &            NELEMENT, .TRUE., .TRUE., elements, addivl,
     &            'kT (keV)', intivl, ntsteps, 1, ntsteps, parval,
     &            nenrgy+1, efilear(0), iunit, ierr, status, contxt)
      WRITE(contxt,'(a,i2)') 'Failure in WFTBMD - code ', ierr
      IF ( status .NE. 0 ) GOTO 999

c Set the initial element abundances to a small number except for H
 
      frac(1) = 1.
      DO i = 2 , NELEMENT+1
         frac(i) = 1E-7
      ENDDO
 
c Loop round the required temperatures using log steps.

      DO j = 1 , ntsteps
 
         WRITE (*,'(a,f10.5,a)') '+** calculating model with T/keV= ' , 
     &                           parval(j)

c Write the parameter value column

         CALL ftpcle(iunit, 1, j, 1, 1, parval(j), status)
         contxt = 'Failed to write PARAMVAL column'
         IF ( status .NE. 0 ) GOTO 999

c Calculate zero-metal continuum

         CALL fmekal(efilear(0), efilear(1), work1, nenrgy, 1., 
     &               density, parval(j), frac, ed, status)
         contxt = 'Failure in MEKAL code'
         IF ( status .NE. 0 ) GOTO 999

         DO l = 1, nenrgy
            work1(l) = work1(l)*enorm(l) / ed
         ENDDO
 
c Write the zero-metal continuum

         CALL ftpcle(iunit, 2, j, 1, nenrgy, work1, status)
         contxt = 'Failed to write INTPSPEC column'
         IF ( status .NE. 0 ) GOTO 999
 
c Create the heavy metal spectra
 
         DO i = 1 , NELEMENT

            frac(i+1) = 1.
            CALL fmekal(efilear(0), efilear(1), work2, nenrgy, 1., 
     &                  density, parval(j), frac, ed, status)
            frac(i+1) = 1.E-7
            DO l = 1, nenrgy
               work2(l) = (enorm(l)*work2(l)/ed) - work1(l)
            ENDDO

            CALL ftpcle(iunit, i+2, j, 1, nenrgy, work2, status)
            WRITE(contxt,'(a,i3)') 'Failed to write column ', i+2
            IF ( status .NE. 0 ) GOTO 999

         ENDDO
 
c End loop over temperatures
 
      ENDDO
 
      CALL ftclos(iunit, status)
      CALL FRELUN(iunit)

 999  CONTINUE
      IF ( status .NE. 0 ) THEN
         CALL fcerr(contxt)
         CALL fcerrm(status)
      ENDIF
 
      RETURN
      END

