      subroutine fstepr2(unit,hdunum,radin,radout,temp,pres,nrhdimj,
     $                idat1,rdat1,kdat1,nptrs,npnxt,npfi,npar,
     $                nplin,nlsvn,rcem,oplin,tau0,nloopctl,
     $                lun11,status)
C
C     Append a FITS extension binary table containing
C     nrhs columns and at most nrhdimj rows
c     author: T. Bridgman
C
C     Parameters:
C        unit    integer            File unit number
C        hdunum  integer            Number of last HDU written
C        radin   real               inner radius of shell
C        radout  real               outer radius of shell
C        temp    real               temperature of shell
C        pres    real               pressure in shell
C        nrhdimj  integer            Maximum number of rows
C        idat1   integer(nidat1)    Needed by the atomic database
C        rdat1   real(nidat1)       Needed by the atomic database
C        kdat1   char*nidat1        Needed by the atomic database
C        nptrs                      Needed by the atomic database
C        npnxt                      Needed by the atomic database
C        npfi                       Needed by the atomic database
C        npfirst                    Needed by the atomic database
C        npcon                      Needed by the atomic database
C        npconi                     Needed by the atomic database
C        npcon2                     Needed by the atomic database
C        xilev   real(nrhdimj)       Fractional level population array
C        cemab   real(2,nrhdimj)     Recombination emission
C        opakab  real(nrhdimj)       Opacity
C        tauc    real(2,nrhdimj)     Optical depth
C        poptol  real               Tolerance for population level
C        nloopctl integer           Loop control variable
C        nzone   integer            Pass number through iteration process
C        status  integer            Returned status code
C
      implicit none
      include './PARAM'
      integer mllz, nrhdimj

      integer nptmpdim
      parameter (nptmpdim=200000)
c
C     Allocation for passed parameters
      real tau0(2,nrhdimj), rcem(2,nrhdimj)
      real rdat1(nrdat1)
      real radin, radout, temp, pres
      integer unit,hdunum, nrows, status, nloopctl
      integer idat1(nidat1),nptrs(nptt,ndat2)
c     line opacities
      real oplin(nnnl)

c     pointers to master data
      integer npnxt(ndat2)
      integer npfi(ntyp,nni),npar(ndat2)
      integer nplin(nnnl)

C     Internal work areas
      real rwrk1(50000)
      integer ntptr(50000)
      character*10 kion(50000)
      character*20 klevl(50000),klevu(50000),kblnk20
      integer tfields,varidat
      character*16 ttype(10),tform(10),tunit(10)
      integer colnum,frow,felem,hdutype,ll, ltyp
      integer lrtyp, lcon, nrdt, nidt, mm, lun11, lpril,lpri
      integer jkk, nlev, mlpar
      integer nlplmx,ln,lnn,ml,nilin2,nlpl,lmm,kltmpn,kltmpo,
     $         llo,lup,llofnd,lupfnd,nlines,nlsvn,
     $         k,kl2,kk,nilin,lm,mlm
      real eliml,elimh,elin,elmmtpp,elcomp
      character*33 extname
      character*20 ktmp2
      character*1 kdat1(nkdat1)

C     Database manipulation quantities
      real rdat(nptmpdim)
      integer idat(nptmpdim), nkdt
      character*1 kdat(nptmpdim),kblnk,kdtmp(200)
      integer kltmp(50000)
      real elsv(50000)
      logical done

      data kblnk/' '/
      data kblnk20/'                    '/
c
      data tform/'1J','1E','8A','20A','20A','1E','1E','1E',
     $ '1E','1E'/

      data ttype/'index','wavelength','ion',
     $ 'lower_level','upper_level','emis_inward',
     $ 'emis_outward','opacity','tau_in','tau_out'/

      data tunit/' ','A',' ',' ',' ','erg/cm^3/s',
     $ 'erg/cm^3/s','/cm',' ',' '/

      varidat=0
c
      lpri=0
      lpril=lpri
c

C     Move to the last HDU (hdunum) in the file
      if (lpri.ne.0)
     $ write(lun11,*)'fstepr: Moving to end-of-FITS file'
      call ftmahd(unit,hdunum,hdutype,status)
      if (status .gt. 0)call printerror(lun11,status)

C     append a new empty extension after the last HDU
      if (lpri.ne.0)
     $ write (lun11,*)'fstepr: Create the new extension'
      call ftcrhd(unit,status)
      if (status .gt. 0)call printerror(lun11,status)

