
      SUBROUTINE kbmodl

c  Runs Kazik Borkowski's NEI/SNR MEKAL code to produce tabulated data 
c  files for use by XSPEC.
c
C    v1.00  5/26/99   kaa
C    v1.10 10/24/00   kaa  Updated for the new models which include a lower
C                          ionization age parameter. This is set to zero for
C                          the grids. Also changed the names of the models
C                          to pshock and npshock to match those now used in
C                          XSPEC.
c
 

C Pointers for dynamic arrays

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


      REAL tstart, tend, taust, tauend, testrt, teend

      INTEGER nset, ntsteps, ntaust, nenrgy, status, nret, ntestp

      CHARACTER filename*72, efilenm*72, contxt*72, kbmdnm*6

      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

      INTEGER ldatdir
      CHARACTER datdir*128

      CHARACTER trfnam*128
      EXTERNAL trfnam

      CHARACTER*40 taskname
      COMMON /task/ taskname

      CALL fcecho('KBMODELS vers 1.10    24 Feb 2000.')

      taskname = 'kbmodels v1.10'

C Get the parameters

      CALL gpkbmd(kbmdnm, filename, tstart, tend, ntsteps, 
     &            efilenm, taust, tauend, ntaust, testrt, teend, 
     &            ntestp, qclob, status)
      contxt = 'Error getting parameters'
      IF ( status .NE. 0 ) GOTO 999

C Load the array that specifies where the data files are to be
C found. This is either the standard XSPEC data directory or a
c special directory defined by the environmental variable 
c XSPEC_MDATA_DIR.

      CALL trlog('XSPEC_MDATA_DIR', 15, datdir, ldatdir)
      IF ( ldatdir .EQ. 0 ) THEN
         datdir = trfnam('$XANADU,spectral/xspec/manager,')
      ENDIF
      CALL fpdatd(datdir)

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 KBMODELS model for ', 
     &                         nenrgy, ' energy bins'
      CALL fcecho(contxt)

C Get the memory that will be required

      CALL udmget(nenrgy+1, 6, iefilear, status)
      contxt = 'Insufficient memory for efilear 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(3*MAX(ntsteps,ntaust,ntestp), 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, MEMR(iestep), 
     &            MEMR(iemin), MEMI(insteps), MEMR(iefilear))

C Call the routine that actually does the work

      CALL dokbmd(kbmdnm, filename, nset, ntsteps, ntaust, 
     &            MAX(ntsteps, ntaust, ntestp), nenrgy, tstart, 
     &            tend, taust, tauend, testrt, teend, ntestp,
     &            MEMR(iestep), MEMR(iemin), MEMI(insteps), 
     &            MEMR(iefilear), 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 gpkbmd(kbmdnm, filename, tstart, tend, 
     &                  ntsteps, efilenm, taust, tauend, ntaust, 
     &                  testrt, teend, ntestp, qclob, status)

      REAL tstart, tend, taust, tauend, testrt, teend

      INTEGER ntsteps, ntaust, ntestp, status

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

      LOGICAL qclob

C Gets the parameters from the raysmith parameter file
C      filename       c       r: Output filename
C      tstart         r       r: Low shock temperature
C      tend           r       r: High shock temperature
C      ntsteps        i       r: The number of (logarithmic) temperature steps
C      efilenm        c       r: Filename specifying energy bins
C      taust          r       r: Low tau
C      tauend         r       r: High tau
C      ntaust         i       r: The number of (logarithmic) tau steps
C      testrt         r       r: Low electron temperature
C      teend          r       r: High electron temperature
C      ntestp         i       r: The number of electron temperatures
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 model type required

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

      IF ( kbmdnm .NE. 'pshock' .AND. kbmdnm .NE. 'npshock' .AND. 
     &     kbmdnm .NE. 'sedov' ) THEN
         CALL fcerr('Modelname should be pshock, npshock, or sedov')
         CALL fcerrm(0)
      ENDIF

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 lower ionization timescale

      CALL uclgsr('taustart', taust, status)
      contxt = 'Cannot get the taustart parameter'
      IF ( status .NE. 0 ) GOTO 999
      
C Get the upper ionization timescale

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

      CALL uclgsi('ntaust', ntaust, status)
      contxt = 'Cannot get the ntaust parameter'
      IF ( status .NE. 0 ) GOTO 999
      
      IF ( kbmdnm .EQ. 'sedov' .OR. kbmdnm .EQ. 'nshock' ) THEN

C Get the lower electron temperature

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

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

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

      ELSE

         ntestp = 0

      ENDIF
      
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, estep, emin, 
     &                  nsteps, efilear)

      INTEGER nset, nenrgy

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

      INTEGER nsteps(nset)

      CHARACTER*(*) efilenm

