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),fug(nsp,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 50 igo = 0 ipl = 1 fs2 = -9999.*tentoe/2. elag = 0. c initialize species fractions/volumes do 60 i = 1, nsp 60 xs(i) = 0. c get the users choice of EoS: call rfluid (irk, ifug, 6) write (*,*) 'plot (log) species fractions (y/n)?' write (*,*) '(if you answer no fugacities are plotted).' read (*,1220) yes if (yes.eq.'y'.or.yes.eq.'Y') ipl = 0 igo = 0 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 write (*,2900) read (*,*) p, t write (*,2910) vname read (*,*) xmn, xmx, xonc write (*,2920) read (*,*) fmn, fmx 2900 format (' Enter pressure (bar), temperature (K):') 2910 format (' Enter min, max, and increments for the', * ' independent variable (',a8,'):') 2920 format (' Enter min and max limits for the dependent variables,',/, * ' in the plot file, values that exceed these limits will' * ' be set to the',/,' limits established here:') 2930 format (' 8',a8) if (jgo.eq.0) then 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 this is the max X(O), min X(O), c max log(Y), min log(Y) for c psvdraw write (10,*) xmx, xmn, fmx, fmn write (10,2930) vname write (10,*) ' 4log(Y) ' end if write (*,3000) vname,(specie(ins(i)),i=1,isp) write (*,3040) if (ipl.eq.0) then write (*,3030) else write (*,3010) end if 3000 format (/,' values which follow are ',a8, * 'and the mole fractions of:',/,1x,10(a3,1x)) 3030 format (/,' the plot file includes the log of the ' * ,' C-O-H-S atomic proportions and ', * /,' the log of the species fractions.'/) 3010 format (/,' the plot file includes the log of the ' * ,' C-O-H-S atomic proportions and ', * /,' the log of the species fugacities.'/) 3040 format (' and log(fO2).',/) xo = xmn if (xo.eq.0.d0) xo = 0.00000001 c call fluid routine: 10 call cfluid (fo2, fs2) write (*,2000) xo,(xs(ins(i)),i=1,isp),fo2/tentoe k = k + 1 do 90 i = 1, isp j = ins(i) if (xs(j).gt.1.d0.or.xs(j).lt.0.d0) xs(j) = 1.d0 fug(j,k) = dlog10(xs(j)*g(j)*p) x(j,k) = dlog10(xs(j)) if (fug(j,k).lt.fmn) then fug(j,k) = fmn else if (fug(j,k).gt.fmx) then fug(j,k) = fmx end if if (x(j,k).lt.fmn) then x(j,k) = fmn else if (x(j,k).gt.fmx) then x(j,k) = fmx end if 90 continue f(k) = fo2/tentoe if (f(k).lt.fmn) f(k) = fmn if (f(k).gt.fmx) f(k) = fmx y(k) = xo ns(k) = xs(6) + xs(8) + xs(9) + xs(10) no(k) = xs(1) + xs(2)*2.+ xs(3) + xs(7)*2. * + xs(8)*2. + xs(9) + xs(10)*3. nc(k) = xs(2) + xs(3) + xs(4) + xs(9) nh(k) = (xs(1) + xs(5) + xs(6))*2. + xs(4)*4. tot = ns(k) + no(k) + nh(k) + nc(k) ns(k) = ns(k) / tot no(k) = no(k) / tot nc(k) = nc(k) / tot nh(k) = nh(k) / tot if (ns(k).eq.0.) ns(k) = 1. if (no(k).eq.0.) no(k) = 1. if (nh(k).eq.0.) nh(k) = 1. if (nc(k).eq.0.) nc(k) = 1. xfmx = 10.**fmx xfmn = 10.**fmn if (ns(k).lt.fmn) ns(k) = xfmn if (ns(k).gt.fmx) ns(k) = xfmx if (no(k).lt.fmn) no(k) = xfmn if (no(k).gt.fmx) no(k) = xfmx if (nh(k).lt.fmn) nh(k) = xfmn if (nh(k).gt.fmx) nh(k) = xfmx if (nc(k).lt.fmn) nc(k) = xfmn if (nc(k).gt.fmx) nc(k) = xfmx if (xo.ge.0.9999999d0.and.xo.le.xmx) goto 20 xo = xo + xonc goto (20), igo if (xo.gt.xmx) then xo = xmx if (xo.eq.1.d0) xo = 0.9999999d0 igo = 1 end if goto 10 20 do 140 i = 1, isp j = ins(i) write (10,2010) 2*k,j,1,1,1,1,specie(j) if (ipl.eq.0) then write (10,2020) (y(l),x(j,l),l=1,k) else write (10,2020) (y(l),fug(j,l),l=1,k) end if 140 continue c O2 fugacity write (10,2010) 2*k,1,1,1,1,1 write (10,2020) (y(l),f(l),l=1,k) c bulk fractions: write (10,2010) 2*k,1,1,1,1,1,'C' write (10,2020) (y(l), dlog10(nc(l)),l=1,k) write (10,2010) 2*k,1,1,1,1,1,'H' write (10,2020) (y(l),dlog10(nh(l)),l=1,k) write (10,2010) 2*k,1,1,1,1,1,'O' write (10,2020) (y(l),dlog10(no(l)),l=1,k) write (10,2010) 2*k,1,1,1,1,1,'S' write (10,2020) (y(l),dlog10(ns(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 EoS, buffer, or graphite activity (y/n)?',/) 1220 format (a1) 99 write (*,1210) read (*,1220) yes if (yes.eq.'y'.or.yes.eq.'Y') then jgo = 1 goto 50 end if 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