C----------------------------------------------------------------
C
C     Extracting data from the Atomic Database here
C
      if (lpri.ne.0)
     $ write (lun11,*)' '
c      lpril=0
      kltmpo=0
c
c     print important lines
      if (lpri.ne.0)
     $ write (lun11,*)'emission line luminosities (erg/sec/10**38))'
         nlplmx=50000
c
         eliml=0.1
         elimh=1.0e10
c
c        find the strongest lines.
         do  lm=1,nlplmx
            kltmp(lm)=0
            enddo
c
         nlpl=1
         do lnn=1,nlsvn
c
            ln=lnn
            ml=nplin(ln)
            mlm=ml-1
            call dread(ltyp,lrtyp,lcon,
     $        nrdt,rdat,nidt,idat,nkdt,kdat,mlm,
     $        idat1,rdat1,kdat1,nptrs,0,lun11)
            elin=abs(rdat(1))
c            if (lrtyp.eq.14)  elin=rdat(nrdt)
            if (lrtyp.ne.14) then
              nilin=npar(ml)
              mlm=nilin-1
              call dread(ltyp,lrtyp,lcon,
     $          nrdt,rdat,nidt,idat,nkdt,kdat,mlm,
     $          idat1,rdat1,kdat1,nptrs,0,lun11)
              nilin2=idat(nidt)
              elmmtpp=(rcem(2,ln)+rcem(1,ln))/2.
              if (lpril.ne.0)
     $         write (lun11,*)lnn,elin,nilin,elmmtpp,oplin(ln),ln,ml
              if ((ln.gt.0).and.(ln.le.nnnl).and.
     $           (elin.gt.eliml).and.(elin.lt.elimh).and.
     $           (elin.lt.8.9e+4).and.(elmmtpp.gt.1.e-36).and.
     $           (oplin(ln).gt.1.e-36).and.(nilin2.gt.0).and.
     $           (nilin2.le.nni)) then
c
                lmm=0
                elcomp=1.e+10
                do while ((lmm.lt.nlpl).and.(elmmtpp.lt.elcomp))
                  lmm=lmm+1
                  kl2=kltmp(lmm)
                  elcomp=0.
                  if (kl2.gt.0)
     $              elcomp=(rcem(2,kl2)+rcem(1,kl2))/2.
                  enddo
c
                if (lpril.ne.0)
     $            write (lun11,8516)ln,elin,elmmtpp
8516            format (1h ,i4,2e12.4)
                kltmpo=ln
                do  k=lmm,min(nlplmx,nlpl)
                  if (lpril.ne.0)
     $             write (lun11,*)'in 557 loop',k,kltmp(k),kltmpo
                  kltmpn=kltmp(k)
                  kltmp(k)=kltmpo
                  kltmpo=kltmpn
                  enddo
c
                nlpl=min(nlplmx,nlpl+1)
c
                if (lpril.ne.0)
     $           write (lun11,*)'done with 557 loop',lmm,nlpl
c
                endif
              endif
            enddo
         if (nlpl.gt.0) kltmp(nlpl)=kltmpo
c
c
c         nlpl=nlpl-1
c
         if (lpri.ne.0)
     $    write (lun11,959)
959      format (1x,'index, ion, wavelength, transmitted, reflected')
         do  kk=1,nlpl
           if (lpril.ne.0)
     $       write (lun11,*)'kk=',kk
           ln=kltmp(kk)
           if (ln.ne.0) then
             ml=nplin(ln)
             ntptr(kk)=ln
             klevl(kk)=kblnk20
             klevu(kk)=kblnk20
             if (ml.ne.0) then
               if (lpril.ne.0)
     $         write (lun11,*)'   ',ln,ml
               mlm=ml-1
               call dread(ltyp,lrtyp,lcon,
     $           nrdt,rdat,nidt,idat,nkdt,kdat,mlm,
     $           idat1,rdat1,kdat1,nptrs,0,lun11)
               llo=idat(1)
               lup=idat(2)
               elsv(kk)=abs(rdat(1))
               nilin=npar(ml)
               mlm=nilin-1
               call dread(ltyp,lrtyp,lcon,
     $           nrdt,rdat,nidt,idat,nkdt,kdat,mlm,
     $           idat1,rdat1,kdat1,nptrs,0,lun11)
               do mm=1,nkdt
                 kdtmp(mm)=kdat(mm)
                 enddo
               do mm=nkdt+1,9
                 kdtmp(mm)=kblnk
                 enddo
