c Copyright (c) 1990 by James A. D. Connolly, Institute for Mineralogy and c Petrography, Swiss Federal Insitute of Technology, CH-8092 Zurich, c SWITZERLAND. All rights reserved. c Please do not distribute this source. c------------------------------------------------------------------------ c FRENDLY - a thermodynamic calculator for petrologic problems. c------------------------------------------------------------------------ program FRNDLY implicit double precision (a-g,o-y),integer (h-m,z) parameter (l2=5,k6=7,h5=7) character*80 vname*8,gname*78,xname*18,uname*8,y*1,rxny*1 character*20 fname double precision pp(2),tt(2),xx(2),gg(2),ee(2),ss(2),vv(2),ccp(2), * uu(2),aa(2) common/ cst5 /p,t,xco2,u1,u2,tr,pr,r,ps * / cst9 /vmax(l2),vmin(l2),dv(l2)/ cst113 /ifproj * / cst100 /dlnfo2,elag,gz,gy,gx,ibuf,hu,hv,hw,hx * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn * / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5 * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / csta2 /xname(k6),vname(l2),gname(l2) c----------------------------------------------------------------------- c assign data files call fopen idiag = 0 idid = 0 ifproj = 0 write (*,1040) read (*,1000) uname write (*,1050) uname read (*,1020) y if (y.ne.'y'.and.y.ne.'Y') then write (*,1060) uname goto 99 end if 20 write (*,1030) read (*,*,iostat=ier) icopt call rerror (ier,*20) if (icopt.eq.5) goto 99 call input2 (icopt,rxny,uname) if (icopt.eq.1) then c select variables and set up graphics c file: if (idiag.eq.0) then call vars (loop) end if idiag = 1 idid = 1 10 call eqrxn (loop) write (*,1070) read (*,1020) rxny if (rxny.eq.'y'.or.rxny.eq.'Y') then call input2 (icopt,rxny,uname) goto 10 end if else if (icopt.eq.2) then write (*,1180) uname read (*,1020) y if (y.ne.'y'.and.y.ne.'Y') goto 6 5001 write (*,1190) read (*,*,iostat=ier) pmin, pmax, pinc call rerror (ier,*5001) 5002 write (*,1200) read (*,*,iostat=ier) tmin, tmax, tinc call rerror (ier,*5002) if (ifyn.eq.0) then 5003 write (*,1210) read (*,*,iostat=ier) xmin, xmax, xinc call rerror (ier,*5003) loop3 = idint ((xmax-xmin)/xinc) + 1 else loop3 = 1 end if loop2 = idint ((tmax-tmin)/tinc) + 1 loop1 = idint ((pmax-pmin)/pinc) + 1 fname = ' ' write (*,1220) read (*,1230) fname write (*,1240) open (n4, file=fname) write (n4,*) 10 write (n4,*) tinc, pinc, tmin, pmin write (n4,*) loop2, loop1, loop3 t = tmin do 70 i = 1, loop2 p = pmin do 80 j = 1, loop1 xco2 = xmin do 90 k = 1, loop3 call props (g,e,u,s,v,cp) write (n4,1250) p,t,xco2,g,e,s,v,cp,-g/r/t, * -g/r/t/2.302585093 90 xco2 = xco2 + xinc 80 p = p + pinc 70 t = t + tinc close (n4) goto 20 6 write (*,1080) uname,uname,uname,uname 5 write (*,1100) read (*,*,iostat=ier) p,t call rerror (ier,*5) if (p.eq.0.0) goto 20 if (ifyn.eq.0) then 51 write (*,1110) read (*,*,iostat=ier) xco2 call rerror (ier,*51) end if c call props (g,e,u,s,v,cp) c write (*,1120) t,p,g/1000.,e/1000., * (g-p*v)/1000.,u/1000.,s,v,cp, * -g/r/t,-g/r/t/2.302585093 c write (*,1090) read (*,1020) y if (y.ne.'y'.and.y.ne.'Y') goto 5 call change c goto 5 else if (icopt.eq.3) then c write (*,1140) c 40 do 41 i = 1, 2 52 write (*,1150) i,i read (*,*,iostat=ier) pp(i),tt(i) call rerror (ier,*52) p = pp(i) t = tt(i) c if (pp(i).eq.0.0) goto 20 c if (ifyn.eq.0) then 53 write (*,1160) read (*,*,iostat=ier) xx(i) call rerror (ier,*53) xco2 = xx(i) end if c call props * (gg(i),ee(i),uu(i),ss(i),vv(i),ccp(i)) aa(i) = gg(i)-pp(i)*vv(i) 41 continue write (*,1170) (gg(2)-gg(1))/1000.,(ee(2)-ee(1))/1000., * (aa(2)-aa(1))/1000., * (uu(2)-uu(1))/1000., * ss(2)-ss(1), * vv(2)-vv(1),ccp(2)-ccp(1) write (*,1090) read (*,1020) y if (y.ne.'y'.and.y.ne.'Y') goto 40 call change goto 40 else if (icopt.eq.4) then call nentry end if goto 20 99 write (*,1130) uname if (idid.eq.1) write (n4,1010) 1,1,1,1,1,1.,0,0 stop 1000 format (a8) 1010 format (5(1x,i1),1x,f3.1,/,2(1x,i1)) 1020 format (a1) 1030 format (/,' Choose from following options:', * //,2x,'1) calculate equilibrium coordinates for a reaction.', * /,2x,'2) calculate thermodynamic properties for a phase or', * ' reaction relative to',/,5x,'the reference state.', * /,2x,'3) calculate change in thermodynamic properties ', * ' from one p-t-x',/,5x,'condition to another.', * /,2x,'4) create new thermodynamic data file entries.', * /,2x,'5) quit.', * //,' With options 1-3 you may also modify', * ' thermodynamic data,',/,' modified data can then', * ' be stored as a new entry in the thermodynamic data ',/, * ' file.',/) 1040 format (/,' Hi! i am a user freindly program.',/, * ' what is your name? ') 1050 format (/,' Hi ',a8,'!, i like you and i hope we have fun!',/, * ' can you say "therm-o-dy-nam-ics" (y/n)? ') 1060 format (/,' Weeell ', a8,' why dont you go read Gibbs and try me', * ' again later.') 1070 format (' Calculate a different equilibrium (y/n)? ') 1080 format (/,' now ',a8,' there is one thing i want you to remember', * ' and that is that the',/,' values for G and H that i', * ' calculate for you are "apparent" free energies and', * /,' enthalpies. ',a8,' if you dont know what that means,' * ,' read Helgeson et al. 1978.',//,1x,a8, * ' now i am going to prompt you for conditions', * ' at which i should',/,' calculate thermodynamic', * ' properties, when you want to stop just enter', * ' zeroes.',//, * ' ok, here we go, it is fun time ',a8,'!!',//) 1090 format (' Modify or output', * ' thermodynamic parameters of a phase (y/n)? ') 1100 format (' Enter p(bars) and t(k) (zeroes to quit): ') 1110 format (' Enter X(CO2/O) in fluid phase: ') 1120 format (/,' At ',g13.6,'k and ',g13.6, * 'bar:',/,6x,'G(kj) =',g15.8,/,6x,'H(kj) =',g15.8,/, * 6x,'A(kj) =',g15.8,/, * 6x,'U(kj) =',g15.8,/, * 6x,'S(j/k) =',g13.6,/,6x,'V(j/bar)=',g13.6,/, * 6x,'Cp(j/k) =',g13.6,/, * 6x,'loge K =',g13.6,/, * 6x,'log10 K =',g13.6,/) 1130 format (/,' Have a nice day ',a8,'!',/) 1140 format (/,' Option to calculate change in ', * 'thermodynamic properties from',/, * ' p(1)-T(1)-X(CO2/O)(1) to p(2)-T(2)-X(CO2/O)(2)',/, * ' enter zeros to quit.',//) 1150 format (' Enter p(',i1,') (bars) and T(',i1,') (K) (0 to quit):') 1160 format (' Enter X(CO2/O)(',i1,'): ') 1170 format (//,6x,'delta G(kj) =',g15.8,/, * 6x,'delta H(kj) =',g15.8,/, * 6x,'delta A(kj) =',g15.8,/, * 6x,'delta U(kj) =',g15.8,/, * 6x,'delta S(j/k) =',g13.6,/, * 6x,'delta V(j/bar)=',g13.6,/, * 6x,'delta cp(j/k) =',g13.6,/, * 6x,' loge K =',g13.6,/, * 6x,' log10 K =',g13.6,/) 1180 format (1x,a8,' write a properties table (Y/N)?',/) 1190 format (' Enter min, max, and increment for P(bar):',/) 1200 format (' Enter min, max, and increment for T(K):',/) 1210 format (' Enter min, max, and increment for X(CO2/O moles):',/) 1220 format (' Enter file name for table (<20 characters):',/) 1230 format (a20) 1240 format (/,' Table row entries will be:',/, * ' P, T, X(CO2/O), G(j),', * ' H(j), S(j/K), V(J/bar), Cp(j/K), ln K, log K',/) 1250 format (10(g13.6,1x)) end block data c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) character*28 list common/ cst5 /p,t,xco2,u1,u2,tr,pr,r,ps/ cst85 /pp,tt,yy,rr * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst206 /list(18) c----------------------------------------------------------------------- data r,rr/8.3144126d00,83.14/ c the following data statement assigns c logical unit numbers for i/o. data n1,n2,n3,n4,n5,n6,n7,n8,n9/21,22, 6,24,25,26, 5, 6,29/ c data list /'standard free energy g (j) ', * 'standard entropy s (j/k) ', * 'standard volume v (j/bar) ', * 'heat capacity coefficient a ', * 'heat capacity coefficient b ', * 'heat capacity coefficient c ', * 'heat capacity coefficient d ', * 'heat capacity coefficient e ', * 'heat capacity coefficient f ', * 'heat capacity coefficient g ', * 'volumetric coefficient b2 ', * 'volumetric coefficient b3 ', * 'volumetric coefficient b4 ', * 'volumetric coefficient b5 ', * 'volumetric coefficient b6 ', * 'volumetric coefficient b7 ', * 'activity ', * 'reaction coefficient '/ end subroutine chptx c---------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (l2=5,k6=7) character*8 vname*8,gname*78,xname*18 common/ cst63 /delv(l2) * / csta2 /xname(k6),vname(l2),gname(l2) * / cst5 /v(l2),tr,pr,r,ps/ cst9 /vmax(l2),vmin(l2),dv(5) * / cst24 /ipot,jv(l2),iv(l2) write (*,1020) do 10 i = 1, ipot j = iv(i) 20 write (*,1000) vname(j),vmin(j),vmax(j) read (*,*,iostat=ier) vmin(j),vmax(j) if ((i.eq.3.and.vmin(j).lt.0.0).or. * (i.eq.3.and.vmax(j).gt.1.0).or. * (j.ne.3.and.vmin(j).ge.vmax(j)).or.ier.ne.0) then write (*,1010) goto 20 end if v(j) = vmin(j) delv(j) = (vmax(j) - vmin(j)) / 40. 10 dv(j) = (delv(j)) / 40. call concrt 1000 format (/,' Enter new min/max values for ',a8,' (' * 'old values were ',g12.3,',',g12.3,')',/) 1010 format (/,' try again geek.',/) 1020 format (/,' This option does not change plot limits!' * ,' To do this, modify default plot options', * /,' while running PSVDRAW.',/) end subroutine input2 (icopt,rxny,uname) c---------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k0=25,k1=1200,k4=16,k5=12,l2=5,h5=7) parameter (m9=50,k9=50,m6=3,m7=9,m8=9,h8=100) character*5 cmpnt, dname*40 character*8 name, exname, uname character*1 y, rxny dimension germ(k4), gerd(m8), gerl(m7,m6) common/ cst1 /thermo(k4,k1),vf(2,k1),vs(h5,k1),uf(2),us(h5) * / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5 * / cst25 /vnu(11),idr(11),ivct/ cst205 /act(11),idf(2) * / cst201 /vuf(2),vus(h5),iffr,isr/ oxy /io2 * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn * / cst20 /comps(k1,k5) * / cst6 /icomp,istct,iphct,icp/ cst36 /exname(h8) * / cst5 /p,t,xco2,u1,u2,tr,pr,r,ps/ csta6 /name * / cst37 /iprct,ixct,iexyn,istbyn * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / csta5 /dname(8),cmpnt(k0)/ cst42 /ic(k5),idbase * / cst43 /therm(k4),comp(k0),atwt(k0),idh2o,idco2,ikind * / cst202 /tm(m7,m6),td(m8),ilam,idiso,lamin,idsin * / cst203 /therdi(m8,m9),therlm(m7,m6,k9) * / cst204 /ltyp(k1),lmda(k1),idis(k1) * / cst88 /vu(2,k1),jprct,jmct,jmyn c--------------------------------------------------------------------- rewind n2 call topn2 (1,icmpn) jmyn = 1 isyn = 1 c this check is here to avoid c redimensioning comps, but cause c the maximum number of components c counted by frendly to be k5, if you c want as many components in frendly as c in the data base you must set k0 = k5. jcmpn = icmpn icomp = icmpn if (icmpn.gt.k5) then icomp = k5 jcmpn = k5 write (*,3020) end if do 1 i = 1, icomp 1 ic(i) = i c list database phases: write (*,1030) read (*,5000) y if (y.ne.'y'.and.y.ne.'Y') goto 5 write (*,3000) (cmpnt(zi),zi=1,jcmpn) 300 call getphi (name,ibase,icmpn,*5) c reject wrong data base cards: if (idbase.ne.ibase) goto 300 c acceptable data, count the phase write (*,3010) name, (comp(zi),zi=1,jcmpn) goto 300 c begin the real work: c initialize: 5 if (icopt.eq.4) goto 99 do 60 k = 1, k5 60 comps(1,k)=0. do 62 k = 1, 11 lmda(k) = 0 idis(k) = 0 vnu(k) = 0. act(k) = 0. ltyp(k) = 0 62 continue do 80 k = 1, m7 do 90 l = 1, m6 tm(k,l) = 0.d0 90 gerl(k,l) = 0.d0 80 continue do 81 k = 1, m8 81 gerd(k) = 0.d0 iphct=0 io2 = 0 ifyn=1 jlam=1 jdis=1 vuf(1)=0.0d0 vuf(2)=0.0d0 idf(1)=0 idf(2)=0 if (icopt.eq.1) then rxny='y' 5002 write (*,4010) read (*,*,iostat=ier) isct call rerror (ier,*5002) else write (*,4000) read (*,5000) rxny if (rxny.eq.'y'.or.rxny.eq.'y') then 5007 write (*,4010) read (*,*,iostat=ier) isct call rerror (ier,*5007) else isct = 1 end if end if c make the user feel important: jj = isct+1 do 70 k = 1, k4 thermo(k,jj) = 0.d0 70 germ(k) = 0.d0 c do 10 i=1,isct 30 if (rxny.ne.'y'.and.rxny.ne.'Y') then write (*,4020) else write (*,4030) i end if read (*,4060) exname(1) c general input data for main program rewind n2 call eohead (n2) 220 call getphi (name,ibase,icmpn,*999) c reject wrong data base cards: if (idbase.ne.ibase) goto 220 c acceptable data, count the phase if (exname(1).ne.name) goto 220 c set special flag if O2 iphct = iphct + 1 if (exname(1).eq.'O2 ') io2 = iphct c store thermodynamic parameters: if (rxny.eq.'Y'.or.rxny.eq.'y') then 5003 write (*,1040) exname(1) read (*,*,iostat=ier) vvv call rerror (ier,*5003) else vvv = 1 end if if (exname(1).eq.cmpnt(idh2o)) then vuf(1)=vvv idf(1)=iphct ifyn=0 else if (exname(1).eq.cmpnt(idco2)) then vuf(2)=vvv idf(2)=iphct ifyn=0 end if c get activity coefficients: 5004 write (*,4070) exname(1) read (*,*,iostat=ier) act(iphct) call rerror (ier,*5004) c reaction coefficient: vnu(iphct)=vvv c sum super function: call loadit (iphct) c dumb function: do 100 k = 1, k4 100 germ(k) = germ(k) + vvv*therm(k) if (idiso.ne.0) then jdis=0 do 110 k = 1, m8 110 gerd(k) = gerd(k) + vvv*td(k) end if if (ilam.ne.0) then jlam=0 do 120 k= 1, m7 do 130 l= 1, m6 130 gerl(k,l) = gerl(k,l) + vvv*therlm(k,l,lmda(iphct)) 120 continue end if c next phase: goto 10 999 write (*,4050) uname goto 30 10 continue goto (61),ifyn c select equation of state for the c saturated phase. call rfluid (1,ifug,n3) c output standard state parm summations: 61 write (*,2000) write (*,2010) (germ(k),k=1,3) write (*,2020) write (*,2010) (germ(k),k=4,8) write (*,2040) write (*,2010) (germ(k),k=9,10) write (*,2050) write (*,2010) (germ(k),k=11,12) write (*,2060) write (*,2010) (germ(k),k=13,14) write (*,*) goto (140),jdis write (*,1050) write (*,2010) gerd 140 goto (160),jlam write (*,1060) do 150 k=1,3 150 write (*,*) (gerl(l,k),l=1,9) 160 do 170 j=1,k4 do 170 k=1,iphct 170 thermo(j,jj)=thermo(j,jj)+vnu(k)*thermo(j,k) write (*,2080) write (*,2010) (thermo(k,jj),k=1,k4) write (*,*) if (rxny.eq.'y'.or.rxny.eq.'Y') then 55 do 56 k=1,icmpn 56 comps(jj,k)=0.0 do 52 l=1,iphct do 53 k=1,icmpn 53 comps(jj,k)=comps(jj,k)+vnu(l)*comps(l,k) 52 continue itic=0 do 50 k=1,icmpn if (comps(jj,k).eq.0.) goto 50 itic=itic+1 50 continue if (itic.eq.0) goto 99 do 54 k=1,icmpn if (comps(jj,k).eq.0.d0) goto 54 write (*,2460) comps(jj,k),cmpnt(k) 54 continue write (*,1070) read (*,5000) y if (y.ne.'y'.and.y.ne.'Y') goto 99 call stoich goto 55 end if c position pointer to last record of c data file to allow writing: 99 read (n2,1170,end=9999) name goto 99 1000 format (3(a8,18x)) 1010 format (i2,a78) 1020 format (a8,i2,i2,i2,i2) 1030 format (/,' List database phases (y/n)? ') 1040 format (' Enter reaction coefficient for: ',a8, * ' products (+), reactants (-): ') 1050 format (' temperature dependent disordering:') 1060 format (' lambda transitions:') 1070 format (/,' Change reaction coefficients (y/n)? ') 1170 format (a40) 2000 format (/,1x,'standard state properties follow:',//, * 1x,'g(j/mole) s(j/mole k) v(j/mole bar)') 2010 format (5(1x,g15.8)) 2020 format (/,1x,'heat capacity (j/mole k) function:',//, * 1x,' a b*t c/t**2 ', * 'd*t**2 e/t**1/2') 2040 format (/,1x,' f/t g/t**3') 2050 format (/,1x,'volumetric (j/mole bar) function:',//, * 1x,' b2*(tr-t) b3*(exp(-t/300)-exp(tr/300) ') 2060 format (/,1x,' b4*(p-pr) b5*(exp(-p/35.d3)-exp(-pr/35.d3)') 2080 format (/,' super function for g (j/mole):',/) 2460 format (' Warning ** reaction does not balance by ',g13.6, * ' moles of ',a5) 3000 format (t30,'composition',/,' phase',4x,12(a5,1x)) 3010 format (1x,a8,1x,12(f5.2,1x)) 3020 format (/,' too many components, only first 12 ', * ' will be listed.',/) 4020 format (' Calculate thermodynamic properties for phase: ') 4030 format (' Enter phase or species number ',i2, * ' in your reaction: ') 4060 format (a8) 4050 format (' sorry ',a8,', that name is invalid, try again.',/) 4070 format (' Enter activity of: ',a8, * ' (enter 1.0 for H2O or CO2): ') 4000 format (/,' Calculate thermodynamic ', * 'properties for a reaction (y/n)? ') 4010 format (/,' How many phases or species in the ', * 'reaction? ') 5000 format (a1) 9999 end subroutine eqrxn (loop) c------------------------------------------------------------------------ implicit double precision (a-g,o-y),integer (h-m,z) parameter (l2=5) common/ cst5 /v(l2),tr,pr,r,ps/ cst9 /vmax(l2),vmin(l2),dv(l2) * / cst24 /ipot,jv(l2),iv(l2) c------------------------------------------------------------------------ c search for an equilibrium point c on the x-y coordinate frame. v(iv(3)) = vmin(iv(3)) do 10 i = 1, loop call newhld 10 v(iv(3)) = v(iv(3)) + dv(iv(3)) end subroutine vars (loop) c---------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k0=25,l2=5,k6=7,h5=7,h8=100) character*8 exname,title*72,y*1,n4name*14, * dname*40,cmpnt*5,vname*8,gname*78,xname*18 common/ cst6 /icomp,istct,iphct,icp/ cst79 /isoct * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / csta8 /title/ cst63 /delv(l2) * / csta2 /xname(k6),vname(l2),gname(l2)/ cst83 /ig(l2) * / cst5 /v(l2),tr,pr,r,ps/ cst9 /vmax(l2),vmin(l2),dv(5) * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn * / cst24 /ipot,jv(l2),iv(l2)/ cst36 /exname(h8) * / cst37 /iprct,ixct,iexyn,istbyn * / csta5 /dname(8),cmpnt(k0) ipot = 2 ier = 0 do 16 i = 1, 5 iv(i) = i vmin(i) = 0.0 vmax(i) = 0.0 16 dv(i) = 1.0 c components of saturated phase: iv(3)=3 if (ifyn.eq.0) then iv(3)=3 ipot = ipot + 1 end if write (*,*) ' Do you want to generate a graphics file (y/n)? ' read (*,1050) y io4 = 1 if (y.eq.'y'.or.y.eq.'Y') then write (*,1180) read (*,1190) n4name open (n4,file=n4name) io4 = 0 end if c select the x variable (iv(1)): 5006 write (*,2130) write (*,2140) (j,vname(iv(j)), j = 1, ipot) write (*,*) read (*,*,iostat=ier) ic call rerror (ier,*5006) ix = iv(1) iv(1) = iv(ic) iv(ic) = ix 5007 write (*,2150) vname(iv(1)) read (*,*,iostat=ier) vmin(iv(1)),vmax(iv(1)) call rerror (ier,*5007) c select the x variable (iv(2)): if (ipot.gt.2) then 5008 write (*,2110) write (*,2140) (j,vname(iv(j)), j = 2, ipot) write (*,*) read (*,*,iostat=ier) ic call rerror (ier,*5008) else ic=2 end if ix = iv(2) iv(2) = iv(ic) iv(ic) = ix 5009 write (*,2150) vname(iv(2)) read (*,*,iostat=ier) vmin(iv(2)),vmax(iv(2)) call rerror (ier,*5009) c define default variable increments: do 13 i = 1, 2 13 dv(iv(i)) = (vmax(iv(i)) - vmin(iv(i))) / 40.0 c--------------------------------------------------------------------- c specify sectioning variables (iv(3)): if (ipot.gt.2) then c check if multiple sections are desired: write (*,2160) read (*,1050) y c if (y.eq.'y'.or.y.eq.'Y') then c 5010 write (*,2150) vname(iv(3)) read (*,*,iostat=ier) vmin(iv(3)),vmax(iv(3)) call rerror (ier,*5010) 5011 write (*,1060) read (*,*,iostat=ier) ic call rerror (ier,*5011) dv(iv(3)) = (vmax(iv(3))-vmin(iv(3)))/(ic-1.) c specify remaining sectioning constraints: c only one section is to be calculated: else do 17 j=3,ipot 5012 write (*,2180) vname(iv(j)) read (*,*,iostat=ier) vmin(iv(j)) call rerror (ier,*5012) vmax(iv(j)) = vmin(iv(j)) 17 dv(iv(j)) = 1. end if c end if c--------------------------------------------------------------------- c compute interval range of variables do 117 i=1,3 117 delv(i) = vmax(i)-vmin(i) c calculate number of sections: loop = idint ((delv(iv(3)))/dv(iv(3))) + 1 c set convergence criteria for routine c univeq: call concrt goto (99),io4 c write a title card: write (*,1070) read (*,1170) title write (n4,*) 1 write (n4,*) 0, 0 write (n4,1160) title,' ',' ',' ' write (n4,*) ipot write (n4,1030) (vmax(iv(z)),vmin(iv(z)),z=1,ipot) write (n4,1040) (ig(iv(z)),gname(iv(z)),z=1,ipot) 1030 format (6(g11.5,1x)) 1040 format (i2,a78) 1050 format (a1) 1060 format (/,' Enter number of sections you want (min. 2): ') 1070 format (/,' Enter one-line calculation title:',/) 1160 format (a72,/,a1,/,a1,/,a1) 1170 format (a72) 1180 format (/,' Enter graphics file name, ' * ,'< 15 characters, left justified: ') 1190 format (a14) 2130 format (' Enter number identifying x-axis variable: ') 2140 format (10x,'(',i1,') ',a8) 2150 format (' Enter minimum and maximum values for ',a8,': ') 2110 format (' Enter number identifying y-axis variable: ') 2180 format (' Specify sectioning value for: ',a8) 2160 format (' Calculate sections as a function of a', * ' third variable (y/n)? ') 99 end subroutine fopen c-------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) character*14 n2name,yes*1 common/ cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 common/ cst41a /n2name c--------------------------------------------------------------------- c first the thermo data file 1 write (n8,1000) read (*,1030) n2name open (n2,file=n2name,iostat=ierr,status='old') if (ierr.ne.0) then c system could not find the file write (*,1010) n2name write (*,1040) read (*,1020) yes if (yes.ne.'y'.and.yes.ne.'Y') goto 999 goto 1 c try again end if return 999 write (*,1050) stop 1000 format (/,' Enter the thermo data file name, left justified,' * ,' < 15 characters',/,' (e.g. hp90ver.dat): ') 1010 format (/,' **error ver191** fopen file ',a14, * ' could not be opened',/) 1020 format (a1) 1030 format (a14) 1040 format (/,' Try again (y/n)? ') 1050 format (/,' o.k., then i am quitting also.',/) end subroutine newhld c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (l2=5,k6=7,h5=7) character*8 vname,gname*78,xname*18,y*1 common/ cst23 /a(10,10),b(11),ipvt(10),idv(10), * iophi,idphi,iiphi,iflg1 * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn * / cst6 /icomp,istct,iphct,icp * / cst5 /v(l2),tr,pr,r,ps * / cst9 /vmax(l2),vmin(l2),dv(l2) * / csta2 /xname(k6),vname(l2),gname(l2) * / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5 * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 c----------------------------------------------------------------------- c initialization: c ifuk = 0 10 write (n3,1160) write (n3,1170) vname(iv1),vname(iv2) c write potential variable sectioning c constraints: if (ipot.gt.2) write (n3,1180) vname(iv3), v(iv3) c set starting values for search v(iv1)=vmin(iv1) v(iv2)=vmin(iv2) c test stability along an edge of the c diagrams coordinate frame: 30 call search (ivi,ivd,div,ier) if (ier.eq.1) then write (*,1010) goto 20 end if call trace (ivd,ivi,div,igo) c this would permit continuation of the c search on failures from trace if iste c were not initialized in search (but here). c if (igo.eq.1.and.ifuk.eq.0) then c ifuk = 1 c goto 30 c else if (igo.eq.1.and.ifuk.eq.1) then c write (*,1030) c end if 20 write (*,1040) read (*,1000) y if (y.eq.'y'.or.y.eq.'Y') then call chptx goto 10 end if write (*,1020) read (*,1000) y if (y.ne.'y'.and.y.ne.'Y') goto 99 call change goto 10 1000 format (a1) 1010 format (/,' Equilibrium is not in specified', * ' coordinate frame.',/) 1020 format (' Modify data and', * ' recalculate the equilibrium (y/n)? ') 1030 format (' Oops, I cant find the equilibrium,', * ' maybe you should change PTX limits,',/, * ' i.e., answer yes to next question.',/) 1040 format (' Change PTX limits (y/n)?',/) 1160 format (/,'-------------------------------------------------' * ,'---------------',/) 1170 format (' The ',a8,'-',a8,' loci of the univariant field' * ,' follows:') 1180 format (' (subject to the constraint ',a8,'=',g10.4,')',/) 99 end subroutine search (ivi,ivd,div,ier) c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (l2=5,k1=1200) character*8 names common/ cst23 /a(10,10),b(11),ipvt(10),idv(10), * iophi,idphi,iiphi,iflg1 * / cst5 /v(l2),x(4)/ cst6 /icomp,istct,iphct,icp * / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5 * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst9 /vmax(l2),vmin(l2),dv(l2)/ cst8 /names(k1) c----------------------------------------------------------------------- c initialization ier=0 v(iv1) = vmin(iv1) v(iv2) = vmin(iv2) gst = grxn(ier) do 60 i = 1, 4 c set default dependent and independent c variables and increments for each edge c to be searched. goto (10,20,30,40),i c traverse 1. 10 ivi=iv2 ivd=iv1 ddv=dv(ivd) div=dv(ivi) v(ivi)=vmin(ivi) goto 80 c traverse 2. 20 ivi=iv1 ivd=iv2 ddv=dv(ivd) div=-dv(ivi) v(ivi)=vmax(ivi) goto 80 c traverse 3. 30 ivi=iv2 ivd=iv1 ddv=-dv(ivd) div=-dv(ivi) v(ivi)=vmax(ivi) goto 80 c traverse 4. 40 ivi=iv1 ivd=iv2 ddv=-dv(ivd) div=dv(ivi) v(ivi)=vmin(ivi) c begin search: 80 v(ivd)=v(ivd)+ddv c out of range?: goto (110,110,120,120),i 110 if (v(ivd).gt.vmax(ivd)) then v(ivd)=vmax(ivd) else if (i.eq.1) then if (v(ivd).gt.vmin(ivd)) goto 130 ddv=dabs(ddv)/2.d0 v(ivd)=vmin(ivd) iflg1=0 goto 80 end if goto 130 c 120 if (v(ivd).lt.vmin(ivd)) then v(ivd)=vmin(ivd) else if (i.eq.1) then if (v(ivd).lt.vmin(ivd)) goto 130 ddv=-dabs(ddv)/2.d0 v(ivd)= vmin(ivd) iflg1=0 goto 80 end if c calculate phase energies: 130 if (grxn(i)*gst.lt.0.0d0) goto 99 c check if search is in range: goto (140,140,150,150),i 140 if (v(ivd).ge.vmax(ivd)) goto 60 goto 80 150 if (v(ivd).le.vmin(ivd)) goto 60 goto 80 c next traverse. 60 continue c set this ier flag to indicate c the reaction wasn't found 90 ier=1 c done: 99 end subroutine trace (iovd,iovi,odiv,igo) c---------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (l2=5) common/ cst23 /a(10,10),b(11),ipvt(10),idv(10), * iophi,idphi,iiphi,iflg1 * / cst25 /vnu(11),idr(11),ivct * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst5 /v(l2),tr,pr,r,ps/ cst32 /ptx(450),logx,ipt2 c----------------------------------------------------------------------- ivi = iovi ivd = iovd igo = 0 50 call univeq (ivd,ier) c if univeq fails on a bounding edge c write error message and return: goto (9000,9000),ier c set the increment for the iv div=odiv c initialize counters ipt2=0 icter=0 c assign the 1st point call assptx c follow the equilibrium 60 call sfol1 (ivd,ivi,ier,div) c sfol1 returns ier=1 goto (70,70),ier ivi=iovi ivd=iovd goto 9999 70 call switch (div,ivi,ivd,jer) goto (75),jer icter=icter+1 if (icter.lt.4) goto 60 75 call warn (10, v(ivi), igo, 'TRACE') call outrxn ivi=iovi ivd=iovd c return on error goto 9999 c error in univeq: 9000 call warn (79, v(ivi), ird, 'TRACE') write (*,*) ' failed at P=',v(1),' T=',v(2),' XCO2 =',v(3) goto (9999), igo ivi = iovd ivd = iovi igo = 1 goto 50 9999 end subroutine sfol1 (ivd,ivi,ier,dv) c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (l2=5) common/ cst5 /v(l2),tr,pr,r,ps * / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5 * / cst9 /vmax(l2),vmin(l2),ddv(l2) * / cst6 /icomp,istct,iphct,icp/ cst32 /ptx(450),logx,ipt2 * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 c----------------------------------------------------------------------- c begin traverse: 10 v(ivi)=v(ivi)+dv c is search in range? if (v(ivi).gt.vmax(ivi)) then v(ivi)=vmax(ivi) else if (v(ivi).lt.vmin(ivi)) then v(ivi)=vmin(ivi) end if c solve for the equilibrium conditions: call univeq (ivd,ier) c on error return: c calling routine will switch variables. goto (9999,9999),ier c iflag=0: if (ipt2.gt.449) goto 9000 c dependent v in range? if c greater than the maximum value for v c or less than the minimum value for v c reset conditions, and c switch independent/dependent variables if (v(ivd).gt.vmax(ivd)) then v(ivd)=vmax(ivd) else if (v(ivd).lt.vmin(ivd)) then v(ivd)=vmin(ivd) else call assptx if ((v(ivi).eq.vmax(ivi)).or.(v(ivi).eq.vmin(ivi))) * goto 9000 goto 10 end if c solve for the equilibrium with c the switched variables: call univeq (ivi,ier) goto (9000,9000),ier c assign final value call assptx c output the traversed equilibrium: 9000 call outrxn ier=0 9999 end subroutine outrxn c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (l2=5,k1=1200) character*8 names common/ cst25 /vnu(11),idr(11),ivct/ cst8 /names(k1) * / cst6 /icomp,istct,iphct,icp * / cst32 /ptx(450),logx,ipt2/ cst5 /v(l2),q(4) * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 c----------------------------------------------------------------------- if (iphct.gt.4) goto 20 write (n3,1050) 1,(vnu(zl),names(zl),zl=1,iphct) goto 30 20 write (n3,1050) 1,(vnu(zl),names(zl),zl=1,4) write (n3,1060) (vnu(zl),names(zl),zl=5,iphct) 30 write (n3,*) write (n3,1000) (ptx(zi),zi=1,ipt2) write (n3,*) goto 10 10 goto (99), io4 if (ipt2.eq.0) goto 99 write (n4,*) ipt2,0,1,iphct,(i,i=1,iphct) write (n4,1020) (vnu(zl),zl=1,iphct) write (n4,1000) (ptx(zi),zi=1,ipt2) 1000 format (3(1x,g10.4,1x,g10.4,3x)) 1020 format (10(g9.3,1x)) 1050 format (/,' (',i3,')',4(1x,g9.3,1x,a6)) 1060 format (6x,4(1x,g9.3,1x,a6),/,6x,4(1x,g9.3,1x,a6)) 99 end function grxn (i) c-------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (h5=7) common/ cst5 /p,t,xco2,u1,u2,tr,pr,r,ps/ oxy /io2 * / cst201 /vuf(2),vus(h5),iffr,isr/ cst205 /act(11),idf(2) * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn/ cst11 /fh2o,fco2 * / cst6 /icomp,istct,iphct,icp/ cst25 /vnu(11),idr(11),ivct c--------------------------------------------------------------------- c compute free energy change of the rxn grxn = 0.d0 fh2o = 0.d0 fco2 = 0.d0 call updumy if (ifyn.eq.0) call cfluid (fo2,fs2) do 20 j = 1,iphct 20 grxn = grxn + vnu(j) * (gphase(j) + r * t * dlog(act(j))) grxn = grxn + vuf(1) * fh2o*r*t + vuf(2)*fco2*r*t if (io2.ne.0.and.ifyn.eq.0) grxn = grxn + vnu(io2)*r*t*fo2 end subroutine outlam (id) c--------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,m9=50,k9=50,m6=3,m7=9,m8=9) common/ cst203 /therdi(m8,m9),therlm(m7,m6,k9) * / cst100 /dlnfo2,elag,gz,gy,gx,ibuf,hu,hv,hw,hx * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst204 /ltyp(k1),lmda(k1),idis(k1) c----------------------------------------------------------------------- c routine only handles h&p and ubc if (ltyp(id).ne.0) then ilam = ltyp(id) jd = lmda(id) if (ilam.le.3) then do 10 j = 1, ilam 10 write (n2,1000) dsqrt (therlm(1,j,jd)), * dsqrt (therlm(2,j,jd)), * (therlm(k,j,jd), k = 3, 8), 0. else if (ilam.gt.9.and.ilam.lt.13) then do 20 j = 1, ilam - 9 20 write (n2,1010) therlm (1,j,jd), therlm(2,j,jd) end if end if if (idis(id).eq.0) goto 99 jd = idis(id) write (n2,1000) (therdi(j,jd), j = 1, m8) 1000 format (5(g13.7,1x)) 1010 format (f7.1,1x,g13.7,7(' 0.')) 99 end subroutine change c--------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (l2=5,k1=1200,k4=16,k5=12) parameter (m9=50,k9=50,m6=3,m7=9,m8=9,h5=7) character*8 names,list*28,y*1,n2name*14 common/ cst1 /thermo(k4,k1),vf(2,k1),vs(h5,k1),uf(2),us(h5) * / cst5 /p,t,xco2,u1,u2,tr,pr,r,ps/ cst8 /names(k1) * / cst44 /dlogt,dsqrtt,dexpt,dexpp/ cst206 /list(18) * / cst203 /therdi(m8,m9),therlm(m7,m6,k9) * / cst24 /ipot,jv(l2),iv(l2) * / cst205 /act(11),idf(2) * / cst100 /dlnfo2,elag,gz,gy,gx,ibuf,hu,hv,hw,hx * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn * / cst6 /icomp,istct,iphct,icp/ cst42 /ic(k5),idbase * / cst25 /vnu(11),idr(11),ivct / cst20 /comps(k1,k5) * / cst41a /n2name/ cst201 /vuf(2),vus(h5),iffr,isr * / cst204 /ltyp(k1),lmda(k1),idis(k1) c----------------------------------------------------------------------- write (*,1110) read (*,1060) y ier = 0 if (y.ne.'y'.and.y.ne.'Y') goto 20 write (*,1120) read (*,1060) y if (y.ne.'y'.and.y.ne.'Y') then 5013 write (*,1000) (z,names(z),z=1,iphct) read (*,*,iostat=ier) id call rerror (ier,*5013) write (*,1040) names(id) read (*,1050) names(id) call unver * (thermo(1,id),thermo(2,id),thermo(3,id),thermo(4,id), * thermo(5,id),thermo(6,id),thermo(7,id),thermo(8,id), * thermo(9,id),thermo(10,id),thermo(11,id), * thermo(12,id),thermo(13,id),thermo(14,id), * thermo(15,id),thermo(16,id)) c get lamda and dis codes if (ltyp(id).gt.3.and.ltyp(id).lt.10) goto 91 write (n2,1100) names(id),idbase,idis(id),ltyp(id),0 write (n2,1080) (comps(id,i),i=1,icomp) c add in activity correction gact = thermo(1,id) + tr * r * dlog (act(id)) sact = thermo(2,id) - r * dlog (act(id)) write (n2,1090) gact,sact,(thermo(i,id),i=3,k4) call outlam (id) goto 99 else id = 20 idis(id) = 0 ltyp(id) = 0 jdis = 0 klam = 0 write (*,1130) read (*,1140) names(id) do 15 i = 1, k4 15 thermo (i,id) = 0. do 11 i = 1, icomp 11 comps (id,i) = 0. do 12 j = 1, iphct do 13 i = 1, k4 13 thermo (i,id) = thermo(i,id) + vnu(j) * thermo(i,j) do 14 i = 1, icomp 14 comps(id,i) = comps(id,i) + vnu(j) * comps(j,i) if (idis(j).ne.0) then idis(id) = m9 jdis = jdis + 1 jd = idis(j) if (jdis.gt.1) goto 91 do 19 i = 1, 7 19 therdi(i,m9) = vnu(j) * therdi(i,jd) therdi(8,m9) = therdi(8,jd) therdi(9,m9) = therdi(9,jd) end if if (ltyp(j).ne.0) then ltyp(id) = ltyp(j) lmda(id) = k9 jd = lmda(j) klam = klam + 1 if (ltyp(j).eq.10) then if (klam.gt.3) goto 91 ltyp(id) = 9 + klam therlm(1,klam,k9) = therlm(1,1,jd) therlm(2,klam,k9) = vnu(j) * therlm(2,1,jd) else if (ltyp(j).lt.4) then if (klam.gt.1) goto 91 do 17 i = 1, ltyp(j) therlm(1,i,k9) = vnu(j) * therlm(1,i,jd) therlm(2,i,k9) = vnu(j) * therlm(2,i,jd) therlm(5,i,k9) = vnu(j) * therlm(5,i,jd) therlm(6,i,k9) = vnu(j) * therlm(6,i,jd) therlm(3,i,k9) = therlm(3,i,jd) therlm(4,i,k9) = therlm(4,i,jd) therlm(7,i,k9) = therlm(7,i,jd) 17 therlm(8,i,k9) = therlm(8,i,jd) else goto 91 end if end if 12 continue call unver * (thermo(1,id),thermo(2,id),thermo(3,id),thermo(4,id), * thermo(5,id),thermo(6,id),thermo(7,id),thermo(8,id), * thermo(9,id),thermo(10,id),thermo(11,id), * thermo(12,id),thermo(13,id),thermo(14,id), * thermo(15,id),thermo(16,id)) write (n2,1100) names(id),idbase,idis(id),ltyp(id),0 write (n2,1080) (comps(id,i),i=1,icomp) c add in activity correction gact = 0.d0 sact = 0.d0 do 16 i = 1, iphct gact = gact + vnu(i) * tr * r * dlog (act(i)) 16 sact = sact - vnu(i) * r * dlog (act(i)) write (n2,1090) thermo(1,id) + gact, * thermo(2,id) + sact,(thermo(i,id),i=3,k4) call outlam (id) goto 99 end if 20 write (*,1000) (z,names(z),z=1,iphct) read (*,*,iostat=ier) id call rerror (ier,*20) ichk = 0 write (*,1010) c 10 write (*,1020) (z,list(z),z=1,18) read (*,*,iostat=ier) kv call rerror (ier,*10) if (kv.eq.0) then if (ichk.eq.0) goto 30 c write entry to permanant file: write (*,1070) names(id),n2name read (*,1060) y c if (y.eq.'y'.or.y.eq.'Y') then write (n2,1100) names(id),idbase,ltyp(id),idis(id),0 write (n2,1080) (comps(id,i),i=1,icomp) c add in activity correction gact = thermo(1,id) + tr * r * dlog (act(id)) sact = thermo(2,id) - r * dlog (act(id)) write (n2,1090) gact,sact,(thermo(i,id),i=3,k4) call outlam (id) end if c call conver * (thermo(1,id),thermo(2,id),thermo(3,id),thermo(4,id), * thermo(5,id),thermo(6,id),thermo(7,id),thermo(8,id), * thermo(9,id),thermo(10,id),thermo(11,id), * thermo(12,id),thermo(13,id),thermo(14,id), * thermo(15,id),thermo(16,id)) goto 30 end if c if (ichk.eq.0) then call unver * (thermo(1,id),thermo(2,id),thermo(3,id),thermo(4,id), * thermo(5,id),thermo(6,id),thermo(7,id),thermo(8,id), * thermo(9,id),thermo(10,id),thermo(11,id), * thermo(12,id),thermo(13,id),thermo(14,id), * thermo(15,id),thermo(16,id)) write (*,1040) names(id) read (*,1050) names(id) ichk = 1 end if if (kv.eq.17) goto 50 if (kv.eq.18) goto 60 5014 write (*,1030) list(kv),names(id),thermo(kv,id) read (*,*,iostat=ier) thermo(kv,id) call rerror (ier,*5014) goto 10 50 write (*,1030) list(kv),names(id),act(id) read (*,*,iostat=ier) act(id) call rerror (ier,*50) goto 10 c 60 write (*,1030) list(kv),names(id),vnu(id) read (*,*,iostat=ier) vnu(id) call rerror (ier,*60) if (id.eq.idf(1)) then vuf(1) = vnu(id) else if (id.eq.idf(2)) then vuf(2) = vnu(id) end if goto 10 30 write (*,1150) read (*,1060) y if (y.eq.'y'.or.y.eq.'Y') goto 20 goto (99),ifyn write (*,1160) read (*,1060) y if (y.ne.'y'.and.y.ne.'Y') goto 99 call rfluid (1,ifug,n3) goto 99 91 call warn (9,t,ifyn,names(id)) 1000 format (/,' Select phase to modify', * ' or output:',9(/,6x,i2,') ',a8)) 1010 format (' Thermodynamic properties are calculated from the ', * 'reference state constants'/,' G, S, V, and an ', * 'activity coefficient together ', * 'with the isobaric heat',/,' capacity function:',//,6x, * 'cp(pr,t)=a+b*t+c/(t*t)+d*t*t+e/t**(1/2)+f/t+g/t**3',//, * ' and the volumetric function:',//,6x, * 'v(p,t)=v(pr,tr)+b2*(t-tr)+b3*(exp(-t/300)-exp(tr/300))', * '+b4*(p-pr)',/,12x,' +b5*(exp(-p/35000)-exp(-pr/35000))', * '+b6*(p-pr)**2+b7*(t-tr)**2') 1020 format (/,' Enter the number of the parameter to be modified:' * ,//,9(1x,i2,') ',a28,3x,i2,') ',a28,/),/, * ' enter zero when you are finished: ') 1030 format (/,' Old value for ',a28,' of ',a8,' was ',g15.8,/, * ' Enter new value: ') 1040 format (/,' Enter a name (<8 characters left justified) to', * ' distinguish the modified',/,' version of ',a8,'.', * ' WARNING: if you intend to store modified data',/, * ' in the data file this name must be unique.',/) 1050 format (a8) 1060 format (a1) 1070 format (/,' Store ',a8,' as an entry', * ' in data file ',a14,' (y/n)? ') 1080 format (12(f5.2,1x)) 1090 format (5(g13.7,1x)) 1100 format (a8,i2,i2,i2,i2) 1110 format (/,' Do you only want to output data (y/n)? ') 1120 format (/,' Output reaction properties (y/n)? ') 1130 format (/,' Enter an 8 character name for the reaction: ') 1140 format (a8) 1150 format (/,' Modify properties of another phase (y/n)? ') 1160 format (/,' Change fluid equation of state (y/n)? ') 99 end subroutine nentry c---------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k0=25,k4=16,k5=12) character*28 list,cmpnt*5,dname*40,name*8,y*1 common/ cst5 /p,t,xco2,u1,u2,tr,pr,r,ps * / csta5 /dname(8),cmpnt(k0)/ cst42 /ic(k5),idbase * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst206 /list(18)/ cst6 /icomp,istct,iphct,icp * / cst43 /therm(k4),comp(k0),atwt(k0),idh2o,idco2,ikind c----------------------------------------------------------------------- ier = 0 write (*,1000) dname(idbase),tr,pr write (*,1120) read (*,1100) y if (y.eq.'c') then fac = 4.184 else fac = 1.0 end if write (*,1130) read (*,1100) y isv = 0 if (y.eq.'v') isv = 1 50 write (*,1010) read (*,1020) name do 10 i = 1, icomp 5016 write (*,1030) cmpnt(i),name read (*,*,iostat=ier) comp(i) 10 call rerror (ier,*5016) do 20 i = 1, 2 5020 write (*,1040) list(i),name read (*,*,iostat=ier) therm(i) call rerror (ier,*5020) 20 therm(i) = fac * therm(i) 5017 write (*,1040) list(3),name read (*,*,iostat=ier) therm(3) call rerror (ier,*5017) write (*,1050) do 30 i = 4, 10 5018 write (*,1040) list(i),name read (*,*,iostat=ier) therm(i) call rerror (ier,*5018) 30 therm(i) = fac * therm(i) write (*,1060) do 40 i = 11, 16 5019 write (*,1040) list(i),name read (*,*,iostat=ier) therm(i) call rerror (ier,*5019) 40 if (isv.eq.1) therm(i) = therm(3) * therm(i) write (n2,1070) name,idbase,0,0,0 write (n2,1080) (comp(i),i=1,icomp) write (n2,1090) therm write (*,1110) read (*,1100) y if (y.eq.'y'.or.y.eq.'Y') goto 50 1000 format (/,' Your entry will be formatted for the ',a40,/, * ' data base with a t=',g13.6,'(k) p=',g13.6,'(bar) ',/, * ' reference state (all energy terms must be in joules).',/) 1010 format (' Enter name for your entry, <8 characters, left', * ' justified.',/,' WARNING: this name must not duplicate', * ' an entry already',/,' in the data file!') 1020 format (a8) 1030 format (' Enter number of moles of ',a5,' in ',a8,': ') 1040 format (' Enter ',a28,' for ',a8,': ') 1050 format (/,' Thermodynamic properties are calculated from', * ' the heat capacity function:',/,6x, * 'cp(pr,t)=a+b*t+c/(t*t)+d*t*t+e/t**(1/2)+f/t+g/t**3',/) 1060 format (/,' Thermodynamic properties are calculated from the', * ' volumetric function (j/bar):',/,6x, * 'v(p,t)=v(pr,tr)+b2*(t-tr)+b3*(exp(-t/300)-exp(tr/300))', * '+b4*(p-pr)',/,12x,' +b5*(exp(-p/35000)-exp(-pr/35000))', * '+b6*(p-pr)**2+b7*(t-tr)**2',/) 1070 format (a8,i2,i2,i2,i2) 1080 format (12(f5.2,1x)) 1090 format (5(g13.7,1x)) 1100 format (a1) 1110 format (/,' Make another entry (y/n)? ') 1120 format (/,' Enter a "c" to enter G, S, V, and a-g in calories: ') 1130 format (/,' Enter a "v" to scale b2-b7 by' * ,' standard state volume: ') end subroutine stoich c---------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,h5=7) character*8 names,y*1 common/ cst205 /act(11),idf(2)/ cst25 /vnu(11),idr(11),ivct * / cst5 /p,t,xco2,u1,u2,tr,pr,r,ps/ cst8 /names(k1) * / cst6 /icomp,istct,iphct,icp * / cst201 /vuf(2),vus(h5),iffr,isr c----------------------------------------------------------------------- ier = 0 20 write (*,1000) (z,names(z),vnu(z),z=1,iphct) write (*,*) read (*,*,iostat=ier) id call rerror (ier,*20) 5019 write (*,1030) names(id),vnu(id) read (*,*,iostat=ier) vnu(id) call rerror (ier,*5019) if (id.eq.idf(1)) then vuf(1) = vnu(id) else if (id.eq.idf(2)) then vuf(2) = vnu(id) end if write (*,1020) read (*,1010) y if (y.ne.'Y'.and.y.ne.'y') goto 99 goto 20 1000 format (/,' Enter number of phase to be modified:', * 9(/,6x,i2,') ',a8,' reaction coeff.=',f8.4)) 1010 format (a1) 1020 format (/,' Modify coefficient of another phase (y/n)? ') 1030 format (/,' Old coefficient for ',a8,' was ',f8.4, * ' enter new value: ') 99 end subroutine props (g,e,u,s,v,cp) c---------------------------------------------------------------------- c props calculates thermodynamic properties by finite difference c approximations from the gibbs function. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) common/ cst5 /p,t,xco2,u1,u2,tr,pr,r,ps c finite difference increments: dp = 0.1 dt = 0.01 p0 = p t0 = t g = grxn (i) p = p0 + dp gp= grxn (i) p = p0 t = t0 + dt gt= grxn(i) c entropy, volume, enthalpy: s = -(gt - g)/dt v = (gp -g)/dp e = s * t0 + g u = e - p0 * v c heat capacity, this should be centered c on t0, but it isn't: p = p0 t = t0 + 2.*dt gtt = grxn(i) s2 = -(gtt - gt) / dt e2 = s2 * (t0+dt) + gt cp = (e2 - e) / dt p = p0 t = t0 end