c Copyright (c) 1990 by James A. D. Connolly, Institute for Mineralogy c & Petrography, Swiss Federal Insitute of Technology, CH-8092 Zurich, c SWITZERLAND. All rights reserved. c Please do not distribute this source. c RLIB - a library of subprograms called by VERTEX, FRENDLY, ACTCOR c BUILD, ADIABAT, and ISOCHOR. subroutine rmodel (tname,istot) c--------------------------------------------------------------------- c rmodel - reads solution models from LUN n9. c--------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) parameter (k1=1000,j9=999,h9=18,m0=5) character*10 mname*8,tname,fname,names*8 common/ cst108 /wg(27,11),xmn(3,3),xmx(3,3),xnc(3,3),xfx(3), * iend(4,3,3),isub(27,9,2),imd(3),imsol(4,3,3), * ijk(3),insp(3,4),jsp(3),jsitp(3),iosp(3), * jst(3),jfix,idfc(3),ist(3),isp(3),isite, * ispct(3,4),jnsp(3,4),jsite,iterm,iord common/ cst107 /a0(5,5),acoef(5,5,m0),smult(5),ecoef(3,3,3), * nsite,nspm1(5),nterm(5,5), * nsub1(5,5,m0,4),nsub2(5,5,m0,4),nttyp(5,5,m0) common/ cst8 /names(k1)/ cst60 /iasmbl(j9),ipoint,imyn * / csta7 /fname(h9)/ cst18a /mname(4,3,3)/ cst79 /isoct * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 istot = 0 10 read (n9, 1000, end = 99) tname c if blank, assume comment. if (tname.eq.' ') goto 10 c read number of independent sites: read (n9,*,err=90) isite do 1 i = 1, 3 isp(i) = 1 1 ist(i) = 1 c read number of species on each site: read (n9,*,err=90) (isp(i),ist(i),i=1,isite) c total number of endmembers: istot = isp(1)*isp(2)*isp(3) c read endmember names: read (n9,1010) (((mname(i,j,k), i= 1,isp(1)), * j= 1,isp(2)), k= 1,isp(3)) c read endmember flags: read (n9,*,err=90) ((( iend(i,j,k), i= 1,isp(1)), * j= 1,isp(2)), k= 1,isp(3)) c read composition limits: do 20 i = 1, isite read (n9,*,err=90) (xmn(i,j),xmx(i,j),xnc(i,j), * j= 1, isp(i)-1),imd(i) c valid subdivision choice?: 20 if (imd(i).gt.10) call error (323,a0(1,1),i,'PRIME9') c read total number of terms, and order c of highest order term, in mixing model c ideal models: iterm = iord = 0. read (n9,*,err=90) iterm, iord if (iterm.eq.0) goto 40 c read parameters for margules c expression of the excess free c energy of the solution. do 30 h = 1, iterm read (n9,*,err=90) ((isub(h,i,j), j = 1, 2), i = 1, iord) 30 read (n9,*,err=90) (wg(h,k), k = 1, 11) c expansion for S(configurational) 40 read (n9,*,err=90) nsite if (nsite.ne.0) call smodel (tname) c permit fixed activity corrections: read (n9,*,err=90) jfix if (jfix.ne.0) read (n9,*,err=90) * (idfc(i), xfx(i),i = 1, jfix) goto 99 90 call error (20,wg(1,1),i,tname) 1000 format (a10) 1010 format (3(a8,1x)) 99 end subroutine zeroi (iarray,index,ivalue,n) integer iarray(n), n, index, ivalue do 10 i = 1, index 10 iarray(i) = ivalue end subroutine cmodel (im,idsol,tname,istot,jstot,ibuild,x1,x2) c--------------------------------------------------------------------- c cmodel - checks to see if solution models contain valid endmembers c and, if not, whether they can be recast as simpler models. c--------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) parameter (k1=1000,k2=800,k3=800,j9=999,h9=18) character*8 mname,tname*10,fname*10,names,missin(36),x1,x2 integer imiss common/ cst108 /wg(27,11),xmn(3,3),xmx(3,3),xnc(3,3),xfx(3), * iend(4,3,3),isub(27,9,2),imd(3),imsol(4,3,3), * ijk(3),insp(3,4),jsp(3),jsitp(3),iosp(3), * jst(3),jfix,idfc(3),ist(3),isp(3),isite, * ispct(3,4),jnsp(3,4),jsite,iterm,iord common/ cst18a /mname(4,3,3) common/ cst8 /names(k1)/ cst60 /iasmbl(j9),ipoint,imyn * / cst204 /ltyp(k1),lmda(k1),idis(k1)/ csta7 /fname(h9) * / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar * / cst79 /isoct jstot = 0 c if called by build (ibuild = 1) skip the c name check: if (ibuild.eq.1) goto 60 c check that solution name matches a c name requested in input from n1 do 50 i = 1, isoct if (tname.ne.fname(i)) goto 50 c got a match, goto 60: idsol = i im = im + 1 goto 60 c 50 continue c didn't find a match, read a new name: goto 99 c check that the endmembers match with data c from n2: c disable fixed activity corrections: 60 if (jfix.ne.0) then call warn (32,xfx(1),i,tname) goto 99 end if itl = 0 itd = 0 ith = 0 do 1 i = 1, 3 jst(i) = 1 do 2 j = 1, 4 ispct(i,j) = 0 do 3 k = 1, 3 3 imsol(j,i,k) = 0 2 insp(i,j) = 1 1 continue do 70 h = 1, ipoint do 80 i = 1, isp(1) do 90 j = 1, isp(2) do 100 k = 1, isp(3) ijk(1) = i ijk(2) = j ijk(3) = k c in build, check for forbidden choices: if (ibuild.eq.1) then if (mname(i,j,k).eq.x1.or. * mname(i,j,k).eq.x2) then write (*,1010) tname, x1, x2 jstot = 0 goto 99 end if end if if (names(h).ne.mname(i,j,k)) goto 100 c check for funny transitions, kill c pgm execution if they occur: if (idis(h).ne.0) then itd = itd + 1 if (itd.gt.1) call error (61,xfx(1),i,tname) end if if (ltyp(h).eq.10) then ith = ith + 1 if (ith.gt.3) call error (61,xfx(1),i,tname) else if (ltyp(h).gt.0) then itl = itl + 1 if (itl.gt.1) call error (61,xfx(1),i,tname) end if c got a valid match, goto 70 or 110 jstot = jstot + 1 imsol(i,j,k) = h c do 101 i1 = 1, isite do 103 j1 = 1, isp(i1) 103 if (ijk(i1).eq.j1) ispct(i1,j1)=ispct(i1,j1)+1 101 continue if (iend(i,j,k).eq.0) then ikp(h) = idsol else if (iend(i,j,k).eq.2) then c set kill flag: ikp(h) = -1 end if c found all endmembers: if (jstot.eq.istot) goto 99 goto 70 100 continue 90 continue 80 continue 70 continue c missing endmember warnings: if (jstot.eq.0) then call warn (25,xfx(1),i,tname) else if (jstot.eq.1) then jstot = 0 call warn (34,xfx(1),i,tname) else imiss = 0 do 10 i = 1, isp(1) do 20 j = 1, isp(2) do 30 k = 1, isp(3) if (imsol(i,j,k).ne.0) goto 30 imiss = imiss + 1 missin(imiss) = mname(i,j,k) 30 continue 20 continue 10 continue write (*,1000) tname,(missin(i),i=1,imiss) end if 1000 format (' **warning ver114** the following endmembers', * ' are missing for ',a10,//,4(8(2x,a8),/)) 1010 format (' **warning ver113** ',a10,' is not a valid model', * ' because component ',a5,' or ',a5, * ' is constrained',/) 99 end subroutine emodel (tname,istot) c--------------------------------------------------------------------- c emodel - attempts to recast solution models as simpler models. c--------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) parameter (k1=1000,k2=800,k3=800) character*8 tname*10,names integer lmn(3) common/ cst108 /wg(27,11),xmn(3,3),xmx(3,3),xnc(3,3),xfx(3), * iend(4,3,3),isub(27,9,2),imd(3),imsol(4,3,3), * ijk(3),insp(3,4),jsp(3),jsitp(3),iosp(3), * jst(3),jfix,idfc(3),ist(3),isp(3),isite, * ispct(3,4),jnsp(3,4),jsite,iterm,iord common/ cst8 /names(k1) * / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar do 1 i = 1, 3 iosp(i) = 0 do 2 j = 1, 4 2 insp(i,j) = 0 1 continue c look for sites where there c are enough endmembers, assuming c that the count of a good species c should be if (isite.gt.1) then do 410 i = 1, isp(1) ijk(1) = i do 420 j = 1, isp(2) ijk(2) = j do 430 k = 1, isp(3) ijk(3) = k if (imsol(i,j,k).ne.0) then c check if the endmember can c make a binary: do 440 i1 = 1, isp(1) lmn(1) = i1 do 450 j1 = 1, isp(2) lmn(2) = j1 do 460 kk = 1, isp(3) lmn(3) = kk if (imsol(i1,j1,kk).ne.0) then do 470 i2 = 1, isite if (ijk(i2).ne.lmn(i2)) then do 480 i3 = 1, isite if (i3.eq.i2) goto 480 if (ijk(i3).ne.lmn(i3)) goto 460 c species ijk(i2) and lmn(i2) are c valid on site i2: insp(i2,ijk(i2)) = 1 insp(i2,lmn(i2)) = 1 480 continue end if 470 continue end if 460 continue 450 continue 440 continue end if 430 continue 420 continue 410 continue else do 610 j = 1, isp(1) 610 if (ispct(1,j).gt.0) insp(1,j) = 1 end if c count the number of endmembers c left: 600 do 601 i = 1, 3 jsp(i) = 0 do 602 j = 1, 4 602 jnsp(i,j) = 1 601 continue c and make a set of valid species do 490 i = 1, isite do 500 j = 1, isp(i) if (insp(i,j).eq.0) goto 500 jsp(i) = jsp(i) + 1 jnsp(i,j) = jsp(i) 500 continue 490 continue jstot = 0 do 510 i = 1, isp(1) do 520 j = 1, isp(2) do 530 k = 1, isp(3) if (imsol(i,j,k).eq.0) goto 530 ijk(1) = i ijk(2) = j ijk(3) = k do 535 l = 1, isite if (jsp(l).ne.0.and.insp(l,ijk(l)).eq.0) then c endmember is not valid: if (iend(i,j,k).eq.1) ikp(imsol(i,j,k)) = 0 imsol(i,j,k) = 0 end if 535 continue jstot = jstot + 1 530 continue 520 continue 510 continue kstot = 1 do 25 i = 1, 3 if (jsp(i).eq.0) jsp(i) = 1 kstot = kstot * jsp(i) 25 jsitp(i) = i if (kstot.lt.2) then call warn (34,wg(1,1),i,tname) istot = 0 end if c if jsp(1)*jsp(2)*jsp(3) > jstot c it's possible to make several c different solutions, eliminate c a species on first site possible c and try again: if (kstot.gt.jstot) then call warn (33,wg(1,1),i,tname) imin = 100 do 540 i = 1, isite do 550 j = 1, isp(i) if (jsp(i).gt.1.and.insp(i,j).ne.0.and. * ispct(i,j).lt.imin) then imin = ispct(i,j) ik = i jk = j end if 550 continue 540 continue insp(ik,jk) = 0 goto 600 end if c make a pointer for old site c index to new: jsite = 0 jmax = isite do 30 i = 1, isite if (jsp(i).eq.1) then jsitp(i) = jmax jst(jmax) = 1 jmax = jmax - 1 else jsite = jsite + 1 jsitp(i) = jsite jst(jsite) = ist(i) end if 30 continue c locate non-zero species id's do 80 i = 1, isp(1) do 90 j = 1, isp(2) do 100 k = 1, isp(3) if (imsol(i,j,k).eq.0) goto 100 ijk(1) = i ijk(2) = j ijk(3) = k do 110 l = 1, 3 110 if (isp(l).gt.1.and.jsp(l).eq.1 * .and.insp(l,ijk(l)).eq.0) iosp(l) = ijk(l) 100 continue 90 continue 80 continue c if this works, then the next loop c could be lot simpler c relocate endmembers: do 40 i = 1, isp(1) if (i.ne.iosp(1).and.insp(1,i).eq.0) goto 40 ijk(jsitp(1)) = jnsp(1,i) do 50 j = 1, isp(2) if (j.ne.iosp(2).and.jsite.gt.1.and.insp(2,j).eq.0) * goto 50 ijk(jsitp(2)) = jnsp(2,j) do 60 k = 1, isp(3) if (k.ne.iosp(3).and.jsite.gt.2.and.insp(3,k).eq.0) * goto 60 ijk(jsitp(3)) = jnsp(3,k) if (imsol(i,j,k).eq.0) goto 60 imsol(ijk(1),ijk(2),ijk(3)) = imsol(i,j,k) 60 continue 50 continue 40 continue c relocate compositional limits: do 620 i = 1, isite if (jsp(i).gt.1) then do 630 j = 1, jsp(i) - 1 xmx(jsitp(i),j) = xmx(i,j) xmn(jsitp(i),j) = xmn(i,j) 630 xnc(jsitp(i),j) = xnc(i,j) if (jsp(i).eq.2.and.isp(i).eq.3.and.imd(i).gt.0) then imd(jsitp(i)) = 0 if (xnc(jsitp(i),j).eq.0.or.xnc(jsitp(i),j).gt. * xmx(jsitp(i),j)-xmn(jsitp(i),j)) xnc(jsitp(i),j) = * (xmx(jsitp(i),j)-xmn(jsitp(i),j))/10. write (*,1140) tname,jsitp(i),i, * imd(jsitp(i)),xmx(jsitp(i),1), * xmn(jsitp(i),1),xnc(jsitp(i),1) else imd(jsitp(i)) = imd(i) end if end if 620 continue c redefine the other counters: istot = kstot do 71 i = 1, 3 71 isp(i) = 1 do 70 i = 1, isite ist(i) = jst(jsitp(i)) 70 isp(i) = jsp(jsitp(i)) isite = jsite c write warning: write (*,2000) tname write (*,2010) (((i,j,k,names(imsol(i,j,k)), * i = 1, isp(1)),j = 1, isp(2)), k = 1, isp(3)) write (*,*) 1140 format (/,' **warning ver117** the subdivision scheme for ' * ,a10,' will be interpreted',/, * ' as a binary model on site ',i1,' (originally site ',i1 * ,') with:',/,t12,'imod=', * 1x,i1,/,t12,'xmax= ',f7.3,/,t12,'xmin= ',f7.3,/,t12, * 'xinc= ',f7.3,/,' To change this, or if an' * ' error (probably from routine CART) follows,',/, * ' modify the solution model', * ' (documentation sect. 4.1)',/) 2000 format (/,' **warning ver501** ',a10, * ' will be recast with endmembers:',/) 2010 format (10x,5(i1,i1,i1,' - ',a8,' ')) end subroutine nmodel c---------------------------------------------------------------------- c nmodel - nmodel reassigns solution model parameters to be c consistent with the endmembers chosen by subroutine emodel. c---------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) parameter (k1=1000,k2=800,k3=800,m0=5) character*8 names common/ cst108 /wg(27,11),xmn(3,3),xmx(3,3),xnc(3,3),xfx(3), * iend(4,3,3),isub(27,9,2),imd(3),imsol(4,3,3), * ijk(3),insp(3,4),jsp(3),jsitp(3),iosp(3), * jst(3),jfix,idfc(3),ist(3),isp(3),isite, * ispct(3,4),jnsp(3,4),jsite,iterm,iord common/ cst8 /names(k1) * / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar common/ cst107 /a0(5,5),acoef(5,5,m0),smult(5),ecoef(3,3,3), * nsite,nspm1(5),nterm(5,5), * nsub1(5,5,m0,4),nsub2(5,5,m0,4),nttyp(5,5,m0) c next eliminate unused terms c in margules and configurational c entropy expressions. if (iterm.ne.0) then jterm = 0 do 10 h = 1, iterm do 20 k = 1, iord if (isub(h,k,1).eq.0) then goto 20 else if (insp(isub(h,k,1),isub(h,k,2)).eq.0) then goto 10 end if 20 continue c term contains no invalid species c count it: jterm = jterm + 1 do 30 i = 1, iord if (isub(h,i,1).eq.0) then isub(jterm,i,1) = 0 isub(jterm,i,2) = 0 else isub(jterm,i,1) = jsitp(isub(h,i,1)) isub(jterm,i,2) = jnsp(isub(h,i,1),isub(h,i,2)) end if 30 continue do 40 j = 1, 11 40 wg(jterm,j) = wg(h,j) 10 continue iterm = jterm end if c now entropy terms: if (nsite.ne.0) then lsite = 0 do 50 i = 1, nsite lsp = 0 do 60 j = 1, nspm1(i) l2 = 0 l4 = 0 if (a0(i,j).ne.0) then c the z function contains a c constant if (lsp.eq.0) lsite = lsite + 1 lsp = lsp + 1 l2 = 1 l4 = 1 a0(lsite,lsp) = a0(i,j) smult(lsite) = smult(i) end if c look if it contains any c valid terms lterm = 0 l3 = 0 do 70 k = 1, nterm(i,j) lin = 0 do 80 l = 1, nttyp(i,j,k) n1 = nsub1(i,j,k,l) n2 = nsub2(i,j,k,l) if (insp(n1,n2).eq.0.and.n2.ne.iosp(n1)) goto 70 if (n2.eq.iosp(n1)) lin = lin + 1 80 continue c the term contributes to the c site fraction if (lsp.eq.0) then lsite = lsite + 1 lsp = 1 l2 = 1 l3 = 1 a0(lsite,lsp) = 0.0 smult(lsite) = smult(i) else if (l3.eq.0) then if (l4.eq.0) then lsp = lsp + 1 a0(lsite,lsp) = 0.0 smult(lsite) = smult(i) end if l3 = 1 end if if (lin.eq.nttyp(i,j,k)) then c the term contains only c constants: a0(lsite,lsp) = a0(lsite,lsp) + acoef(i,j,k) else c the term contains at least one X. lterm = lterm + 1 acoef(lsite,lsp,lterm) = acoef(i,j,k) lin = 0 do 90 l = 1, nttyp(i,j,k) n1 = nsub1(i,j,k,l) n2 = nsub2(i,j,k,l) if (n2.ne.iosp(n1)) then lin = lin + 1 nsub1(lsite,lsp,lterm,lin) = jsitp(n1) nsub2(lsite,lsp,lterm,lin) = jnsp(n1,n2) end if 90 continue nttyp(lsite,lsp,lterm) = lin end if 70 continue if (lterm.gt.0.or.l2.gt.0) nterm(lsite,lsp) = lterm 60 continue if (lsp.gt.0) nspm1(lsite) = lsp 50 continue nsite = lsite end if end subroutine smodel (tname) c----------------------------------------------------------------- c smodel reads models for configurational entropy c of solutions with dependent site mixing or disordered c endmembers, see documentation section 1.3.1, equation (8b) c----------------------------------------------------------------- implicit double precision (a-g,o-z),integer (h-n) parameter (m0=5) character*10 tname common/ cst107 /a0(5,5),acoef(5,5,m0),smult(5),ecoef(3,3,3), * nsite,nspm1(5),nterm(5,5), * nsub1(5,5,m0,4),nsub2(5,5,m0,4),nttyp(5,5,m0) common/ cst108 /wg(27,11),xmn(3,3),xmx(3,3),xnc(3,3),xfx(3), * iend(4,3,3),isub(27,9,2),imd(3),imsol(4,3,3), * ijk(3),insp(3,4),jsp(3),jsitp(3),iosp(3), * jst(3),jfix,idfc(3),ist(3),isp(3),isite, * ispct(3,4),jnsp(3,4),jsite,iterm,iord * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 c----------------------------------------------------------------- if (nsite.gt.5) call error (31,a0(1,1),isite,tname) c if (nsite.lt.isite) call error (30,a0(1,1),isite,tname) c for each site do 10 i = 1, nsite c read # of species, and site c multiplicty. read (n9,*,err=90) nsp,smult(i) nspm1(i) = nsp - 1 c for each species, read c function to define the c site fraction of the species: do 20 j = 1, nspm1(i) c read # of terms in the c site fraction function and a0. read (n9,*,err=90) nterm(i,j), a0(i,j) if (nterm(i,j).gt.m0) call error (33,a0(1,1),m0,tname) c for each term: do 30 k = 1, nterm(i,j) c read term type: read (n9,*,err=90) nttyp(i,j,k) if (nttyp(i,j,k).lt.5) then read (n9,*,err=90) acoef(i,j,k), * (nsub1(i,j,k,l), * nsub2(i,j,k,l), * l = 1, nttyp(i,j,k) ) else call error (29,a0(1,1),nttyp(i,j,k),tname) end if 30 continue 20 continue 10 continue return 90 call error (20,a0(1,1),i,tname) end subroutine snorm (tname) c---------------------------------------------------------------------- implicit double precision (a-g,o-z),integer (h-n) parameter (m0=5) character*10 tname double precision x(3,3) common/ cst107 /a0(5,5),acoef(5,5,m0),smult(5),ecoef(3,3,3), * nsite,nspm1(5),nterm(5,5), * nsub1(5,5,m0,4),nsub2(5,5,m0,4),nttyp(5,5,m0) * / cst5 /p,t,xc,u1,u2,tr,pr,r,ps common/ cst108 /wg(27,11),xmn(3,3),xmx(3,3),xnc(3,3),xfx(3), * iend(4,3,3),isub(27,9,2),imd(3),imsol(4,3,3), * ijk(3),insp(3,4),jsp(3),jsitp(3),iosp(3), * jst(3),jfix,idfc(3),ist(3),isp(3),isite, * ispct(3,4),jnsp(3,4),jsite,iterm,iord c get normalization constants: do 40 i = 1, isp(1) c set x's to endmember values: x(1,i) = 1. do 50 i1 = 1, isp(1) 50 if (i1.ne.i) x(1,i1) = 0. do 60 j = 1, isp(2) x(2,j) = 1. do 70 j1 = 1, isp(2) 70 if (j1.ne.j) x(2,j1) = 0. do 80 k = 1, isp(3) x(3,k) = 1. do 90 k1 = 1, isp(3) 90 if (k1.ne.k) x(3,k1) = 0. c x's now set, evaluate ln W: call omega (x,dlnw,tname) ecoef(i,j,k) = dlnw 80 continue 60 continue 40 continue write (*,2000) tname write (*,2010) (((i,j,k,r*ecoef(i,j,k), * i = 1, isp(1)),j = 1, isp(2)), k = 1, isp(3)) write (*,*) 2000 format (/,' Endmember configurational entropies (', * 'doc. eq. 8.2) for ',a10,' are:',/) 2010 format (5(2x,i1,i1,i1,' - ',f8.5,2x)) end subroutine omega (x,dlnw,tname) c---------------------------------------------------------------------- c subroutine to evaluate the log of the number of site configurations c per mole of a solution (dlnw) with composition x. c---------------------------------------------------------------------- implicit double precision (a-g,o-z),integer (h-n) parameter (m0=5) character*10 tname double precision x(3,3) common/ cst107 /a0(5,5),acoef(5,5,m0),smult(5),ecoef(3,3,3), * nsite,nspm1(5),nterm(5,5), * nsub1(5,5,m0,4),nsub2(5,5,m0,4),nttyp(5,5,m0) dlnw = 0.0 c for each site do 10 i = 1, nsite zt = 0. c get site fractions do 20 j = 1, nspm1(i) z = a0(i,j) if (nterm(i,j).ne.0) then c for each term: do 30 k = 1, nterm(i,j) zad = acoef(i,j,k) do 40 l = 1, nttyp(i,j,k) 40 zad = zad * x(nsub1(i,j,k,l),nsub2(i,j,k,l)) z = z + zad 30 continue end if call zchk (z,tname) zt = zt + z if (z.gt.0.0.and.z.lt.1.0) then dlnw = dlnw - smult(i) * z * dlog (z) end if 20 continue z = 1. - zt call zchk (z,tname) if (z.gt.0.0.and.z.lt.1.0) * dlnw = dlnw - smult(i) * z * dlog (z) 10 continue end subroutine subdiv (tname) c--------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n) character*10 tname parameter (k1=1000) common/ cst86 /xy(2,k1),y(k1,3,2),ntot,npairs common/ cst108 /wg(27,11),xmn(3,3),xmx(3,3),xnc(3,3),xfx(3), * iend(4,3,3),isub(27,9,2),imd(3),imsol(4,3,3), * ijk(3),insp(3,4),jsp(3),jsitp(3),iosp(3), * jst(3),jfix,idfc(3),ist(3),isp(3),isite, * ispct(3,4),jnsp(3,4),jsite,iterm,iord c--------------------------------------------------------------------- c do the first site: call cart (xmn(1,1),xmx(1,1),xnc(1,1), * xmn(1,2),xmx(1,2),xnc(1,2),isp(1),imd(1),tname) c ntot = npairs np1 = npairs do 10 h = 1, np1 do 20 i = 1, isp(1) - 1 20 y(h,1,i) = xy(i,h) 10 continue if (isite.eq.1) goto 99 c do the second site: call cart (xmn(2,1),xmx(2,1),xnc(2,1), * xmn(2,2),xmx(2,2),xnc(2,2),isp(2),imd(2),tname) c there will be a total of c (npairs-1)*np1 compositions, c copy the first site 2 distribution c into the first np1 compositions: do 30 h = 1, np1 do 40 i = 1, isp(2) - 1 40 y(h,2,i) = xy(i,1) 30 continue do 50 h = 2, npairs c for each site 2 composition, c duplicate the range of site 1 c compositions. do 60 i = 1, np1 ntot = ntot + 1 do 70 j = 1, isp(1) - 1 70 y(ntot,1,j) = y(i,1,j) c put in the new site 2 compositions: do 80 j = 1, isp(2) - 1 80 y(ntot,2,j) = xy(j,h) 60 continue 50 continue c do the third site: if (isite.eq.2) goto 99 np1 = (npairs-1) * np1 call cart (xmn(3,1),xmx(3,1),xnc(3,1), * xmn(3,2),xmx(3,2),xnc(3,2),isp(3),imd(3),tname) c copy the first site 3 distribution c into the first np1*np2 compositions: do 90 h = 1, ntot do 100 i = 1, isp(3) - 1 100 y(h,3,i) = xy(i,1) 90 continue c for each site 3 composition, c duplicate the range of site 1 and c site 2 compositions. do 110 h = 2, npairs do 120 i = 1, np1 ntot = ntot + 1 do 130 j = 1, isp(1) - 1 130 y(ntot,1,j) = y(i,1,j) do 140 j = 1, isp(2) - 1 140 y(ntot,2,j) = y(i,1,j) do 150 j = 1, isp(3) - 1 150 y(ntot,3,j) = xy(j,h) 120 continue 110 continue 99 end subroutine cart (xmin1,xmax1,xinc1,xmin2,xmax2,xinc2, * iscp,imod,sname) c--------------------------------------------------------------------- c xmax2 |\ cartesian subdivision onto 1- & 2- c | \ dimensional simplexes. iscp is number of c | \ vertices, in 1-dimensional case c xinc2 | \ ordinates of vertices are xmin1 & xmax1, c | \ in 2-dimensional case coordinates c | \ are (xmin1, xmin2), (xmax1, xmin2) & c xmin2 | \ (xmin1, xmax2). c --------------- c xmin1 xinc1 xmax1 c--------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) parameter (k1=1000) real vco1(673), vco2(673) integer ibary(4) character sname*10 save vco1, vco2, ibary common/ cst86 /xy(2,k1),yy(k1,3,2),ntot,npairs c ibary(j) is the number of stars after c jth barycentric subdivision. data ibary/7,25,121,673/ c (vco1(i), vco2(i)) are the simplicial c coordinates of the stars generated by c barycentric subdivsion. data (vco1(zi),zi = 1,3)/ 0., 1., 0./ data (vco2(zi),zi = 1,3)/ 1., 0., 0./ data (vco1(zi),zi = 4,121)/.500,.500,.0 ,.333, c 4 after subdivision: 1 *.750,.611,.250,.417,.278,.0 ,.167,.111,.0 ,.167,.111,.250, *.167,.278,.750,.667,.417,.611, c 22 after subdivision: 2 *.875,.787,.625,.681,.620,.556,.509,.514,.454,.472,.537,.806, *.639,.759,.375,.343,.125,.264,.176,.139,.148,.222,.259,.375, *.306,.343,.458,.389,.347,.398,.0 ,.037,.0 ,.056,.037,.056, *.093,.139,.204,.250,.222,.204,.083,.056,.139,.093,.0 ,.037, *.0 ,.056,.037,.056,.093,.139,.204,.250,.222,.204,.083,.056, *.139,.093,.125,.176,.375,.264,.343,.389,.398,.347,.343,.250, *.306,.259,.083,.139,.222,.148,.875,.787,.625,.681,.620,.458, *.556,.509,.375,.514,.454,.500,.472,.537,.833,.806,.639,.759/ c 118 after subdivision: 3 data (vco1(zi),zi = 122,337)/ *.938,.887,.813,.831,.804,.769,.739,.734,.693,.699,.735,.894, *.796,.864,.688,.665,.563,.623,.582,.560,.559,.588,.596,.646, *.616,.637,.715,.685,.650,.684,.489,.484,.461,.463,.480,.512, *.545,.583,.560,.559,.528,.505,.532,.522,.415,.414,.387,.394, *.420,.463,.512,.562,.532,.526,.465,.435,.484,.461,.457,.519, *.568,.602,.614,.588,.596,.542,.574,.540,.403,.435,.505,.448, *.864,.796,.753,.653,.713,.688,.625,.699,.670,.708,.685,.725, *.903,.880,.782,.855,.438,.406,.313,.359,.323,.296,.285,.303, *.295,.310,.336,.421,.366,.410,.188,.184,.063,.150,.100,.088, *.105,.157,.198,.271,.227,.239,.257,.213,.220,.230,.077,.116, *.133,.157,.179,.185,.216,.208,.213,.188,.069,.074,.144,.096, *.225,.255,.281,.296,.299,.282,.281,.250,.269,.253,.194,.213, *.241,.216,.354,.350,.396,.359,.378,.380,.369,.345,.323,.292, *.310,.309,.319,.338,.324,.327,.479,.452,.437,.428,.424,.382, *.407,.387,.312,.373,.341,.333,.338,.355,.444,.449,.394,.429, *.0 ,.012,.0 ,.019,.012,.019,.031,.046,.068,.074,.068,.019, *.046,.031,.0 ,.012,.0 ,.019,.012,.019,.031,.046,.068,.083, *.074,.068,.028,.019,.046,.031,.059,.088,.114,.130,.133,.116/ data (vco1(zi),zi = 338,553)/ *.114,.083,.102,.086,.028,.046,.074,.049,.207,.227,.262,.269, *.253,.213,.179,.125,.157,.151,.153,.185,.171,.170,.292,.262, *.208,.227,.207,.185,.170,.171,.151,.167,.157,.179,.278,.269, *.213,.253,.042,.059,.125,.088,.114,.153,.130,.133,.125,.116, *.114,.083,.102,.086,.028,.046,.074,.049,.0 ,.012,.0 ,.019, *.012,.019,.031,.046,.068,.074,.068,.019,.046,.031,.0 ,.012, *.0 ,.019,.012,.019,.031,.046,.068,.083,.074,.068,.028,.019, *.046,.031,.059,.088,.114,.130,.133,.116,.114,.083,.102,.086, *.028,.046,.074,.049,.207,.227,.262,.269,.253,.213,.179,.125, *.157,.151,.153,.185,.171,.170,.292,.262,.208,.227,.207,.185, *.170,.171,.151,.167,.157,.179,.278,.269,.213,.253,.042,.059, *.125,.088,.114,.153,.130,.133,.125,.116,.114,.083,.102,.086, *.028,.046,.074,.049,.063,.100,.188,.150,.184,.213,.230,.220, *.239,.227,.198,.088,.157,.105,.313,.323,.438,.359,.406,.421, *.410,.366,.336,.271,.310,.295,.257,.296,.303,.285,.452,.428, *.424,.407,.387,.373,.341,.333,.338,.355,.444,.449,.394,.429, *.378,.359,.350,.338,.327,.324,.309,.312,.310,.323,.382,.380, *.345,.369,.292,.281,.208,.255,.225,.213,.216,.241,.253,.292/ data (vco1(zi),zi = 554,673)/ *.269,.281,.319,.296,.282,.299,.042,.077,.125,.116,.133,.194, *.157,.179,.250,.185,.216,.208,.213,.188,.069,.074,.144,.096, *.938,.887,.813,.831,.804,.769,.739,.734,.693,.699,.735,.894, *.796,.864,.688,.665,.563,.623,.582,.560,.559,.588,.596,.646, *.616,.637,.715,.685,.650,.684,.479,.489,.437,.484,.461,.463, *.480,.512,.545,.583,.560,.559,.528,.505,.532,.522,.396,.415, *.354,.414,.387,.394,.420,.463,.512,.562,.532,.526,.465,.435, *.484,.461,.417,.457,.583,.519,.568,.602,.614,.588,.596,.542, *.574,.540,.403,.435,.505,.448,.917,.864,.750,.796,.753,.653, *.713,.688,.625,.699,.670,.708,.685,.725,.903,.880,.782,.855/ c 670 after subdivision: 4 data (vco2(zi),zi = 4,121)/.500,.0 ,.500,.333, c 4 after subdivision: 1 *.250,.278,.750,.417,.611,.750,.667,.611,.250,.417,.278,.0 , *.167,.111,.0 ,.167,.167,.111, c 22 after subdivision: 2 *.125,.176,.375,.264,.343,.389,.398,.347,.343,.306,.259,.139, *.222,.148,.625,.620,.875,.681,.787,.806,.759,.639,.537,.375, *.472,.454,.458,.556,.514,.509,.875,.787,.625,.681,.620,.556, *.509,.514,.454,.500,.472,.537,.833,.806,.639,.759,.375,.343, *.125,.264,.176,.139,.148,.222,.259,.375,.306,.343,.458,.389, *.347,.398,.0 ,.037,.0 ,.056,.037,.056,.093,.139,.204,.250, *.222,.204,.083,.056,.139,.093,.0 ,.037,.0 ,.056,.037,.083, *.056,.093,.250,.139,.204,.250,.222,.204,.083,.056,.139,.093/ c 118 after subdivision: 3 data (vco2(zi),zi = 122,337)/ *.063,.100,.188,.150,.184,.213,.230,.220,.239,.227,.198,.088, *.157,.105,.313,.323,.438,.359,.406,.421,.410,.366,.336,.271, *.310,.295,.257,.296,.303,.285,.452,.428,.424,.407,.387,.373, *.341,.333,.338,.355,.444,.449,.394,.429,.378,.359,.350,.338, *.327,.324,.309,.312,.310,.323,.382,.380,.345,.369,.281,.255, *.225,.213,.216,.241,.253,.292,.269,.281,.319,.296,.282,.299, *.077,.116,.133,.194,.157,.179,.250,.185,.216,.208,.213,.188, *.069,.074,.144,.096,.563,.582,.688,.623,.665,.685,.684,.650, *.637,.616,.596,.560,.588,.559,.813,.804,.938,.831,.887,.894, *.864,.796,.735,.646,.699,.693,.715,.769,.734,.739,.864,.796, *.753,.713,.688,.699,.670,.708,.685,.725,.903,.880,.782,.855, *.568,.519,.457,.435,.448,.505,.540,.625,.574,.596,.653,.602, *.588,.614,.354,.387,.396,.414,.415,.435,.461,.484,.526,.542, *.532,.512,.403,.394,.463,.420,.479,.489,.437,.484,.461,.465, *.463,.480,.562,.512,.545,.583,.560,.559,.528,.505,.532,.522, *.938,.887,.813,.831,.804,.769,.739,.734,.693,.699,.735,.894, *.796,.864,.688,.665,.563,.623,.582,.560,.559,.588,.596,.646, *.616,.637,.715,.685,.650,.684,.489,.484,.461,.463,.480,.512/ data (vco2(zi),zi = 338,553)/ *.545,.583,.560,.559,.528,.505,.532,.522,.415,.414,.387,.394, *.420,.463,.512,.562,.532,.526,.465,.435,.484,.461,.417,.457, *.583,.519,.568,.602,.614,.588,.596,.542,.574,.540,.403,.435, *.505,.448,.917,.864,.750,.796,.753,.653,.713,.688,.625,.699, *.670,.708,.685,.725,.903,.880,.782,.855,.438,.406,.313,.359, *.323,.296,.285,.303,.295,.310,.336,.421,.366,.410,.188,.184, *.063,.150,.100,.088,.105,.157,.198,.271,.227,.239,.257,.213, *.220,.230,.077,.116,.133,.157,.179,.185,.216,.208,.213,.188, *.069,.074,.144,.096,.225,.255,.281,.296,.299,.282,.281,.250, *.269,.253,.194,.213,.241,.216,.354,.350,.396,.359,.378,.380, *.369,.345,.323,.292,.310,.309,.319,.338,.324,.327,.479,.452, *.437,.428,.424,.382,.407,.387,.312,.373,.341,.333,.338,.355, *.444,.449,.394,.429,.0 ,.012,.0 ,.019,.012,.019,.031,.046, *.068,.074,.068,.019,.046,.031,.0 ,.012,.0 ,.019,.012,.019, *.031,.046,.068,.083,.074,.068,.028,.019,.046,.031,.059,.088, *.114,.130,.133,.116,.114,.083,.102,.086,.028,.046,.074,.049, *.207,.227,.262,.269,.253,.213,.179,.125,.157,.151,.153,.185, *.171,.170,.292,.262,.208,.227,.207,.185,.170,.171,.151,.167/ data (vco2(zi),zi = 554,673)/ *.157,.179,.278,.269,.213,.253,.042,.059,.125,.088,.114,.153, *.130,.133,.125,.116,.114,.083,.102,.086,.028,.046,.074,.049, *.0 ,.012,.0 ,.019,.012,.019,.031,.046,.068,.074,.068,.019, *.046,.031,.0 ,.012,.0 ,.019,.012,.019,.031,.046,.068,.083, *.074,.068,.028,.019,.046,.031,.042,.059,.125,.088,.114,.130, *.133,.116,.114,.083,.102,.086,.028,.046,.074,.049,.208,.207, *.292,.227,.262,.269,.253,.213,.179,.125,.157,.151,.153,.185, *.171,.170,.292,.262,.208,.227,.207,.185,.170,.171,.151,.167, *.157,.179,.278,.269,.213,.253,.042,.059,.125,.088,.114,.153, *.130,.133,.125,.116,.114,.083,.102,.086,.028,.046,.074,.049/ c 670 after subdivision: 4 c--------------------------------------------------------------------- npairs = 0 if (iscp.eq.3) goto 30 c one-dimensional simplex: if (imod.eq.0) then c cartesian imod = 0: call xchk (xmin1, xmax1, sname) loop1 = idint ((xmax1-xmin1) / xinc1 + 1.0000001) if (xmax1-xmin1.le.0.d0) loop1 = 1 x = xmin1 c do 40 j = 1, loop1 npairs = npairs+1 c generate compositions: xy(1,npairs) = x 40 x = x + xinc1 c else if (imod.eq.1) then c symmetric stretching imod = 1: b = xmin1 if (b.lt.1.d0) call error (207,b,j,sname) bm1 = b - 1. bp1 = b + 1. dx = 1. / (xinc1 - 1.) x = dx loop = idint (xinc1) - 1 c first loop generates phases c with xx(1) = 0 through 0.5 . do 50 j = 1, loop npairs = npairs+1 c generate compositions: t1 = (bp1/bm1)**(1.-x) xy(1,npairs) = 0.5 * (bp1-bm1*t1)/(1.+t1) 50 x = x + dx c second loop generates phases c with y = 0.5 through 1 - xmin. loop = loop - 1 jpairs = loop do 60 j = 1,loop npairs = npairs+1 c generate compositions: xy(1,npairs) = 1. - xy(1,jpairs) 60 jpairs = jpairs - 1 c add endmember comps just to c be safe, this is only for two c site models. npairs = npairs + 2 xy(1,npairs ) = 0.0 xy(1,npairs-1) = 1.0 c else if (imod.eq.2) then c assymmetric stretching imod = 2: b = xmin1 if (b.lt.1.d0) call error (207,b,j,sname) bm1 = b - 1. bp1 = b + 1. dx = 1. / (xinc1 - 1.) x = dx loop = idint (xinc1) - 1 c with x = xmin through xmax do 70 j = 1, loop npairs = npairs+1 c generate compositions: t1 = (bp1/bm1)**(1.-x) xy(1,npairs) = xmax1 * (bp1-bm1*t1)/(1.+t1) 70 x = x + dx c value of imod is illegal, c write error msg: else if (imod.gt.3) then call error (169,b,imod,'CART') end if goto 99 c two dimensional simplex: c check valid limits: 30 call xchk (xmin2, xmax2, sname) call xchk (xmin1, xmax1, sname) if (xmax1.gt.(1. - xmin2)) then call warn (170, x, imod, 'CART') xmax1 = 1. - xmin2 end if if (xmax2.gt.(1.-xmin1)) then call warn (171, x, imod, 'CART') xmax2 = 1. - xmin1 end if if (imod.eq.0) then c cartesian subdivision: loop1 = idint ((xmax1-xmin1)/xinc1 + 1.000001) loop2 = idint ((xmax2-xmin2)/xinc2 + 1.000001) x = xmin1 do 10 i = 1, loop1 y = xmin2 do 20 j = 1, loop2 c generate compositions: if (x+y.gt.1.0) goto 10 npairs = npairs + 1 if (npairs.gt.k1) call error (180,b,k1,'CART') xy(1,npairs) = x xy(2,npairs) = y 20 y = y + xinc2 10 x = x + xinc1 else if (imod.gt.0.and.imod.lt.5) then c barycentric subdivision of ternary c solutions, imod = # of iterations. npairs = ibary(imod) if (npairs.gt.k1) call error (180,b,k1,'CART') dx1 = xmax1 - xmin1 dx2 = xmax2 - xmin2 do 100 j = 1, npairs xy(1,j) = xmin1 + ( vco1(j) - xmin1 ) * dx1 100 xy(2,j) = xmin2 + ( vco2(j) - xmin2 ) * dx2 else if (imod.gt.4) then c illegal value of imod, write error msg: call error (169,b,imod,'CART') end if 99 end subroutine xchk (xmin, xmax, tname) implicit double precision (a-g,o-y),integer (h-n,z) character tname*10 if (xmax.gt.1.0) then call warn (21,xmax,1,tname) xmax = 1. end if if (xmin.lt.0.0) then call warn (22,xmin,1,tname) xmin = 0. end if if (xmax.lt.xmin) then call warn (23,xmax,1,tname) xmax = 1. xmin = 0. end if end subroutine zchk (z, tname) implicit double precision (a-g,o-z),integer (h-n) character tname*10 if (z.gt.0.0.and.z.le.1.d-3) then z = 0.0 else if (z.ge.0.999.and.z.lt.1.0) then z = 1.0 else if (z.ge.0.0.and.z.le.1.0) then return else if (z.gt.1.0.and.z.lt.1.01) then z = 1.0 else if (z.gt.-.01) then z = 0.0 else call error (125,z,1,tname) end if end subroutine soload (z,isolct,tname) implicit double precision (a-g,o-z),integer (h-n) parameter (k1=1000,k2=800,k3=800,k4=16,k5=12,k6=7,m0=5) parameter (m9=50,k9=50,m6=3,m7=9,m8=9,h5=7,h6=500,nb=20) character*10 tname, names*8 double precision z(3,3) common / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 common/ cst108 /wg(27,11),xmn(3,3),xmx(3,3),xnc(3,3),xfx(3), * iend(4,3,3),isub(27,9,2),imd(3),imsol(4,3,3), * ijk(3),insp(3,4),jsp(3),jsitp(3),iosp(3), * jst(3),jfix,idfc(3),ist(3),isp(3),isite, * ispct(3,4),jnsp(3,4),jsite,iterm,iord common/ cst107 /a0(5,5),acoef(5,5,m0),smult(5),ecoef(3,3,3), * nsite,nspm1(5),nterm(5,5), * nsub1(5,5,m0,4),nsub2(5,5,m0,4),nttyp(5,5,m0) common/ cst1 /thermo(k4,k1),vf(2,k1),vs(h5,k1),uf(2),us(h5) * / cst202 /tm(m7,m6),td(m8),ilam,idiso,lamin,idsin * / cst203 /therdi(m8,m9),therlm(m7,m6,k9) * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn * / cst20 /comps(k1,k5)/ cst8 /names(k1) * / cst204 /ltyp(k1),lmda(k1),idis(k1) * / cst6 /icomp,istct,iphct,icp * / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar * / cst5 /p,t,xco2,u1,u2,tr,pr,r,ps/ cst12 /cp(k6,k1) * / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2 * / cst88 /vu(2,k1),jprct,jmct,jmyn * / cst300 /cblk(nb,k5),atw(k1),vol(k1),iblk,jcby,jcth c eliminate end-member compositions do 10 i = 1, isp(1) call zchk (z(1,i),tname) do 20 j = 1, isp(2) call zchk (z(2,j),tname) do 30 k = 1, isp(3) call zchk (z(3,k),tname) 30 if (z(1,i)*z(2,j)*z(3,k).gt.0.999) goto 999 20 continue 10 continue c eliminate out of range compositions: c is this necessary? do 50 i = 1, isite if (imd(i).gt.0) goto 50 do 60 j = 1, isp(i) - 1 60 if (z(i,j).gt.xmx(i,j).or.z(i,j).lt.xmn(i,j)) goto 999 50 continue c the composition is acceptable. iphct = iphct + 1 jltyp = -1 if (iphct.gt.k1) call error (180,z(1,1),k1,'PRIME9') ikp(iphct) = isolct idis(iphct) = 0 lmda(iphct) = 0 c encode a name if (isite.eq.1) then if (isp(1).eq.2) then if (z(1,1).eq.1.0) then write (names(iphct),1050) names(imsol(1,1,1)),100 else if (z(1,1).ge.0.04.and.z(1,1).le.0.96) then write (names(iphct),1030) names(imsol(1,1,1)), * idint(1.d2*z(1,1)) else write (names(iphct),1070) names(imsol(1,1,1)), * 100.*z(1,1) end if else c ternary solutions write (names(iphct),1040) names(imsol(1,1,1)), * idint(z(1,1)*1.d2),names(imsol(2,1,1)),idint(z(1,2)*1.d2) end if else if (isite.eq.2.and.isp(1).eq.2.and.isp(2).eq.2) then int1 = idint (1.d2*z(1,1)) int2 = idint (1.d2*z(2,1)) if (int1.gt.99) int1 = 99 if (int2.gt.99) int2 = 99 write (names(iphct),1020) names(imsol(1,1,1)),int1,int2 else c this is a pretty shitty name c but for the time being your going to c have live with it. write (names(iphct),1010) ((idint(1.d1*z(i,j)), * j = 1, isp(i) - 1), i = 1, isite) write (n3,1000) names(iphct),(((z(1,i)*z(2,j)*z(3,k), * names(imsol(i,j,k)),i=1,isp(1)), * j=1,isp(2)), * k=1,isp(3)) end if c get blanks out of name: call unblnk (names(iphct),8) c initialize constants: smix = 0.0 esum = 0. atw(iphct) = 0. vol(iphct) = 0. do 130 i = 1, icp 130 cp(i,iphct) = 0. do 140 i = 1, k4 140 thermo(i,iphct) = 0. do 150 i = 1, icomp 150 comps(iphct,i) = 0. c load constants: do 160 i = 1, isp(1) if (z(1,i).eq.0.) goto 160 do 170 j = 1, isp(2) if (z(2,j).eq.0.) goto 170 do 180 k = 1, isp(3) if (z(3,k).eq.0.) goto 180 zp = z(1,i)*z(2,j)*z(3,k) id = imsol(i,j,k) esum = esum + zp * ecoef(i,j,k) atw(iphct) = atw(iphct) + zp * atw(id) vol(iphct) = vol(iphct) + zp * vol(id) do 190 l = 1, icp 190 cp(l,iphct) = cp(l,iphct) + zp * cp(l,id) do 200 l = 1, icomp 200 comps(iphct,l) = comps(iphct,l) + zp * comps(id,l) do 210 l = 1, k4 210 thermo(l,iphct) = thermo(l,iphct) * + zp * thermo(l,id) c t dependent disordering: if (idis(id).ne.0) then idsin = idsin + 1 idis(iphct) = idsin do 330 l = 1, 6 330 therdi(l,idsin) = zp * therdi(l,id) do 340 l = 7, 9 340 therdi(l,idsin) = therdi(l,id) end if c lamda-like transitions: ld = lmda(id) if (ld.ne.0) then if (jltyp.eq.-1) lamin = lamin + 1 jltyp = jltyp + 1 ltyp(iphct) = ltyp(id) lmda(iphct) = lamin ld = lmda(id) if (ltyp(id).eq.1) then c ubc-approach: do 350 l = 1, 3 therlm(1,l,lamin) = therlm(1,l,ld)*zp therlm(2,l,lamin) = therlm(2,l,ld)*zp therlm(3,l,lamin) = therlm(3,l,ld) therlm(4,l,lamin) = therlm(4,l,ld) therlm(5,l,lamin) = therlm(5,l,ld)*zp therlm(6,l,lamin) = therlm(6,l,ld)*zp therlm(7,l,lamin) = therlm(7,l,ld) therlm(8,l,lamin) = therlm(8,l,ld)*zp 350 therlm(9,l,lamin) = therlm(9,l,ld)*zp c else if (ltyp(id).lt.8) then c helgeson approach: if (ltyp(id).gt.4.and.ltyp(id).lt.7) then jlam = ltyp(id) - 3 else jlam = 1 end if do 365 m = 1, jlam therlm(1,m,lamin) = therlm(1,m,ld) therlm(2,m,lamin) = therlm(2,m,ld) do 360 l = 3, 8 360 therlm(l,m,lamin) = therlm(l,m,ld)*zp 365 continue else if (ltyp(id).eq.10) then c holland and powell landau model: h = jltyp + 1 therlm(1,h,lamin) = therlm(1,1,ld) do 222 m = 2, 5 222 therlm(m,h,lamin) = therlm(m,1,ld)*zp end if end if 180 continue 170 continue 160 continue c sort phase into subcomposition call sort (ier) c if ier = 1, the phase consists c entirely of saturated components: if (ier.ne.0) call satsrt c compute ideal configurational negentropy: if (nsite.eq.0.and.ist(1).ne.0) then do 220 i = 1, isite do 230 j = 1, isp(i) 230 if (z(i,j).ne.0.0) smix = * smix + float(ist(i))*z(i,j)*dlog(z(i,j)) 220 continue else if (nsite.ne.0) then call omega (z,dlnw,tname) smix = -dlnw + esum end if c save it: thermo(2,iphct) = thermo(2,iphct) + r * smix goto (290),isyn c load stoichiometric coefficients of c saturated components: do 300 k = 1, isat 300 vs(k,iphct) = comps(iphct,icp+k) c load stoichiometric coefficients of c mobile components: 290 goto (310),jmyn do 320 i = 1, jmct 320 vu(i,iphct) = comps(iphct,jprct+i) 310 if (iff(1).ne.0) vf(1,iphct) = comps(iphct,iff(1)) if (iff(2).ne.0) vf(2,iphct) = comps(iphct,iff(2)) c load excess terms: if (iterm.ne.0) then do 240 i = 1, iterm zp = 1.0 do 250 j = 1, iord if (isub(i,j,1).eq.0) goto 250 zp = zp * z(isub(i,j,1),isub(i,j,2)) 250 continue do 260 j = 1, 11 260 thermo(j,iphct) = thermo(j,iphct) + zp * wg(i,j) 240 continue end if 1000 format (' icky name ',a8,' = ',4(4(f5.3,1x,a4,1x),/)) 1010 format (8i1) 1020 format (a2,i2,'-',i2) 1030 format (a2,i2) 1040 format (a1,i2,a1,i2) 1050 format (a2,i3) 1070 format (a2,f4.1) 999 end subroutine assptx c--------------------------------------------------------------------- c assptx assigns the current value of v1 and v2 to the array ptx and c increments the ordinate counter ipt2. c--------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (l2=5) common/ cst32 /ptx(450),logx,ipt2 * / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5 * / cst5 /v(l2),tr,pr,r,ps/ cst63 /delv(l2) * / cst109 /xxx(450),x1,x2,jpt2 c delete near replicates 10 ipt2=ipt2+2 ptx(ipt2-1) = v(iv1) ptx(ipt2) = v(iv2) xxx(ipt2 - 1) = x1 xxx(ipt2) = x2 99 end subroutine concrt c--------------------------------------------------------------------- c concrt determines convergence criteria and limits for reasonable c solutions for the routine univeq. the array delt is also used in c the routine slope. c input: dv - array of default search increments. c vmax,vmin - arrays containing the maximum and minimum values c of the independent and dependent intensive variables. c iv - array indexes variables in vmax, vmin, and output arrays c output: ulim, blim - arrays with upper and lower limits for reasonable c solutions (vmax,vmin+/-4*dv) c delt - array containing the finite difference increments c (vmax-vmin/10**5), delt also serves as the test value c for convergence. c--------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (l2 = 5) common/ cst62 /blim(l2),ulim(l2),dgr * / cst9 /vmax(l2),vmin(l2),dv(l2) c--------------------------------------------------------------------- do 10 i = 1,5 if (dv(i).lt.0.d0) call error (34,dv(i),i,'CONCRT') ulim(i) = vmax(i)+dv(i) blim(i) = vmin(i)-dv(i) ddv = vmax(i)-vmin(i) 10 if (ddv.lt.0.d0) call error (35,ddv,i,'CONCRT') end subroutine conver (g,s,v,a,b,c,d,e,f,gg,b2,b3,b4,b5,b6,b7) c--------------------------------------------------------------------- c convert thermodynamic equation of state from a pr tr reference state c to get polynomial coefficients c--------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) common / cst5 /p,t,xco2,u1,u2,tr,pr,r,ps save c1, c2 data c1, c2/ -300.d0, -35000.d0/ c first compute G* g = g c add the integral of sdt at tr,pr: * + s * tr - a * tr - b * tr * tr / 2.d0 + c / tr * - d * tr**3 / 3.d0 - 2.d0 * e * sqrt(tr) * - f * dlog(tr) + gg / tr / tr / 2.d0 + f c add the constant components of the c integral -vdp at (t,pr): * - v * pr + b2 * tr * pr + b3 * exp( tr / c1) * pr * + b4 * pr * pr / 2.d0 + b5 * pr * dexp( pr / c2) * - b5 * c2 * dexp ( pr / c2 ) - b6 * pr**3/3.d0 * + b7 * tr * tr * pr c S* (T) s = a - b2 * pr - s + a * dlog(tr) + b * tr - c / tr / tr /2.d0 * + d * tr * tr / 2.d0 - 2.d0 * e / sqrt(tr) - f / tr * - gg / tr**3 / 3.d0 c V* (P) v = v - b2 * tr - b3 * exp(tr/c1) * - b4 * pr - b5 * dexp(pr/c2) - b7 * tr * tr c a* (-T log T) c a = a c b* (-T*T) b = b7 * pr + b / 2.d0 c c* (-1/T) c = c / 2.d0 c d* (-T**3) d = d / 6.d0 c e* (sqrt T) e = 4.d0 * e c f* (log T) c f = f c gg* (-1/T**2) gg = gg/6.d0 c b2* (PT) c b2 = b2 c b3* ((P-Pr) exp (T / c1)) c b3 = b3 c b4* (P**2) b4 = b4/2.d0 c b5* (exp (P/ c2)) b5 = c2 * b5 c b6* (P**3) b6 = b6/3.d0 c b7* (P*T**2) c b7 = b7 end subroutine disord (tk,d0,d1,d2,d3,d4,d5,d6,tr,tmax,g) c---------------------------------------------------------------------- c compute t-dependent disorder contribution, g is the c gibbs energy of disordering, tr is the t of onset of c disordering, tmax is the t of complete disorder. c--------------------------------------------------------------------- implicit double precision (a-g,o-z),integer (h-m) common/ cst5 /p,v(4),pr,z(3) g = 0.d0 if (tk .lt. tr) return tt = tk if (tk.gt.tmax) tt = tmax dh = d0 * (tt - tr) * + 2.d0 * d1 * (tt**0.5 - tr**0.5) * - d2 * (1.d0 / tt - 1.d0 / tr) * + d4 * dlog(tt/tr) * + d5 * (tt*tt - tr*tr) / 2.d0 * + d6 * (tt**3 - tr**3) / 3.d0 ds = d0 * (dlog(tt) - dlog(tr)) * - 2.d0 * d1 * (tt**(-0.5) - tr**(-0.5)) * - d2 * (1.d0/tt/tt - 1.d0/tr/tr) / 2.d0 * - d4 * (1.d0/tt - 1.d0 / tr) * + d5 * (tt - tr) * + d6 * (tt*tt - tr*tr) / 2.d0 gss = dh - (tk * ds) if (d3.ne.0.0) gss = gss + dh/d3 * (p - pr) end function gphase(id) c----------------------------------------------------------------------- c gphase computes the gibbs free energy of a phase identified by c the arguement 'id' from the thermochemical parameters stored c in the array 'thermo' which is located in common block 'cst1'. c the parameters are: g(pr,tr) [thermo(1,id)], s(pr,tr) [thermo c (2,id)], v(pr,tr) [thermo(3,id)], the coefficients of the extended c maier-kelley heat capacity equation: c cp(pr,t)=a+b*t+c/(t*t)+d*t*t+e/t**(1/2)+f/t+g/t**3 c [thermo(id,4-10)], and the coefficients of the haas molar volume c equation (modified for a reference state v(pr,tr)): c v(p,t)=v(pr,tr)+b2*(t-tr)+b3*(exp(-t/300)-exp(tr/300)) c +b4*(p-pr)+b5*(exp(-p/35000)-exp(-pr/35000)) c +b6*(p-pr)**2+b7*(t-tr)**2 c [thermo(11-16,id)]. c--------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1000,k4=16,k6=7,m9=50,m6=3,m7=9,k9=50,m8=9, * h5=7,k5=12) common/ cst1 /thermo(k4,k1),vf(2,k1),vs(h5,k1),uf(2),us(h5) * / cst88 /vu(2,k1),jprct,jmct,jmyn/ cst11 /f(2) * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn * / cst5 /p,t,xco2,u1,u2,tr,pr,r,ps * / cst44 /dlogt,dsqrtt,dexpt,dexpp * / cst203 /therdi(m8,m9),therlm(m7,m6,k9) * / cst42 /ic(k5),idbase/ cst12 /cp(k6,k1) * / cst204 /ltyp(k1),lmda(k1),idis(k1) save c1, c2 data c1, c2/ -300.d0, -35000.d0/ gphase = thermo(1,id) c -sdt * + t * (thermo(2,id) - thermo(4,id) * dlogt * - t * (thermo(5,id) + thermo(7,id) * t)) * - (thermo(6,id) + thermo(10,id) / t) / t * + thermo(8,id) * dsqrtt + thermo(9,id)*dlogt c +vdp (vs(p,t)) c the term thermo(12,id) should c actually be multiplied by p-pr, c this won't make much difference c provided p is far from the pr and c pr is large. * + p * (thermo(3,id) + thermo(12,id) * dexp (t/c1) * + (thermo(16,id) * t + thermo(11,id)) * t * + (thermo(15,id) * p + thermo(13,id)) * p) * + thermo(14,id) * dexp (p/c2) c -ndu * - vu(1,id) * u1 - vu(2,id) * u2 c check for lambda transitions: ltype = ltyp(id) if (ltype.ne.0) then if (ltype.lt.4) then c ubc-type transitions call lamubc (p,t,dg,lmda(id),ltype) gphase = gphase + dg else if (ltype.lt.7) then c standard transitions du = - vu(1,id) * u1 - vu(2,id) * u2 call lamhel (p,t,gphase,du,lmda(id),ltype,id) else if (ltype.lt.10) then c supcrt q/coe lambda transition du = - vu(1,id) * u1 - vu(2,id) * u2 call lamqtz (p,t,gphase,du,lmda(id),id) else if (ltype.gt.9 .and. ltype.le.12 ) then c holland and powell (landau model) call lamlan (t,dg,lmda(id),ltype) gphase = gphase + dg end if end if c check for temperature dependent c order/disorder: if (idis(id).ne.0) then jd = idis(id) call disord (t,therdi(1,jd),therdi(2,jd),therdi(3,jd), * therdi(4,jd),therdi(5,jd),therdi(6,jd), * therdi(7,jd),therdi(8,jd),therdi(9,jd),gss) gphase = gphase + gss end if c modification for fluid c projections: if (ifyn.eq.1.and.thermo(3,id).eq.0.0) then c if pr or tr = 0 then c this is a special calculation c and we can skip this part. if (pr.eq.0.0.or.tr.eq.0.0) goto 999 if (cp(1,id).eq.0.0.and.cp(2,id).eq.0.0) goto 999 xco2 = cp(2,id) call cfluid (fo2,fs2) gphase = gphase + r*t*(cp(1,id)*f(1) + cp(2,id)*f(2)) end if 999 end subroutine getphi (name,ibase,icmpn,*) c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k0=25,k4=16,k5=12,m6=3,m7=9,m8=9) character*8 name,oname common/ cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst43 /therm(k4),comp(k0),atwt(k0),idh2o,idco2,ikind * / cst202 /tm(m7,m6),td(m8),ilam,idiso,lamin,idsin * / cst207 /ctrans(k0,k5),ictr(k5),itrans save oname data oname/' '/ 30 read (n2,1020,end=90,err=98) * name, ibase, ikind, ilam, idiso if (name.eq.' ') goto 30 read (n2,*,err=98) (comp(i), i=1, icmpn), therm c do component transformation if c itrans is not zero if (itrans.gt.0) then do 10 i = 1, itrans it = ictr(i) if (comp(it).eq.0.d0) goto 10 c ct is how much of the new c component is in the phase. ct = comp(it) / ctrans(it,i) do 20 j = 1, icmpn 20 comp(j) = comp(j) - ct * ctrans(j,i) comp(it) = ct 10 continue end if if (ilam.ne.0) then c determine number of transitions from c flag ilam: jlam=ilam if (ilam.gt.3) jlam = ilam-3 if (ilam.gt.6) jlam = ilam-6 if (ilam.gt.9) jlam = ilam-9 do 1 i= 1, jlam 1 read (n2,*,err=98) (tm(j,i), j = 1, m7) end if if (idiso.ne.0) read (n2,*,err=98) td oname = name return 90 return 1 98 if (oname.ne.' ') then call error (23,ct,i,oname) else call error (23,ct,i,' none ') end if 1020 format (a8,i2,i2,i2,i2) end subroutine chkphi (ichk,name,icout,ibase,icmpn,*) c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k0=25,k4=16,k5=12,h8=100) character*8 name, exname integer icout(k0) common/ cst43 /therm(k4),comp(k0),atwt(k0),idh2o,idco2,ikind * / cst6 /icomp,istct,iphct,icp/ cst36 /exname(h8) * / cst37 /iprct,ixct,iexyn,istbyn/ cst42 /ic(k5),idbase c reject the data if wrong data base: if (ibase.ne.idbase) goto 90 if (iexyn.eq.0) then c reject the data if excluded in input: do 1 i = 1, ixct 1 if (name.eq.exname(i)) goto 90 end if if (ichk.eq.0) return c if ichk ne 0 check for components: c reject phases which have a component c not included in the input: do 2 i= 1, icmpn 2 if ((icout(i).eq.0).and.(comp(i).ne.0.d0)) goto 90 c reject phases which consist entirely c of constrained components: if (iprct.eq.0) goto 20 do 3 i= 1, icp 3 if (comp(ic(i)).ne.0.d0) goto 20 goto 90 20 return 90 return 1 end subroutine lamlan (t,dg,ld,ltype) c--------------------------------------------------------------------- c calculate the extra energy of a lamdba transition using the c Landau model (Holland and Powell). c input variables c t = temperature in k c ld = pointer to phase in therlm c id = pointer to phase in thermo c returned - dg - difference in free energy due to the transition c--------------------------------------------------------------------- implicit double precision (a-g,o-z),integer (h-m) parameter (m9=50,k9=50,m6=3,m7=9,m8=9) common/ cst203 /therdi(m8,m9),therlm(m7,m6,k9) lct = ltype - 9 dg = 0.0 do 10 i = 1, lct tc = therlm(1,i,ld) smax = therlm (2,i,ld) c subtract the reference state c contribution first: dg = dg + t*therlm(5,i,ld) - therlm(4,i,ld) if (t.gt.tc) then c fully disordered: dg = dg + therlm(3,i,ld) - t * smax else c partially disordered: q = (1.-(t/tc))**0.25 dg = dg + smax * (q*q*(t-tc) + tc*(2.+q**6)/3. - t) end if 10 continue end subroutine lamhel (p,t,g,du,ld,ltype,id) c--------------------------------------------------------------------- c calculate the extra energy of a lamdba transition using model c of helgeson et al 1978 (AJS). c there is something seriously wrong in this routine!!!! c input variables c p = pressure in bars c t = temperature in k c ld = pointer to phase in therlm c g = initial free energy c returned - g - modified phase free energy c--------------------------------------------------------------------- implicit double precision (a-g,o-z),integer (h-m) parameter (m9=50,k9=50,m6=3,m7=9,m8=9,k1=1000,k4=16,h5=7) common/ cst1 /thermo(k4,k1),vf(2,k1),vs(h5,k1),uf(2),us(h5) * / cst203 /therdi(m8,m9),therlm(m7,m6,k9) * / cst5 /v(6),pr,u(2) trt = therlm(1,1,ld) c below 1 bar transition T, don't do c anything except for q/coe case: if (t.lt.trt) goto 99 c find out which transition range j = 1 if (ltype.gt.4) then ilam = ltype - 3 do 10 i = 1, ilam if (t.gt.therlm(1,i,ld)) then j = i trt = therlm(1,i,ld) end if 10 continue end if c T > 1 bar transition T: c integrate cp function to t from tr g = therlm(8,j,ld) - therlm(3,j,ld) * (t-trt) + du * + therlm(5,j,ld) * (t - trt - t * dlog(t/trt)) * - (therlm(7,j,ld) + therlm(6,j,ld) * t * trt * trt) * * (t - trt)**2 / 2.d0 / t / trt / trt * + therlm(9,j,ld) * * ( t**3/3. - trt**3/3. + t*t/2. - trt**2/2.) c add pdv terms if (therlm(2,1,ld).ne.0.0) then trtp = trt + (p-pr)/therlm(2,1,ld) else trtp = trt end if if (t .gt. trtp) then c 1 bar polymorph is stable at p, c p*thermo(3,id) is 0 bar standard c state pdv polymorph contribution & c therlm(4,j,ld) is Delta V of trans c contribution. pstar = (t-trt)*therlm(2,1,ld) + pr g = g + p * thermo(3,id) + (p-pstar)*therlm(4,j,ld) else c the 1 bar polymorph isn't stable. g = g + p * thermo(3,id) * + (t-trt) * therlm(2,1,ld) * therlm(4,j,ld) end if 99 end subroutine lamqtz (p,t,g,du,ld,id) c--------------------------------------------------------------------- c calculate the extra energy of a lamdba transition using model c of helgeson et al 1978 (AJS). eq. (114) (corrected for the c misplaced parenthesis and the sign on c sub alpha), this is C probably bogus for coesite note that only one transition c is allowed. c input variables c p = pressure in bars c t = temperature in k c ld = pointer to phase in therlm c g = initial free energy c constants c Tr = 298.15 K c Pr = 1. bar c S(lope) = 38.5 bar/K c ba = 0.066 j/bar c aa = 549.82 K c ca = -4.973d-6 j/bar/bar c trt = 848. K c returned - g - modified phase free energy c--------------------------------------------------------------------- implicit double precision (a-g,o-z),integer (h-m) parameter (m9=50,k9=50,m6=3,m7=9,m8=9,k1=1000,k4=16,h5=7) common/ cst1 /thermo(k4,k1),vf(2,k1),vs(h5,k1),uf(2),us(h5) * / cst203 /therdi(m8,m9),therlm(m7,m6,k9) save trt, tr, pr, s, aa, ba, ca , vb data trt, tr, pr, s, aa, ba, ca, vb / 848., 298.15, 1.0, * 38.5, 549.82, 0.066, -4.973d-6, 2.372 / c trtp = trt + (p-pr) / s if (t.gt.trt) then ps = (t-trt) * therlm(2,1,ld) + pr else ps = pr end if c if above the ref P transition T c do the cp integration: if (t.gt.trt) g = therlm(8,1,ld) + (p-ps) * thermo(3,id) * - therlm(3,1,ld) * (t-trt) + du * + therlm(5,1,ld) * (t - trt - t * dlog(t/trt)) * - (therlm(7,1,ld) + therlm(6,1,ld) * t * trt * trt) * * (t - trt)**2 / 2.d0 / t / trt / trt c now add in pdv terms to c the free energy, note that c pdv term for the ref polymorph c is already in: pdv = vb * (ps-pr) * - ca * ( (2. * pr * (p-ps) - (p * p - ps * ps) ) / 2. * + s * (t-tr) * (p-ps) ) * + s * (ba + aa*ca*s) * (t-tr) * * dlog ((aa+p/s) / (aa + ps/s)) g = g + pdv end subroutine lamubc (p,t,gspk,k,ltype) c--------------------------------------------------------------------- c calculate the extra energy of a lamdba transition using model c of berman and brown (1985, contribs. miner. petro.) c input variables c p = pressure in bars c t = temperature in k c k = pointer to phase in therlm c returned - gspk - energy because of lambda spike c modified from perkins et al 1986 (computers and geosciences). c--------------------------------------------------------------------- implicit double precision (a-g,o-z),integer (h-m) parameter (m9=50,k9=50,m6=3,m7=9,m8=9) common/ cst203 /therdi(m8,m9),therlm(m7,m6,k9) c--------------------------------------------------------------------- gspk=0.d0 do 10 j = 1, ltype tq1bar = therlm(3,j,k) if (tq1bar.eq.0.d0) goto 10 tr = therlm(7,j,k) teq = therlm(4,j,k) * (p-1.d0) + tq1bar ctrans = tq1bar - teq tr9 = tr - ctrans if (t.lt.tr9) goto 10 aspk2 = therlm(1,j,k) bspk2 = therlm(2,j,k) dvdtr = therlm(5,j,k) dvdp = therlm(6,j,k) dstr = therlm(8,j,k) / tq1bar abspk = therlm(9,j,k) if (t .gt. teq) then t9 = teq else t9 = t end if ct2 = ctrans * ctrans ct3 = ct2 * ctrans a1 = aspk2 * ctrans + 2. * abspk * ct2 + bspk2 * ct3 b1 = aspk2 + 4. * abspk * ctrans + 3. * bspk2 * ct2 c1 = 2. * abspk + 3. * ctrans * bspk2 t92 = t9 * t9 t93 = t9 *t92 tr92 = tr9 * tr9 tr93 = tr92 * tr9 dhspk = a1 * (t9 - tr9) * + b1 * (t92 - tr92) / 2. * + c1 * (t93 - tr93) / 3. * + bspk2 * (t9*t93 - tr93*tr9) / 4. dsspk = a1 * (dlog(t9) - dlog(tr9)) * + b1 * (t9 - tr9) * + c1 * (t92 - tr92) / 2. * + bspk2 * (t93 - tr93) / 3. gspk = gspk - (t9 * dsspk) + dhspk if (t.gt.teq) gspk = gspk - (dstr + dsspk ) * (t-teq) gspk = gspk + dvdtr * (p-1.) * (t9-tr) * + dvdp/2. * (p*p-1.) - dvdp * (p-1.) 10 continue end subroutine loadit (id) c--------------------------------------------------------------------- c loadit loads descriptive data for phases and species (name,comp, c and therm) into the appropriate arrays (names,comps,thermo,vf, c and vs). the arguement 'id' indexes the phase in the arrays. c note that loadit also computes reference state constants which c are dependent on the state function being used and on its c analytical expression. c--------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k0=25,k1=1000,k2=800,k3=800,k4=16,k5=12,k6=7) parameter (m9=50,k9=50,m6=3,m7=9,m8=9,h5=7,h6=500,j9=999) parameter (nb=20) character*8 names, name, cmpnt*5, dname*40 common/ cst60 /iasmbl(j9),ipoint,imyn/ csta6 /name * / cst88 /vu(2,k1),jprct,jmct,jmyn * / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar * / cst1 /thermo(k4,k1),vf(2,k1),vs(h5,k1),uf(2),us(h5) * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn * / cst20 /comps(k1,k5)/ cst8 /names(k1) * / cst5 /p,t,xco2,u1,u2,tr,pr,r,ps/ cst12 /cp(k6,k1) * / cst6 /icomp,istct,iphct,icp * / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2 * / csta5 /dname(8),cmpnt(k0)/ cst42 /ic(k5),idbase * / cst202 /tm(m7,m6),td(m8),ilam,idiso,lamin,idsin * / cst203 /therdi(m8,m9),therlm(m7,m6,k9) * / cst43 /therm(k4),comp(k0),atwt(k0),idh2o,idco2,ikind * / cst300 /cblk(nb,k5),atw(k1),vol(k1),iblk,jcby,jcth * / cst204 /ltyp(k1),lmda(k1),idis(k1) c--------------------------------------------------------------------- c load stoichiometry of saturated c components in the phase. goto (20),isyn do 10 i = 1, isat 10 vs(i,id)=comp(ic(icp+i)) c load stoichiometric coefficients of c the variable components in the saturated c phase: 20 vf(1,id) = comp(idh2o) vf(2,id) = comp(idco2) c load stoichiometric coefficients of c mobile components: goto (50),jmyn do 60 i = 1, jmct 60 vu(i,id) = comp(ic(jprct+i)) c load name and phase flag 50 names(id)=name ikp(id)=ikind atw(id) = 0. c load stoichiometry of components. do 30 i=1,icomp atw(id) = comp(ic(i)) * atwt(ic(i)) + atw(id) 30 comps(id,i) = comp(ic(i)) c load stoichiometry of components c with unconstrained potentials. do 40 i = 1, icp 40 cp(i,id) = comp(ic(i)) c compute reference state constants, c etc.. do 1 i = 1, k4 1 thermo(i,id) = therm(i) c load ref volume for bulkck, assume j/bar vol(id) = therm(3)*10. c if 0 set to very big value to avoid c infinite mode. if (vol(id).eq.0.) vol(id) = 1.d30 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)) if (tr.eq.0) then thermo(1,id) = therm(1) thermo(2,id) = therm(2) end if c lmda transitions: idis(id) = 0 lmda(id) = 0 if (ilam.ne.0) then lamin=lamin+1 if (ilam.eq.10) then c holland and powell, landau model: tc = tm(1,1) smax = tm(2,1) therlm (1,1,lamin) = tc therlm (2,1,lamin) = smax c create a constant which contains c the 0 to Tr contribution to the c free energy: q = (1.-(tr/tc))**0.25 if (tr.gt.tc) q = 0. c maximum enthalpy of disorder: therlm (3,1,lamin) = 2./3.*smax*tc c enthalpy of disorder at tr: therlm (4,1,lamin) = 2.*smax*tc* (q**6/6. - q*q/2. + 1./3.) c entropy of disorder at tr: therlm (5,1,lamin) = smax * (1. - q*q) else if (ilam.le.3) then c ubc: do 90 j = 1, ilam therlm(1,j,lamin)=tm(1,j)*tm(1,j) therlm(2,j,lamin)=tm(2,j)*tm(2,j) therlm(9,j,lamin)=tm(1,j)*tm(2,j) do 80 k = 3, 8 80 therlm(k,j,lamin)=tm(k,j) 90 continue else if (ilam.gt.3.and.ilam.lt.8) then c helgeson: jlam = 1 c set transition type to null c for call to gphase lmda(id) = 0 c set P and T to Transition T and Pr p = pr t = tm(1,1) therlm(1,1,lamin) = tm(1,1) c get G(Transition T, Pr) call updumy therlm(8,1,lamin) = gphase(id) c get S(Transition T, Pr) t = t + 0.1d-2 call updumy c get delta V transition therlm(4,1,lamin) = 0. if (tm(2,1).ne.0.) therlm(4,1,lamin) = tm(2,1) / tm(3,1) therlm(3,1,lamin) = tm(3,1) - * (gphase(id) - therlm(8,1,lamin))/0.1d-2 c load dp/dt therlm(2,1,lamin) = tm(2,1) if (tm(2,1).lt.0.) call warn (12,r,i,name) c load remaining heat capacity terms: do 120 j = 4, 6 120 therlm(j+1,1,lamin) = tm(j,1) therlm(9,1,lamin) = tm(7,1) c if 3 < ilam < 7 load other transitions if (ilam.eq.4.or.ilam.eq.7) goto 130 jlam = ilam - 3 lmda(id)=lamin do 100 j = 2, jlam c increment transition type ltyp(id) = 2 + j c set transition T t = tm(1,j) therlm(1,j,lamin) = tm(1,j) c load dp/dt therlm(2,j,lamin) = tm(2,j) c get delta V transition, this will c lead to an erroneous pdv term. therlm(4,j,lamin) = 0. if (tm(2,j).ne.0.) therlm(4,1,lamin) = tm(2,j) / tm(3,j) if (tm(2,j).ne.0.) call warn (11,r,i,name) c get G at T trans call updumy therlm(8,j,lamin) = gphase(id) c get S at T trans t = t + .001d0 call updumy therlm(3,j,lamin) = tm(3,j) - * (gphase(id) - therlm(8,j,lamin))/0.001d0 c load cp terms, be careful of the c odd numbering. do 110 k = 4, 6 110 therlm(k+1,j,lamin) = tm(k,j) therlm(9,j,lamin) = tm(7,j) 100 continue end if 130 lmda(id)=lamin ltyp(id)=ilam end if c t dependent order: load berman and brown c parameters, (this should be cleaned up) if (idiso.ne.0) then idsin = idsin+1 idis(id) = idsin do 70 j = 1, m8 70 therdi(j,idsin) = td(j) end if end subroutine difgx (x,fact,index,iord,isub) c--------------------------------------------------------------------- c subroutine to evaluate the function P(X)+(1-xi)*diff(P(X),xi) c where P(X) is the product x(isub(1))*x(isub(2))*..., c isub contains the indexes of x in P(X), and index is the c value of i, the value is returned as difgx. iord is the maximum c order possible for P(X). c--------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) double precision x(4) integer isub(4), jnot(3), iord c if (x(index).eq.1.0.or.x(index).eq.0.0) * call error (27,x(1),iord,'IMBECILLE!') c find the order of P(X) = jord, c and count the number of occurences c of x(index) = ixi. ixi = 0 ini = 0 do 10 i = 1, iord if (isub(i).ne.0) then c assume 0 subscripts at end of isub jord = i if (isub(i).eq.index) then ixi = ixi + 1 else ini = ini + 1 jnot(ini) = isub(i) end if end if 10 continue c 2nd order parameter: if (jord.eq.2) then xm = 1. - x(index) if (ixi.eq.0) then fact = -x(isub(1))*x(isub(2)) else if (ixi.eq.1) then fact = x(jnot(1))*x(index) * + (xm*x(jnot(1))-2.*x(index)*x(jnot(1))) / xm else call error (999,fact,jord,'difgx') end if c 3rd order parameter: else if (jord.eq.3) then if (ixi.eq.0) then fact = -2.*x(isub(1))*x(isub(2))*x(isub(3)) else if (ixi.eq.1) then fact = -2.*x(isub(1))*x(isub(2))*x(isub(3)) * + x(jnot(1))*x(jnot(2)) else if (ixi.eq.2) then fact = 2.*(x(index)*x(jnot(1)) - x(index)**2*x(jnot(1))) else call error (999,fact,jord,'difgx') end if else if (jord.eq.4) then c 4th order parameter: if (ixi.eq.0) then call error (999,fact,jord,'difgx') else if (ixi.eq.1) then fact = -2.*x(index)**2*x(jnot(1))*x(jnot(2))*x(jnot(3)) * + 2. * x(isub(1))*x(isub(2))*x(isub(3))*x(isub(4)) else if (ixi.eq.2) then fact = 2.*x(index)*x(jnot(1))*x(jnot(2)) * - 3.*x(index)**2*x(jnot(1))*x(jnot(2)) else if (ixi.eq.3) then fact = 3.*(x(jnot(1))*x(index)**2 * - x(jnot(1))*x(index)**3) end if end if end subroutine sort (ier) c--------------------------------------------------------------------- c sort determines if a phase in the system defined in input may c be classified into a subsystem. i indexes the phase in the arrays. c c referenced by: main c references to: none c input arrays: cp,ibc,itc,iqac c output arrays: x,icct,iuct,ibct,itct,iqact,idu,idb c idt,idq,idqi c on return if ier = 1, the phase is not in the thermodynamic c composition space. c--------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) parameter (k1=1000,k6=7) parameter (l5=400,l6=400,l7=400,l8=200,l9=100) parameter (m1=21,m2=35,m3=35,m4=21,m5=7,h1=100,h2=100) character*8 names common/ cst3 /x(k1,k6)/ cst18 /icct(k1)/ cst12 /cp(k6,k1) * / cst6 /icomp,istct,i,icp/ cst8 /names(k1) * / cst14 /i2s,i3s,i4s,i5s,i6s,i7s common/ cst89 /i2c(2,21),i3c(3,35),i4c(4,35),i5c(5,21),i6c(6,7) * / cst91 /i1ct(k6),i2ct(m1),i3ct(m2),i4ct(m3),i5ct(m4), * i6ct(m5),i7ct/ cst19 /hi,iqth(4),ivchk(k1) common/ cst92 /id1(k6,l5),id2(m1,l6),id3(m2,l7),id4(m3,l8), * id5(m4,l9),id6(m5,h1),id7(h2) * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 c--------------------------------------------------------------------- tot = 0. ier = 0 icct(i) = 0 do 10 j = 1, icp if (cp(j,i).lt.1.d-8) cp(j,i) = 0.0 if (cp(j,i).eq.0.0) goto 10 icct(i) = icct(i) + 1 tot = tot + cp(j,i) 10 continue c if icct(i) = 0, the phase is c not in thermodynamic composition space if (icct(i).eq.0) then ier = 1 goto 99 end if c compute compositions: do 15 j = 1, icp 15 x(i,j) = cp(j,i)/tot c if icct