c               nilin=idat(3)
               write(kion(kk),'(10a1)')(kdtmp(mm),mm=1,10)
               done=.false.
               jkk=1
               do while (.not.done)
                 ml=npfi(13,jkk)
                 if (ml.ne.0) then
                   if ((npar(ml).eq.nilin).or.(jkk.gt.nni)) 
     $               done=.true.
                   endif
                 jkk=jkk+1
                 enddo
               if (jkk.gt.nni) ml=0
               if (ml.ne.0) then
                 mllz=npar(ml)
                 mlpar=npar(ml)
                 lupfnd=0
                 llofnd=0
                 do while ((ml.ne.0).and.(mlpar.eq.mllz)
     $             .and.((llofnd.ne.1).or.(lupfnd.ne.1)))
                   mlm=ml-1
                   call dread(ltyp,lrtyp,lcon,
     $               nrdt,rdat,nidt,idat,nkdt,kdat,mlm,
     $               idat1,rdat1,kdat1,nptrs,0,lun11)
                   nlev=idat(nidt-1)
                   if (lpri.ne.0)
     $               call dprinto(ltyp,lrtyp,lcon,
     $               nrdt,rdat,nidt,idat,nkdt,kdat,lun11)
                   if (lpri.ne.0)
     $               write (lun11,*)nlev,llo,lup,llofnd,lupfnd
                   if (nlev.eq.llo) then
                     do mm=1,20
                       if (mm.le.nkdt) then
                           write(ktmp2(mm:mm),'(a1)')kdat(mm)
                         else
                           write(ktmp2(mm:mm),'(a1)')kblnk
                         endif
                       enddo
c                     write (lun11,*)kk,ktmp2
                     klevl(kk)=ktmp2
                     llofnd=1
                     endif
                   if (nlev.eq.lup) then
                     do mm=1,20
                       if (mm.le.nkdt) then
                           write(ktmp2(mm:mm),'(a1)')kdat(mm)
                         else
                          write(ktmp2(mm:mm),'(a1)')kblnk
                         endif
                       enddo
                     klevu(kk)=ktmp2
                     lupfnd=1
                     endif
                   ml=npnxt(ml)
                   mlpar=0
                   if (ml.ne.0) mlpar=npar(ml)
                   enddo
                 if (lpril.ne.0) then
                   write (lun11,*)ml,nilin,npar(ml)
                   write (lun11,9955)kk,ln,(kdtmp(mm),mm=1,8),elsv(kk),
     $               rcem(1,ln),rcem(2,ln)
                   write (lun11,*)klevu(kk)
                   write (lun11,*)klevl(kk)
                   endif
 9955            format (1x,2i8,1x,8a1,3(1pe11.3))
                 endif
               endif
             endif
           enddo
c
c      if (nlpl.le.0) return
      nlpl=max(nlpl,1)
c

C     End of atomic database extraction
C----------------------------------------------------------------
C     define parameters for the binary table (see the above data statements)
      nrows=nlpl
      if (lpri.ne.0)
     $ write (lun11,*)'before header write'
      tfields=10
C     Build extension name
      extname='XSTAR_RADIAL'
      if(nloopctl.gt.0) then
          write(ktmp2,'(I4.4)')nloopctl
          extname='XSTAR_RADIAL_' // ktmp2
          endif

      if (lpri.ne.0)
     $ write (lun11,*)'fstepr: Write table headers'
C     write the required header parameters for the binary table
      call ftphbn(unit,nrows,tfields,ttype,tform,tunit,extname,
     $              varidat,status)
      if (status .gt. 0)call printerror(lun11,status)

      if (lpri.ne.0)
     $ write (lun11,*)'fstepr: Add some more keywords'

C     Write some model parameters in the extension header
      call ftpcom(unit,'***********************************',status)
      if (status .gt. 0)call printerror(lun11,status)

      call ftpcom(unit,'Model Keywords',status)
      if (status .gt. 0)call printerror(lun11,status)

C     Write values to 3 decimal places
      call ftpkye(unit,'RINNER',radin,3,'[cm] Inner shell radius',
     $ status)
      if (status .gt. 0)call printerror(lun11,status)

      call ftpkye(unit,'ROUTER',radout,3,'[cm] Outer shell radius',
     $ status)
      if (status .gt. 0)call printerror(lun11,status)

      call ftpkye(unit,'TEMPERAT',temp,3,'[10**4K] Shell Temperature',
     $ status)
      if (status .gt. 0)call printerror(lun11,status)

      call ftpkye(unit,'PRESSURE',pres,3,'[dynes/cm**2] Shell Pressure',
     $ status)
      if (status .gt. 0)call printerror(lun11,status)

      if (lpri.ne.0)
     $ write (lun11,*)'after header write'
