program isocoh implicit double precision (a-g,o-y),integer (h-n) parameter (nsp=10,l2=5,h5=7) character specie(nsp)*3, vname*8, y*1 integer ier, igo, ins(nsp), ifug, irk double precision nc, nh, no, ns, nc0, no0, nh0 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 * / cst100 /dlnfo2,elag,gz,gy,gx,ibuf,hu,hv,hw,hx * / cst24 /ipot,jv(l2),iv(l2) data tentoe, fo2, fs2, specie /2.302585093d0, 0.d0, 0d0, * 'H2O','CO2','CO ','CH4','H2 ','H2S','O2 ', * 'SO2','COS','UNK'/ c----------------------------------------------------------------------- open (20,file='iso.dat') open (21,file='spe.dat') write (20,*) ' 1' write (20,*) ' 0 0 ' write (20,*) 'gonk' write (20,*) '' write (20,*) '' write (20,*) '' write (20,*) ' 2' write (20,*) '1000. 200. 10000. 0. ' write (20,*) ' 6 T(K) ' write (20,*) ' 6P(bar) ' write (21,*) ' 1' write (21,*) ' 0 0 ' write (21,*) 'bonk' write (21,*) '' write (21,*) '' write (21,*) '' write (21,*) ' 2' write (21,*) '1000. 200. 1. 0. ' write (21,*) ' 6 T(K) ' write (21,*) ' 6 things ' 50 igo = 0 elag = 0. c initialize species fractions/volumes vname = ' X(CO2) ' c get the users choice of EoS: call rfluid (irk, ifug, 6) vname = ' X(O) ' isp = 5 do 1 i = 1, 5 1 ins(i) = i write (*,*) 'xo, pmin, tmin, tmax' read (*,*) xo, pmin, tmin, tmax write (*,*) 'enter dt, tol' read (*,*) dt, tol c call fluid routine: 10 write (*,*) 'p(bar) at t(k) =',tmax,' ?' read (*,*) p pmax = p if (p.eq.0.) goto 99 t = tmax call cfluid (fo2, fs2) no0 = xs(1) + xs(2)*2.+ xs(3) + xs(7)*2. * + xs(8)*2. + xs(9) nc0 = xs(2) + xs(3) + xs(4) + xs(9) nh0 = (xs(1) + xs(5) + xs(6))*2. + xs(4)*4. d = (no0 * 16. + nc0 * 12.011 + nh0 * 1.008)/vol npt = 1 call inpt (d,npt) write (*,*) p, t, d v0 = vol igo = 0 100 t = t - dt if (igo.eq.1) goto 200 if (t.lt.tmin) then t = tmin igo = 1 end if if (p.lt.pmin) goto 99 c for reduced t c find isochore: dp = (p - pmin) / 2. 110 p = p - dp call cfluid (fo2,fs2) no = xs(1) + xs(2)*2.+ xs(3) + xs(7)*2. * + xs(8)*2. + xs(9) nc = xs(2) + xs(3) + xs(4) + xs(9) nh = (xs(1) + xs(5) + xs(6))*2. + xs(4)*4. gr = no * 16. + nc * 12.011 + nh * 1.008 d = gr/vol c constant O volume: rat = no0/no vc = vol * rat c bulk volume: vt = vc + (nc0 - nc*rat) * 5.3 dv = vt - v0 if (dabs(dv).lt.tol) then c converged write (*,*) p,t,d if (nc0 - nc*rat.lt.0.) write (*,*) 'ate it' npt = npt + 1 call inpt (d,npt) goto 100 else if (vt.lt.v0) then c p to high if (dp.lt.0.) dp = -dp/2. goto 110 else c vt > v0, p too low if (dp.gt.0.) dp = -dp/2. goto 110 end if c now calculate constant speciation isochore: 200 call outpt (npt,0) npt = 0 t = tmin if (ifug.eq.10) then call nogcoh (fo2) else call mrkmix (ins,isp) end if v0 = vol write (*,*) 'p, t, d = ',p,t,gr/vol t = tmax + dt p = 2.*pmax c calculate isochore from tmax to tmin: igo = 0 220 t = t - dt if (igo.eq.1) goto 230 if (t.lt.tmin) then t = tmin igo = 1 end if if (p.lt.pmin) goto 230 c for reduced t c find isochore: dp = (p - pmin) / 2. 210 p = p - dp call nogcoh (fo2) dv = v0 - vol if (dabs(dv).lt.tol) then c converged write (*,*) p,t,gr/vol npt = npt + 1 call inpt (gr/vol,npt) goto 220 else if (vol.lt.v0) then c p to high if (dp.lt.0.) dp = -dp/2. goto 210 else c vt > v0, p too low if (dp.gt.0.) dp = -dp/2. goto 210 end if 230 call outpt (npt,1) goto 10 1000 format (/,' Enter p(bar), T(K), ',a8': ') 1010 format (/,' Enter p(bar), T(K): ') 1110 format (/,' Enter the log10[f(S2)]:',/) 1120 format (/,' Enter a zero for pressure to quit.',/) 1130 format (/,10x,'f(H2O) = ',g12.5,/,10x,'f(CO2) = ',g12.5,/) 1140 format (/,10x,'f(CO2) = ',g12.5,/) 1150 format (/,10x,'f(H2O) = ',g12.5,/ * ,10x,'f(CO2) = ',g12.5,/ * ,10x,'log[f(O2)] = ',g12.5,/) 1160 format (/,10x,'f(H2O) = ',g12.5, * /,10x,'f(H2 ) = ',g12.5, * /,10x,'log[f(O2)] = ',g12.5,/) 1170 format (/,10x,'log[f(O2)] = ',g12.5, * /,10x,'log[f(S2)] = ',g12.5, * /,10x,' a(gph) = ',g12.5,/) 1180 format (10x,4(5x,a3,5x)) 1190 format (5x,'x',4x,4(g12.5,1x)) 1200 format (5x,'f',4x,4(g12.5,1x),/) 1210 format (/,' Change EoS, buffer, or graphite activity (y/n)?',/) 1220 format (a1) 1230 format (22x,'Speciation/Fugacities',/) 1240 format (/,22x,'Atomic Proportions/Volume',//, * 10x,'C',12x,'H',12x,'O',12x,'S',6x,'V(cm3/mol species)') 1250 format (4x,5(g12.5,1x),/) 1260 format (5x,'v',4x,4(g12.5,1x),/) 1270 format (5x,'Back-calculated X(O) = ',g16.9) 1280 format (/,10x,'p(bar) = ',g12.5,/,10x,'T(K) = ',g12.5) 1290 format (5x,'S/H = ',g12.5,/) 1300 format (10x,'V(cm3/mol) = ',g12.5,/) 1310 format (5x,'Back-calculated X(S) = ',g16.9) 1320 format (5x,'a(gph) = ',g12.5,/) 1330 format (/,5x,'COMPOSITION IS SUPERSATURATED WITH ', * 'RESPECT TO GRAPHITE!!',/) 1340 format (/,' Enter p(bar), T(K), X(O), X(S): ') 1350 format (5x,'Back-calculated X(C) = ',g16.9) 1360 format (/,' Enter p(bar), T(K), X(O), X(C): ') 1370 format (/,5x,'Sum of species fractions: ',f14.9,/) 1380 format (/,5x,'INVALID SPECIIATION!!',/) 99 write (20,*) '1 1 1 1 1 1 1 1 1 1 1 1 ' write (20,*) '0 0 0 0 0 ' write (21,*) '1 1 1 1 1 1 1 1 1 1 1 1 ' write (21,*) '0 0 0 0 0 ' close (20) close (21) 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 subroutine inpt (d,npt) implicit double precision (a-g,o-y),integer (h-n) parameter (nsp=10,l2=5,h5=7) integer ier, igo, ins(nsp), ifug, irk common / cst5 /p,t,xo,u(6) * / cstcoh /xs(nsp),g(nsp),v(nsp) common / crap /xi(100,nsp),den(100),pi(100),ti(100),cc(100) den(npt) = d pi(npt) = p ti(npt) = t cc(npt) = xs(2)/(xs(2) + xs(4)) do i = 1, 5 xi(npt,i) = xs(i) end do end subroutine outpt (npt,igo) implicit double precision (a-g,o-y),integer (h-n) parameter (nsp=10) character specie(nsp)*3 integer igo common / crap /xi(100,nsp),den(100),pi(100),ti(100),cc(100) data specie / * 'H2O','CO2','CO ','CH4','H2 ','H2S','O2 ', * 'SO2','COS','UNK'/ npt2 = 2*npt write (20,2010) npt2,1,1,1,1,1,'isokor' write (20,2020) (ti(l),pi(l),l=1,npt) goto (99), igo write (21,2010) npt2,1,1,1,1,1,'density' write (21,2020) (ti(l),den(l),l=1,npt) write (21,2010) npt2,1,1,1,1,1,'density' write (21,2020) (ti(l),den(l),l=1,npt) write (21,2010) npt2,1,1,1,1,1,'ch4/(co2+ch4)' write (21,2020) (ti(l),cc(l),l=1,npt) do i = 1, 5 write (21,2010) npt2,1,1,1,1,1,specie(i) write (21,2020) (ti(l),xi(l,i),l=1,npt) end do 2010 format (i5,5(i2,1x),a10) 2020 format (6(g12.6,1x)) 99 end