      subroutine fstepr(unit,hdunum,radin,radout,temp,pres,nrhdimj,
     $                idat1,rdat1,kdat1,nptrs,npnxt,npfi,
     $                npfirst,npar,npilev,
     $                xilev,bilev,rniss,nloopctl,
     $                lun11,status)
C
C
C     Write Data for each radial zone to an individual extension
c     author: T. Bridgman
C
C     Append a FITS extension binary table containing
C     nrhs columns and at most nrhdimj rows
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 nrhdimj
c
      integer nptmpdim
      parameter (nptmpdim=200000)
c
C     Allocation for passed parameters
      real xilev(nrhdimj),bilev(nrhdimj),rniss(nrhdimj)
      real rdat1(nrdat1)
      real radin, radout, temp, pres
      integer unit,hdunum, nrows, status, nloopctl
      integer idat1(nidat1),nptrs(nptt,ndat2)

c     pointers to master data
      integer npnxt(ndat2),npfirst(ntyp)
      integer npfi(ntyp,nni),npar(ndat2)
      integer npilev(nd,nni)
C     Internal work areas
      real rwrk1(nnml),rwrk2(nnml), elev(nnml)
      integer ntptr(nnml)
      integer natomic(nnml), mllev(nnml),nupper(nnml)
      character*10 kion(nnml)
      character*20 klevt(nnml)
      integer tfields,varidat
      character*16 ttype(9),tform(9),tunit(9)
      integer colnum,frow,felem,hdutype, klel, mlel, jk, ltyp
      integer lrtyp, lcon, nrdt, nidt, mmlv, mm, lun11, lpril,lpri
      integer mllel, klion, mlion, jkk, kl
      integer mt2, mlleltp, nnz, nions
      character*43 extname
      character*30 ktmp2
      character*1 kdat1(nkdat1)
C     Database manipulation quantities
      real rdat(nptmpdim), xeltp
      integer idat(nptmpdim), nkdt
c      integer lpridread
      character*1 kdat(nptmpdim),kdati(nptmpdim),kblnk
      character*1 klev(100,nd)
      real rnissl(nnml)
      real t,xee,xpx
      integer j,nkdti
      real rlev(10,nd)
      integer ilev(10,nd),nlpt(nd),iltp(nd),nlev
      integer mm2,mmtmp,kkkl,lk,mlm
      real eth

      data kblnk/' '/
      data tform/'1J','1I','1E','8A','1I','20A','1E','1E','1I'/
      data ttype/'index','ion_index','e_excitation','ion',
     $ 'atomic_number','ion_level','population','lte',
     $  'upper index'/
      data tunit/' ',' ','eV',' ',' ',' ',' ',' ',' '/
c
      varidat=0
      lpril=0
      lpri=lpril
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
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     Extracting data from the Atomic Database here
c
C
C     lpril is flag for printing debug information
       nions=0
      if (lpril.ne.0) then
        write (lun11,*)'raw data'
        do j=1,nnml
          if (xilev(j).gt.1.e-37)
     $     write (lun11,*)j,xilev(j)
          enddo
        endif
C      initialize line counter
       mmlv=0
C      First look for element data (jk is element index)
        klel=11
        mlel=npfirst(klel)
        jkk=0
        jk=0
        do while (mlel.ne.0)
          jk=jk+1
          mt2=mlel-1
          call dread(ltyp,lrtyp,lcon,
     $      nrdt,rdat,nidt,idat,nkdt,kdat,mt2,
     $      idat1,rdat1,kdat1,nptrs,0,lun11)
          mllel=idat(nidt)
          xeltp=rdat(1)
          nnz=idat(1)
          if (lpril.ne.0)
     $      write (lun11,*)'element:',jk,mlel,mllel,nnz,
     $                  (kdat(mm),mm=1,nkdt),xeltp
C         ignore if the abundance is small
          if (xeltp.lt.1.e-10) then
              jkk=jkk+nnz
            else
c             now step thru ions (jkk is ion index)
              klion=12
              mlion=npfirst(klion)
              jkk=0
              kl=0
              do while ((mlion.ne.0).and.(kl.lt.nnz))
                jkk=jkk+1
C               retrieve ion name from kdati
                mlm=mlion-1
                call dread(ltyp,lrtyp,lcon,
     $            nrdt,rdat,nidt,idat,nkdti,kdati,mlm,
     $            idat1,rdat1,kdat1,nptrs,0,lun11)