C-------------------------------------------------------------------
C     Step through the columns and write them to the file
C
C     set 'global' parameters for writing FITS columns
      frow=1
      felem=1


C     column  1  (Line number)
      colnum=1
      if (lpri.ne.0)
     $ write(lun11,*)'fstepr: Writing Column ',colnum,nlpl
      nlines=nlpl
      call ftpclj(unit,colnum,frow,felem,nlines,ntptr,status)
      if (status .gt. 0)call printerror(lun11,status)


C     column  2  (wavelength)
      colnum=2
      if (lpri.ne.0)
     $ write(lun11,*)'fstepr: Writing Column ',colnum
      call ftpcle(unit,colnum,frow,felem,nlines,elsv,status)
      if (status .gt. 0)call printerror(lun11,status)
      if (status .gt. 0) stop


C     column  3  (Ion)
      colnum=3
      if (lpri.ne.0)
     $ write(lun11,*)'fstepr: Writing Column ',colnum
      call ftpcls(unit,colnum,frow,felem,nlines,kion,status)
      if (status .gt. 0)call printerror(lun11,status)


C     column  4 (lower Level Designation)
      colnum=4
      if (lpri.ne.0)
     $ write (lun11,*)'fstepr: Writing Column ',colnum
      call ftpcls(unit,colnum,frow,felem,nlines,klevl,status)
      if (status .gt. 0)call printerror(lun11,status)

C     column  5 (Level Designation)
      colnum=5
      if (lpri.ne.0)
     $ write(lun11,*)'fstepr: Writing Column ',colnum
      call ftpcls(unit,colnum,frow,felem,nlines,klevu,status)
      if (status .gt. 0)call printerror(lun11,status)

C----------------------------------------------------------------

C     column  6
      colnum=6
      do ll=1,nlines
         rwrk1(ll)=0.
         if (ntptr(ll).ne.0)
     $      rwrk1(ll)=rcem(1,ntptr(ll))
         enddo
      if (lpri.ne.0)
     $ write(lun11,*)'fstepr: Writing Column ',colnum
      call ftpcle(unit,colnum,frow,felem,nlines,rwrk1,status)
      if (status .gt. 0)call printerror(lun11,status)

C     column  7
      colnum=7
      do ll=1,nlines
         rwrk1(ll)=0.
         if (ntptr(ll).ne.0)
     $    rwrk1(ll)=rcem(2,ntptr(ll))
         enddo
      if (lpri.ne.0)
     $ write (lun11,*)'fstepr: Writing Column ',colnum
      call ftpcle(unit,colnum,frow,felem,nlines,rwrk1,status)
      if (status .gt. 0)call printerror(lun11,status)


C     column  8
      colnum=8
      do ll=1,nlines
         rwrk1(ll)=0.
         if (ntptr(ll).ne.0)
     $    rwrk1(ll)=oplin(ntptr(ll))
         enddo
      if (lpri.ne.0)
     $ write (lun11,*)'fstepr: Writing Column ',colnum
      call ftpcle(unit,colnum,frow,felem,nlines,rwrk1,status)
      if (status .gt. 0)call printerror(lun11,status)


C     column  9
      colnum=9
      do ll=1,nlines
         rwrk1(ll)=0.
         if (ntptr(ll).ne.0)
     $    rwrk1(ll)=tau0(1,ntptr(ll))
         enddo
      if (lpri.ne.0)
     $ write (lun11,*)'fstepr: Writing Column ',colnum
      call ftpcle(unit,colnum,frow,felem,nlines,rwrk1,status)
      if (status .gt. 0)call printerror(lun11,status)

C     column  10
      colnum=10
      do ll=1,nlines
         rwrk1(ll)=0.
         if (ntptr(ll).ne.0)
     $    rwrk1(ll)=tau0(2,ntptr(ll))
         enddo
      if (lpri.ne.0)
     $ write (lun11,*)'fstepr: Writing Column ',colnum
      call ftpcle(unit,colnum,frow,felem,nlines,rwrk1,status)
      if (status .gt. 0)call printerror(lun11,status)


c----------------------------------------------------------------
C     Compute checksums
      call ftpcks(unit,status)
      if (status .gt. 0)call printerror(lun11,status)

      return
      end