C Subroutine to read the EFILENM file and set up the energy bin arrays

      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
      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)

         ENDDO

      ENDDO

      CLOSE(ilun)
      CALL frelun(ilun)

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

      SUBROUTINE dokbmd(kbmdnm, filename, nset, ntsteps, ntaust, 
     &                  maxtab, nenrgy, tstart, tend, taust, tauend, 
     &                  testrt, teend, ntestp, estep, emin, nsteps, 
     &                  efilear, parval, qclob, work1, work2)

      INTEGER nset, ntsteps, ntaust, ntestp, nenrgy, maxtab

      REAL efilear(0:nenrgy), parval(maxtab, 3)
      REAL estep(nset) , emin(nset)
      REAL work1(nenrgy) , work2(nenrgy)

      INTEGER nsteps(nset)

      REAL tstart, tend, taust, tauend, testrt, teend

      CHARACTER filename*64, kbmdnm*(*)

      LOGICAL qclob

C Subroutine to calculate the KBMODELS models for temperatures from
C tstart to tend in ntsteps logarithmic steps and ionization timescales
c taust to tauend in ntaust logarithmic steps.


      INTEGER NELEMENT
      PARAMETER (NELEMENT=11)

      REAL addivl(6, NELEMENT), intivl(6,3)
      REAL tstep, t, taustep, testep, param(17)

      INTEGER intntb(3), intmth(3)
      INTEGER ierr, iunit, nmintp
      INTEGER i, j, k, l, irec, status
 
      CHARACTER contxt*72, mdlnam*9
      CHARACTER*8 elements(NELEMENT), intnam(3)

      DATA elements/'He' , 'C ' , 'N ' , 'O ' , 'Ne', 'Mg' ,
     &              'Si' , 'S ' , 'Ca' , 'Fe' , 'Ni'/
      DATA intnam/'T_s', 'tau', 'T_e'/
      DATA intmth/1,1,1/

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) = 2.
      intivl(2,1) = 0.01
      intivl(3,1) = tstart
      intivl(4,1) = tstart
      intivl(5,1) = tend
      intivl(6,1) = tend

      intivl(1,2) = 1.e11
      intivl(2,2) = 1.e8
      intivl(3,2) = taust
      intivl(4,2) = taust
      intivl(5,2) = tauend
      intivl(6,2) = tauend

      IF ( ntestp .GT. 0 ) THEN
         intivl(1,3) = 2.
         intivl(2,3) = 0.01
         intivl(3,3) = testrt
         intivl(4,3) = testrt
         intivl(5,3) = teend
         intivl(6,3) = teend
      ENDIF

      intntb(1) = ntsteps
      intntb(2) = ntaust
      intntb(3) = ntestp

      IF ( ntestp .GT. 0 ) THEN
         nmintp = 3
      ELSE
         nmintp = 2
      ENDIF

