implicit double precision (a-g,o-y),integer (h-n,z) parameter (k1=100,h9=10) character*8 names,fname*10 common/ cst8 /names(k1)/ csta7 /fname(h9) * / cst6 /icomp,istct,iphct,icp * / cst79 /isoct * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 n9 = 10 open (10,file='dummy.dat') open (11,file='output') read (10,*) iphct do i = 1, iphct read (10,1000) names(i) end do read (10,*) isoct do i = 1, isoct read (10,1010) fname(i) end do 1010 format (a10) 1000 format (a8) call input9 end subroutine input9 c----------------------------------------------------------------------- c given a list of solution phase names (fname(h9)) solut searches a c data file (on unit n9) for the relevant data and subdivides the c solutions into pseudo-compounds. c----------------------------------------------------------------------- implicit double precision (a-g,o-z),integer (h-n) parameter (k1=100,k3=1200,k4=16,k5=12,k6=7,j9=600) character*10 tname, uname(2)*8 double precision z(3,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(4),jsitp(3),iosp(3), * jst(3),jfix,idfc(3),ist(3),isp(3),isite, * ispct(3,4),jnsp(3,4),jsite,iterm,iord common/ cst86 /xy(2,k1),y(k1,3,2),ntot,npairs common/ cst107 /a0(5,5),acoef(5,5,5),smult(5),ecoef(3,3,3), * nsite,nspm1(5),nterm(5,5), * nsub1(5,5,5,4),nsub2(5,5,5,4),nttyp(5,5,5) * / cst60 /iasmbl(j9),ipoint,imyn common/ cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst79 /isoct/ cst6 /icomp,istct,iphct,icp c----------------------------------------------------------------------- c set compound counter ipoint. ipoint = iphct c no request for solutions, goto 999: goto (999),io9 c initialize match flag: im = 0 c read the solution name 10 call rmodel (tname,istot) if (istot.eq.0) goto 99 c check the solution model: call cmodel (im,tname,istot,jstot,0,uname) if (jstot.eq.0) goto 10 c could revise im here. if (jstot.lt.istot) then c if ok, determine a valid c endmember subset. call emodel (tname,istot) c edit model parameters: if (istot.eq.0) goto 10 call nmodel end if call omodel (tname) c if nsite ne 0 get "normalization" c constants for entropy model: if (nsite.ne.0) call snorm (tname) c generate pseudocompound compositions. c subdiv returns the total c number of pseudocompounds (ipcps) and the c array y, of which the element y(h,i,j) is c the site fraction of the jth species on the c ith site of the hth pseudocompound. call subdiv (tname) c load default site fractions c into the unused sites of the c site fraction array: if (isite.lt.3) then do 20 i = isite + 1, 3 z(i,1) = 1.0 z(i,2) = 0.0 20 z(i,3) = 0.0 end if c subdiv generates ntot compositions, c generate the compound data for c each solution: do 30 h = 1, ntot c load the composition into c a the site fraction array: do 40 i = 1, isite zt = 0.0 do 50 j = 1, isp(i) - 1 z(i,j) = y(h,i,j) 50 zt = zt + z(i,j) 40 z(i,isp(i)) = 1. - zt c generate the pseudocompound: call soload (z,im,tname) 30 continue if (im.eq.isoct) goto 999 c next solution goto 10 c at least one solution phase referenced c in the input is not included in the c solution phase data file, write warning: 99 call warn (43,r,isoct-im,'INPUT9') 999 end 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=100,k3=1200,j9=600,h9=10) 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(4),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,5),smult(5),ecoef(3,3,3), * nsite,nspm1(5),nterm(5,5), * nsub1(5,5,5,4),nsub2(5,5,5,4),nttyp(5,5,5) common/ cst8 /names(k1)/ cst60 /iasmbl(j9),ipoint,imyn * / csta7 /fname(h9)/ cst18a / mname(3,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,*) isite call zeroi (isp,3,1,3) call zeroi (ist,3,1,3) 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,r,i,'PRIME9') c read the total number of terms, and the order c of the highest order term, in the 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,r,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,tname,istot,jstot,ibuild,xname) 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=100,k2=1200,k3=1200,j9=600,h9=10) character*8 mname,tname*10,fname*10,names,missin(36),xname(2) 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(4),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(3,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 disable fixed activity corrections: if (jfix.ne.0) then call warn (32,xfx(1),i,tname) goto 99 end if c if called by build (ibuild = 1) skip the c name check: if (ibuild.eq.1) goto 60 c check that the solution name matches a name c requested in the input from n1 do 50 i = 1, isoct if (tname.ne.fname(i)) goto 50 c got a match, goto 60: 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: 60 itl = 0 itd = 0 ith = 0 call zeroi (insp,12,1,12) call zeroi (jst,3,1,3) call zeroi (imsol,istot,0,36) call zeroi (ispct,12,0,12) 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 do 40 l = 1, 2 if (mname(i,j,k).eq.xname(l)) then write (*,1010) tname, xname(l) jstot = 0 goto 99 end if 40 continue 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,r,i,tname) end if if (ltyp(h).eq.10) then ith = ith + 1 if (ith.gt.3) call error (61,r,i,tname) else if (ltyp(h).gt.0) then itl = itl + 1 if (itl.gt.1) call error (61,r,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) ikp(h) = im 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,r,i,tname) else if (jstot.eq.1) then jstot = 0 call warn (24,r,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 ver112** the following endmembers', * ' are missing for solution ',a10,/,4(8(2x,a8),/)) 1010 format (' **warning ver113** ',a10,' is not a valid model', * ' because component ',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=100,k2=1200,k3=1200,j9=600,h9=10) 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(4),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 call zeroi (insp,9,0,9) call zeroi (iosp,3,0,3) 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 call zeroi (jsp,3,0,3) call zeroi (jnsp,9,1,9) 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 error (34,r,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,r,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 call zeroi (isp,3,1,3) 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 (*,*) write (*,*) ' iosp ',iosp 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** missing endmembers for solution: ', * a10,/,' the model will be recast with the 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=100,k2=1200,k3=1200) 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(4),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,5),smult(5),ecoef(3,3,3), * nsite,nspm1(5),nterm(5,5), * nsub1(5,5,5,4),nsub2(5,5,5,4),nttyp(5,5,5) 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 isub(jterm,i,1) = jsitp(isub(h,i,1)) 30 isub(jterm,i,2) = jnsp(isub(h,i,1),isub(h,i,2)) 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) character*10 tname common/ cst107 /a0(5,5),acoef(5,5,5),smult(5),ecoef(3,3,3), * nsite,nspm1(5),nterm(5,5), * nsub1(5,5,5,4),nsub2(5,5,5,4),nttyp(5,5,5) 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(4),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----------------------------------------------------------------- c read # of mixing sites if (nsite.gt.5) call error (31,r,i,tname) c if (nsite.lt.isite) call error (30,r,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) 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,r,nttyp(i,j,k),tname) end if 30 continue 20 continue 10 continue return 90 call error (20,r,i,tname) end subroutine snorm (tname) c-------------------------------------------------------------------- implicit double precision (a-g,o-z),integer (h-n) character*10 tname double precision x(3,3) common/ cst107 /a0(5,5),acoef(5,5,5),smult(5),ecoef(3,3,3), * nsite,nspm1(5),nterm(5,5), * nsub1(5,5,5,4),nsub2(5,5,5,4),nttyp(5,5,5) 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(4),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,ecoef(i,j,k), * i = 1, isp(1)),j = 1, isp(2)), k = 1, isp(3)) write (*,*) 2000 format (/,' the epsilon(ijk) coefficients (documentation ', * 'sect. 1.3.1, eq. 8b) for ',a10,' are:',/) 2010 format (10x,5(i1,i1,i1,' - ',f8.5,' ')) 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) character*10 tname double precision x(3,3) common/ cst107 /a0(5,5),acoef(5,5,5),smult(5),ecoef(3,3,3), * nsite,nspm1(5),nterm(5,5), * nsub1(5,5,5,4),nsub2(5,5,5,4),nttyp(5,5,5) 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 error (ier,real,int,char) implicit double precision (a-g,o-y),integer (h-m,z) character char*14 write (*,*) 'error ',ier,real,int,char stop end subroutine warn (ier,real,int,char) implicit double precision (a-g,o-y),integer (h-m,z) character char*14 write (*,*) 'warn ',ier,real,int,char end subroutine omodel (tname) c----------------------------------------------------------------------- c omodel - outputs revised solution models c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) parameter (k1=100,k3=1200,j9=600,h9=10) character*8 tname*10,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(4),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,5),smult(5),ecoef(3,3,3), * nsite,nspm1(5),nterm(5,5), * nsub1(5,5,5,4),nsub2(5,5,5,4),nttyp(5,5,5) common/ cst8 /names(k1) 10 write (11, 1000) tname c read number of independent sites: write (11,*) isite,' isite' c read number of species on each site: write (11,*) (isp(i),ist(i),i=1,isite),' isp, ist' c read endmember names: write (11,1010) (((names(imsol(i,j,k)), i= 1,isp(1)), * j= 1,isp(2)), k= 1,isp(3)) c read endmember flags: write (11, * ) ((( 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 20 write (11,*) (xmn(i,j), xmx(i,j), xnc(i,j), j= 1, isp(i)-1), * imd(i) c read the total number of terms, and the order c of the highest order term, in the mixing model. c ideal models: iterm = iord = 0. write (11,*) iterm, iord,' 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 write (11,*) ((isub(h,i,j), j = 1, 2), i = 1, iord) 30 write (11,*) (wg(h,k), k = 1, 11),' wg' c expansion for S(configurational) 40 write (11,*) nsite if (nsite.ne.0) call osmodel c permit fixed activity corrections: if (isite.eq.1) then write (11,*) jfix,' jfix' if (jfix.ne.0) write (11,*) (idfc(i), xfx(i),i = 1, jfix) end if 1000 format (a10) 1010 format (3(a8,1x)) 99 end subroutine osmodel c----------------------------------------------------------------- c osmodel 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) common/ cst107 /a0(5,5),acoef(5,5,5),smult(5),ecoef(3,3,3), * nsite,nspm1(5),nterm(5,5), * nsub1(5,5,5,4),nsub2(5,5,5,4),nttyp(5,5,5) 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(4),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 for each site do 10 i = 1, nsite c read # of species, and site c multiplicty. write (11,*) nspm1(i)+1,smult(i),' nsp, smult' 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. write (11,*) nterm(i,j), a0(i,j),' nterm, a0' c for each term: do 30 k = 1, nterm(i,j) c read term type: write (11,*) nttyp(i,j,k),' nttyp' if (nttyp(i,j,k).lt.5) then write (11,*) acoef(i,j,k), (nsub1(i,j,k,l), * nsub2(i,j,k,l), * l = 1, nttyp(i,j,k) ) end if 30 continue 20 continue 10 continue end subroutine subdiv (tname) c---------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n) character*10 tname parameter (k1=100) 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(4),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 one and two c | \ dimensional simplexes. iscp is the number of c | \ vertices, in the one-dimensional case the c xinc2 | \ ordinates of the vertices are xmin1 and xmax1, c | \ in the two dimensional case the coordinates c | \ are (xmin1,xmin2), (xmax1, xmin2) and (xmin1, c xmin2 | \ xmax2). c --------------- c xmin1 xinc1 xmax1 c---------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) parameter (k1=100) 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 = ((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,i,sname) bm1 = b - 1. bp1 = b + 1. dx = 1. / (xinc1 - 1.) x = dx loop = 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 else if (imod.eq.2) then c assymmetric stretching imod = 2: b = xmin1 if (b.lt.1.d0) call error (207,b,i,sname) bm1 = b - 1. bp1 = b + 1. dx = 1. / (xinc1 - 1.) x = dx loop = 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,real,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, real, int, 'CART') xmax1 = 1. - xmin2 end if if (xmax2.gt.(1.-xmin1)) then call warn (171, real, int, 'CART') xmax2 = 1. - xmin1 end if if (imod.eq.0) then c cartesian subdivision: loop1 = ((xmax1-xmin1)/xinc1 + 1.0000001) step = 0.0 do 10 i = 1, loop1 loop2 = ((xmax2-xmin2-step)/xinc2 + 1.000001) y = xmin1 + step x = xmin2 do 20 j = 1, loop2 c generate compositions: npairs = npairs+1 xy(1,npairs) = x xy(2,npairs) = y 20 x = x + xinc2 10 step = step + xinc1 else if (imod.gt.0.and.imod.lt.5) then c barycentric subdivision of ternary c solutions, imod = # of iterations. npairs = ibary(imod) 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 value of imod is illegal, write error msg: call error (169,real,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,int,tname) xmax = 1. end if if (xmin.lt.0.0) then call warn (22,xmin,int,tname) xmin = 0. end if if (xmax.lt.xmin) then call warn (23,xmax,int,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.ge.0.0.and.z.le.1.0) then return else if (z.gt.1.0.and.z.lt.1.00001) then z = 1.0 else if (z.gt.-.000001) then z = 0.0 else call error (125,z,i,tname) end if end subroutine soload (z,isolct,tname) implicit double precision (a-g,o-z),integer (h-n) parameter (k1=100,k2=1200,k3=1200,k4=16,k5=12,k6=7) parameter (m9=50,k9=50,m6=3,m7=9,m8=9,h5=5,h6=40,j9=600) double precision z(3,3) character*10 tname, 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(4),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,5),smult(5),ecoef(3,3,3), * nsite,nspm1(5),nterm(5,5), * nsub1(5,5,5,4),nsub2(5,5,5,4),nttyp(5,5,5) common/ cst1 /thermo(k4,k1),vf(2,k1),vs(h5,k1),uf(2),us(h5) * / cst203 /therdi(m8,m9),therlm(m7,m6,k9) * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn * / cst3 /x(k1,k6)/ 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 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,r,i,'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(1d2*z(1,1)) else if (z(1,1).gt.0.96) then write (names(iphct),1070) names(imsol(1,1,1)), * 1d2*z(1,1) else write (names(iphct),1060) names(imsol(1,1,1)), * 1d2*z(1,1) end if else c ternary solutions write (names(iphct),1040) names(imsol(1,1,1)), * idint(z(1,1)*100.),names(imsol(2,1,1)),idint(z(1,2)*100.) end if else if (isite.eq.2.and.isp(1).eq.2.and.isp(2).eq.2) then int1 = 100.*z(1,1) int2 = 100.*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 have live c with it. write (names(iphct),1010) ((idint(10.*z(i,j)), * j = 1, isp(i) - 1), i = 1, isite) write (n3,1000) (((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 initialize constants: smix = 0.0 esum = 0.0 do 130 i = 1, icp x(iphct,i) = 0.0 130 cp(i,iphct) = 0.0 do 140 i = 1, k4 140 thermo(i,iphct) = 0.0 do 150 i = 1, icomp 150 comps(iphct,i) = 0.0 c load constants: do 160 i = 1, isp(1) if (z(1,i).eq.0.0) goto 160 do 170 j = 1, isp(2) if (z(2,j).eq.0.0) goto 170 do 180 k = 1, isp(3) if (z(3,k).eq.0.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) do 190 l = 1, icp x(iphct,l) = x(iphct,l) + zp * x(id,l) 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) c 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 c call sort c compute ideal configurational negentropy: if (nsite.eq.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 (ist(1).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 write (*,*) names(iphct),smix 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) 1060 format (a2,f3.1) 1070 format (a2,f4.1) 999 end