program specis

      implicit double precision (a-g,o-y),integer (h-n)

      parameter (nsp=10,l2=5,h5=7)

      character specie(nsp)*3, vname*8, yes*1

      integer  igo, ins(nsp), ifug, irk

      double precision 
     *       x(nsp,1000),f(1000),y(1000),ns(1000),no(1000),
     *       nh(1000),nc(1000)


      common / cst11 /fh2o,fco2/ cst5 /p,t,xo,u(6)
     *       / cstcoh /xs(nsp),g(nsp),v(nsp)/ cst26 /vol
     *       / cst10 /iff(2),idss(h5),ifug,ifyn,isyn
     *       / cst24 /ipot,jv(l2),iv(l2)
     *       / cst100 /dlnfo2,elag,gz,gy,gx,ibuf,hu,hv,hw,hx

      data tentoe, fo2, fs2, specie /2.302585093d0, 0.d0, 0d0,
     *      'H2O','CO2','CO ','CH4','H2 ','H2S','O2 ',
     *      'SO2','COS','UNK'/

c-----------------------------------------------------------------------

      jgo = 0
      fs2 = -9999.*tentoe/2.
      elag = 0.
c                                  initialize species fractions/volumes
      do 60 i = 1, nsp
60       xs(i) = 0.

      fmin = 1.d60
      fmax = -1.d60
      ymin = 1.d60
      ymax = -1.d60
      open (10,file='species.dat')

      write (10,*) '  1'
      write (10,*) '  0   0 '
      write (10,2020) p, t, xonc   
      write (10,*) ''
      write (10,*) ''
      write (10,*) ''
      write (10,*) '  2'
c                                  get the users choice of EoS:   
      call rfluid (irk, ifug, 6)

50    write (*,*) ' p, t, xmin, xmax, xinc, isp1, isp2?'
      read (*,*) p, t, xmin, xmax, xonc, isp1, isp2

      xo = xmin
      k = 0
c                                  for multispecies fluids set
c                                  up species indices:
      if (ifug.gt.6.and.ifug.lt.13) then
         vname = '  X(O)  '
         if (ifug.eq.7.or.ifug.eq.8) vname = 'log(fo2)'
         isp = 5
         do 1 i = 1, 6
1           ins(i) = i
         if (ifug.eq.12) then
            isp = 8
            ins(7) = 8
            ins(8) = 9
         end if
      else if (ifug.eq.16) then
         vname = '  X(O)  '
         isp = 3
         ins(1) = 1
         ins(2) = 5
         ins(3) = 7
      else if (ifug.eq.17) then
         vname = '  X(O)  '
         isp = 4
         ins(1) = 1
         ins(2) = 5
         ins(3) = 6
         ins(4) = 8
      else if (ifug.eq.13.or.ifug.eq.15) then 
         vname = ' X(H2)  '
      end if        
c                                  call fluid routine:
10    call cfluid (fo2, fs2)

      k = k + 1
      f(k) = p*g(isp1)*xs(isp1)
      y(k) = p*g(isp2)*xs(isp2)

      if (f(k).lt.fmin) fmin = f(k)
      if (f(k).gt.fmax) fmax = f(k)
      if (y(k).lt.ymin) ymin = y(k)
      if (y(k).gt.ymax) ymax = y(k)

      write (*,2000) f(k),y(k),100.*xo

      xo = xo + xonc

      if (xo.gt.xmax) goto 20

      goto 10
c                                  this is the max X(O), min X(O),
c                                  max log(Y), min log(Y) for
c                                  psvdraw
20    if (jgo.eq.0) then 
         write (10,*) fmax,fmin, ymax,ymin
         write (10,*) ' 6',specie(isp1)
         write (10,*) ' 4',specie(isp2)
         jgo = 1
      end if 

      write (10,2010) 2*k,j,1,1,1,1
      write (10,2020) (f(l),y(l),l=1,k)

2000  format (8(g12.6,2x))
2010  format (i5,5(i2,1x),a8)
2020  format (6(g12.6,1x))
1210  format (/,' Change p,t,x, etc (y/n)?',/)
1220  format (a1)

99    write (*,1210) 
      read (*,1220) yes
      if (yes.eq.'y'.or.yes.eq.'Y') goto 50

      write (10,*) ' 1   1   1   1   1 1.0'
      write (10,*) ' 0   0  '

      end 
 
      block data
c-----------------------------------------------------------------------
      implicit double precision (a-g,o-y)
 
      common/ cst5 /p,t,xco2,u1,u2,tr,pr,r,ps/ cst85 /pp,tt,yy,rr
      data r,rr/8.3144126d0,83.14d0/
      end