c Store the temperatures and ionization timescales to be used in PARVAL.

      IF ( ntsteps .GT. 1 ) THEN
         tstep = (ALOG10(tend)-ALOG10(tstart))/(ntsteps-1)
      ELSE
         tstep = 0.
      ENDIF

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

      ENDDO

      IF ( ntaust .GT. 1 ) THEN
         taustep = (ALOG10(tauend)-ALOG10(taust))/(ntaust-1)
      ELSE
         taustep = 0.
      ENDIF

      DO i = 1, ntaust
 
         t = ALOG10(taust) + (i-1)*taustep
         parval(i,2) = 10**t

      ENDDO

      IF ( ntestp .GT. 1 ) THEN
         testep = (ALOG10(teend)-ALOG10(testrt))/(ntestp-1)
      ELSE
         testep = 0.
      ENDIF

      DO i = 1, ntestp
 
         t = ALOG10(testrt) + (i-1)*testep
         parval(i,3) = 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.

      mdlnam = 'tab_'//kbmdnm(1:5)

      CALL wftbmd(filename, ' ', mdlnam, 'ph/cm^2/s', nmintp,
     &            NELEMENT, .TRUE., .TRUE., elements, addivl,
     &            intnam, intivl, intntb, intmth, maxtab, 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 Now do the calculation and write the spectra

      IF ( kbmdnm .EQ. 'pshock' ) THEN

         CALL do2dgr(kbmdnm, iunit, ntsteps, ntaust,
     &               maxtab, parval, nenrgy, efilear, nelement, 
     &               param, work1, work2)


      ELSEIF ( kbmdnm .EQ. 'sedov' .OR. kbmdnm .EQ. 'npshock' ) THEN

         CALL do3dgr(kbmdnm, iunit, ntsteps, ntaust, ntestp, 
     &               maxtab, parval, nenrgy, efilear, nelement, 
     &               param, work1, work2)

      ENDIF
 
      CALL ftclos(iunit, status)
      CALL FRELUN(iunit)

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

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

      SUBROUTINE do3dgr(kbmdnm, iunit, ntsteps, ntaust, ntestp,
     &                  maxtab, parval, nenrgy, efilear, nelement, 
     &                  param, work1, work2)

      INTEGER iunit, ntsteps, ntaust, ntestp, maxtab, nenrgy, nelement

      REAL parval(maxtab, 3), efilear(0:nenrgy), param(17)
      REAL work1(nenrgy), work2(nenrgy)

      CHARACTER kbmdnm*(*)


      INTEGER i, j, k, l, m, irec, status

      CHARACTER contxt*80

c Set the initial element abundances to a small number except for H
 
      param(3) = 1.
      DO i = 4 , nelement+3
         param(i) = 1E-7
      ENDDO
      param(15) = 0.
      param(17) = 0.

c Loop round the required temperatures using log steps.

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

         param(1) = parval(j,1)

c Loop round the required ionization timescale using log steps.

         DO k = 1, ntaust

            WRITE (*,'(a,1pe12.3,a)') 
     &                           '                       and tau= ',
     &                              parval(k,2)

            param(16) = parval(k,2)

c Loop round the required temperature ratios

            DO m = 1, ntestp

               WRITE (*,'(a,1pe12.3,a)') 
     &                           '                   and T_e/keV= ',
     &                              parval(m,3)

               param(2) = parval(m,3)

               irec = irec + 1

c Write the parameter value column

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

c Calculate zero-metal continuum

               IF ( kbmdnm .EQ. 'sedov' ) THEN
                  CALL xssedov(efilear, nenrgy, param, 1, work1)
               ELSEIF ( kbmdnm .EQ. 'npshock' ) THEN
                  CALL xsshock(efilear, nenrgy, param, 1, work1)
               ENDIF

c Write the zero-metal continuum

               CALL ftpcle(iunit, 2, irec, 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

                  param(i+3) = 1.
                  IF ( kbmdnm .EQ. 'sedov' ) THEN
                     CALL xssedov(efilear, nenrgy, param, 1, work2)
                  ELSEIF ( kbmdnm .EQ. 'npshock' ) THEN
                     CALL xsshock(efilear, nenrgy, param, 1, work2)
                  ENDIF
                  param(i+3) = 1.E-7

                  DO l = 1, nenrgy
                     work2(l) = work2(l) - work1(l)
                  ENDDO

                  CALL ftpcle(iunit, i+2, irec, 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 temperature ratios

            ENDDO

c End loop over ionization timescales

         ENDDO
 
c End loop over temperatures
 
      ENDDO

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

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

      SUBROUTINE do2dgr(kbmdnm, iunit, ntsteps, ntaust,
     &                  maxtab, parval, nenrgy, efilear, nelement, 
     &                  param, work1, work2)

      INTEGER iunit, ntsteps, ntaust, maxtab, nenrgy, nelement

      REAL parval(maxtab, 3), efilear(0:nenrgy), param(16)
      REAL work1(nenrgy), work2(nenrgy)

      CHARACTER kbmdnm*(*)


      INTEGER i, j, k, l, irec, status

      CHARACTER contxt*80


c Set the initial element abundances to a small number except for H
 
      param(2) = 1.
      DO i = 3 , nelement+2
         param(i) = 1E-7
      ENDDO
      param(14) = 0.
      param(16) = 0.

c Loop round the required temperatures using log steps.

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

         param(1) = parval(j,1)

c Loop round the required ionization timescale using log steps.

         DO k = 1, ntaust

            WRITE (*,'(a,1pe12.3,a)') 
     &                           '                         and tau= ',
     &                              parval(k,2)

            param(15) = parval(k,2)

            irec = irec + 1

c Write the parameter value column

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

c Calculate zero-metal continuum

            IF ( kbmdnm .EQ. 'pshock' ) THEN
               CALL xsneqs(efilear, nenrgy, param, 1, work1)
            ENDIF

c Write the zero-metal continuum

            CALL ftpcle(iunit, 2, irec, 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

               param(i+2) = 1.
               IF ( kbmdnm .EQ. 'pshock' ) THEN
                  CALL xsneqs(efilear, nenrgy, param, 1, work2)
               ENDIF
               param(i+2) = 1.E-7

               DO l = 1, nenrgy
                  work2(l) = work2(l) - work1(l)
               ENDDO

               CALL ftpcle(iunit, i+2, irec, 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 ionization timescales

         ENDDO
 
c End loop over temperatures
 
      ENDDO

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

