subroutine mpartf(jatom,ion,indmol,t,u)
c     =======================================
c
c     yields partition functions with polynomial data from
c     ref. Irwin, A.W., 1981, ApJ Suppl. 45, 621.
c     ln u(temp)=sum(a(i)*(ln(temp))**(i-1)) 1<=a<=6
c
c     Input:
c       jatom = element number in periodic table
c       ion   = 1 for neutral, 2 for once ionized and 3 for twice ionized
c       indmol= index of a molecular specie (Tsuji index)
c       temp  = temperature
c     Output:
c       u     = partf.(linear scale) for iat,ion, or indmol, and temperature t
c
c
      implicit real*8 (a-h,o-z)
      real*8 a(6,3,92),aa(6),am(6,500)
      dimension indtsu(324),irw(500)
      save iread,a,am
c     data indtsu / 2,  5, 12, 4, 8, 7, 6,
c    *              9, 11, 10, 29, 50, 59, 46, 132, 52, 19,
c    *             13, 42, 38, 39, 37, 44, 36, 14, 118, 33,
c    *              3, 16, 57, 32, 49, 60, 54, 41, 107,  0,
c    *            148, 152, 153, 155, 0, 17, 24, 25, 28, 51,
c    *            112, 119,   0,   0,21, 15, 43, 56,  0, 64,
c    *             47,  65,   0,  61, 0, 62,118, 40, 66/
c     data indtsu / 2,  5, 12, 4, 8, 7, 6,
c    *              9, 11, 10, 29, 50, 59, 46, 132,  52, 19,
c    *             13, 42, 38, 39, 37, 44, 36,  14, 117, 33,
c    *              3, 16, 57, 32, 49, 60, 54,  41, 106,303,
c    *            147, 151, 152, 154, 302, 17,  24,  25, 28, 51,
c    *            111, 118, 102,   0,  21, 15,  43,  56,478, 64,
c    *             47,  65, 413,  61, 190, 62 ,108,  40, 66,214,
c    *             257*0./
      data indtsu / 2,  5, 12, 4, 8, 7, 6,
     *              9, 11, 10, 29, 50, 59, 46, 133,  52, 19,
     *             13, 42, 38, 39, 37, 44, 36,  14, 118, 33,
     *              3, 16, 57, 32, 49, 60, 54,  41, 107,304,
     *            148, 152, 153, 155, 303, 17,  24,  25, 28, 51,
     *            112, 119, 102,   0,  21, 15,  43,  22,478, 64,
     *             47,  65, 414,  61, 191, 62 ,109,  40, 66,214,
     *             120*0,    30, 136*0/

      data iread /0/
c
c     read data if first call:
c
      if(iread.ne.1) then
        if(irwtab.eq.0) then
           open(67,file= './data/irwin_orig.dat',status='old')
           nummol=66
         else
           open(67,file= './data/irwin_bc.dat',status='old')
           nummol=324
        end if
        read(67,*)
        read(67,*)
        do j=1,92
          do i=1,3
            if(j.eq.1.and.i.eq.3) goto 10
            sp=float(j)+float(i-1)/100.
            read(67,*) spec,aa
            do k=1,6
                a(k,i,j)=aa(k)
            end do
   10    continue
         end do
       end do 
c
       read(67,*)
       read(67,*)
       read(67,*)
       do i=1,500
          irw(i)=0
       end do
       do i=1,nummol
          read(67,*,end=15) spec,aa
          indm=indtsu(i)
          if(indm.gt.0) then
             irw(indm)=i
             do j=1,6
                am(j,indm)=aa(j)
             end do
          end if
        end do
   15   continue
        close(67)
        iread=1
      endif
c
c     evaluation of the partition function
c     stop if T is out of limits of Irwin's tables
c
        if(t.lt.1000.) then
          stop 'partf; temp<1000 K'
        else if(t.gt.16000.) then
          stop 'partf; temp>16000 K'
        endif
        tl=log(t)
        u=0.
c
c     atomic species
c
      if(jatom.gt.0.and.ion.gt.0) then
        ulog=    a(1,ion,jatom)+
     *       tl*(a(2,ion,jatom)+
     *       tl*(a(3,ion,jatom)+
     *       tl*(a(4,ion,jatom)+
     *       tl*(a(5,ion,jatom)+
     *       tl*(a(6,ion,jatom))))))
        if(jatom.eq.5.and.ion.eq.3) ulog=1.
        u=exp(ulog)
      end if
c
c     molecular species
c
      if(indmol.gt.0) then
        indm=indmol
        if(irw(indm).gt.0) then
           ulog=    am(1,indm)+
     *          tl*(am(2,indm)+
     *          tl*(am(3,indm)+
     *          tl*(am(4,indm)+
     *          tl*(am(5,indm)+
     *          tl*(am(6,indm))))))
           u=exp(ulog)
c       if(t.gt.5128..and.t.lt.5129)
c    *  write(6,631) t,indmol,indm,u
c 631   format('mpartf',f10.1,2i5,f16.3)
        end if
      end if
      return
      end