C               if not accessing the same element, skip to the next element
                mlleltp=idat(nidt-1)
                if (mlleltp.eq.mllel) then
                  kl=kl+1
                  if (lpril.ne.0)
     $              write (lun11,*)'  ion:',kl,jkk,mlion,mlleltp,
     $                          (kdati(mm),mm=1,nkdti)
                  do mm=nkdti+1,20
                    kdati(mm)=kblnk
                    enddo
                  call func2l(jkk,lpril,lun11,t,xee,xpx,
     $              idat1,rdat1,kdat1,nptrs,
     $              npar,npnxt,npfi,
     $              rnissl,rlev,ilev,
     $              nlpt,iltp,nlev,klev)
                  do mm2=1,nlev
                    mmtmp=npilev(mm2,jkk)
                    if (mmtmp.ne.0) then
                      kkkl=mmtmp
                      mmlv=mmtmp
                      if (xilev(kkkl).gt.-1.e-24) then
                        eth=rlev(1,mm2)
                        nions=nions+1
                        mllev(nions)=idat(nidt-1)
C                       Note that rwrk1 must be written to the file before
C                       it is overwritten in subsequent columns
                        rwrk1(nions)=xilev(mmlv)
                        rwrk2(nions)=rniss(mmlv)
                        elev(nions)=eth
                        ntptr(nions)=kkkl
                        natomic(nions)=nnz
                        nupper(nions)=mm2
                        write(kion(nions),'(10a1)')(kdati(mm),mm=1,10)
                        write(klevt(nions),'(20a1)')
     $                        (klev(mm,mm2),mm=1,20)
                        if (lpri.ne.0) then
                          write (lun11,*)nions,xilev(mmlv),
     $                         rdat(1),nnz,mmlv,kkkl
                          write (lun11,9296)kkkl,(kdati(mm),mm=1,20),
     $                    (klev(lk,mm2),lk=1,20),eth,xilev(kkkl),
     $                     rniss(kkkl),bilev(kkkl)
 9296                     format (1x,i6,1x,(40a1),7(1pe13.5))
                          endif
                        endif
                      endif
                    enddo
                  endif
C               Go to next ion
                mlion=npnxt(mlion)
                enddo
            endif
          mlel=npnxt(mlel)
C         Go to next element
          enddo

c

C     End of atomic database extraction
C----------------------------------------------------------------
C     define parameters for the binary table (see the above data statements)
      nrows=nions
      tfields=9
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)

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
      call ftpclj(unit,colnum,frow,felem,nions,ntptr,status)
      if (status .gt. 0)call printerror(lun11,status)

C     column  2 (Level number of this ion)
      colnum=2
      if (lpri.ne.0)
     $ write(lun11,*)'fstepr: Writing Column ',colnum
      call ftpclj(unit,colnum,frow,felem,nions,mllev,status)
      if (status .gt. 0)call printerror(lun11,status)

C     column  3  (Energy)
      colnum=3
      if (lpri.ne.0)
     $ write(lun11,*)'fstepr: Writing Column ',colnum
      call ftpcle(unit,colnum,frow,felem,nions,elev,status)
      if (status .gt. 0)call printerror(lun11,status)


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


C     column  5  (Atomic Number)
      colnum=5
      if (lpri.ne.0)
     $ write(lun11,*)'fstepr: Writing Column ',colnum
      call ftpclj(unit,colnum,frow,felem,nions,natomic,status)
      if (status .gt. 0)call printerror(lun11,status)


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

C----------------------------------------------------------------
C     column 7 (Level population)
C     rwrk1 can be safely overwritten after this step

      colnum=7
      if (lpri.ne.0)
     $ write(lun11,*)'fstepr: Writing Column ',colnum
      call ftpcle(unit,colnum,frow,felem,nions,rwrk1,status)
      if (status .gt. 0)call printerror(lun11,status)


      colnum=8
      if (lpri.ne.0)
     $ write(lun11,*)'fstepr: Writing Column ',colnum
      call ftpcle(unit,colnum,frow,felem,nions,rwrk2,status)
      if (status .gt. 0)call printerror(lun11,status)


C     column  9 (upper level index)
      colnum=9
      if (lpri.ne.0)
     $ write(lun11,*)'fstepr: Writing Column ',colnum
      call ftpclj(unit,colnum,frow,felem,nions,nupper,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
