c----------------------------------------------------------------------- c RLIB - a library of subprograms called by VERTEX, FRENDLY, ACTCOR c BUILD, ADIABAT, and ISOCHOR. 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----------------------------------------------------------------------- function gexces (id) c----------------------------------------------------------------------- c gexces evaluates the contribution to the gibbs energy of a pseudocompound c arising from configurational entropy and excess properties. c----------------------------------------------------------------------- implicit double precision (a-g,o-z),integer (h-n) parameter (k1=5000,k12=14,k13=40000,h9=30,k5=12) common/ cst304 /sxs(k13),exces(3,k1),jend(h9,k12),ixp(k1),ifp(k1) * / cst5 /p,t,xco2,u1,u2,tr,pr,r,ps * / cst11 /f(2)/ cst12 /cp(k5,k1) c----------------------------------------------------------------------- c two cases: c normal solution: if (ifp(id).eq.0) then gexces = exces(1,id) + t * exces(2,id) + p * exces(3,id) else if (pr.eq.0d0.or.tr.eq.0d0.or. * (cp(1,id).eq.0d0.and.cp(2,id).eq.0d0)) then gexces = 0d0 c if pr or tr = 0 then c this is a special calculation c and we can skip this part. else c modification for fluid c projections xco2 = cp(2,id) call cfluid (fo2,fs2) gexces = r*t*(cp(1,id)*f(1) + cp(2,id)*f(2)) end if end function gzero (id) c----------------------------------------------------------------------- c gzero computes the reference pressure free energy of a compound c with no transitions identified by the arguement 'id'. c see gcpd. c--------------------------------------------------------------------- implicit double precision (a-g,o-z),integer (h-n) parameter (k1=5000,k4=22,h5=7,k10=200) common/ cst1 /thermo(k4,k10),vf(2,k1),vs(h5,k1),uf(2),us(h5) * / cst88 /vu(2,k1),jprct,jmct * / cst5 /p,t,xco2,u1,u2,tr,pr,r,ps/ cst44 /dlogt,dsqrtt gzero = 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 -ndu * - vu(1,id) * u1 - vu(2,id) * u2 end function gphase (id) c----------------------------------------------------------------------- c gphase computes the gibbs free energy of a compound identified by index c id, it is assumed that the variables in cst44 are set prior to the call. c gphase does not assume that the properties of a pseudocompound endmembers c are known (unlike gall) it is thus less efficient than routine gall. c----------------------------------------------------------------------- implicit double precision (a-g,o-z),integer (h-n) parameter (k1=5000,k12=14,k13=40000,h9=30,k2=3000) common/ cst304 /sxs(k13),exces(3,k1),jend(h9,k12),ixp(k1),ifp(k1) common/ cst60 /ipoint,imyn * / cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar if (id.le.ipoint) then c the phase is an endmember compound graw = gcpd(id) else c the phase is a pseudocompound, first c compute endmember properties: graw = 0d0 if (ifp(id).eq.0) then c excess props don't include vdp: do 10 k = 1, jend(ikp(id),1) 10 graw = graw + gcpd(jend(ikp(id),2+k))*sxs(ixp(id)+k) else do 20 k = 1, jend(ikp(id),1) 20 graw = graw + gzero(jend(ikp(id),2+k))*sxs(ixp(id)+k) end if c now get the excess and ideal mixing effect: graw = graw + gexces(id) end if gphase = graw end function gcpd (id) c----------------------------------------------------------------------- c gcpd computes the gibbs free energy of a compound 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(4-10,id)], and the coefficients of the volumetric c equations given in the program documentation (Eqs 2.1-2.3): c b8 = 0 => c v(p,t) = v(pr,tr) + b2*(t-tr) + b4*(p-pr) c + b6*(p-pr)**2 + b7*(t-tr)**2 c b8 < 0 => c v(p,t) = v(pr,tr) * exp[ b3*(t-tr) + b8*(p-pr) ] c b8 > 0 => c v(p,t) = v(pr,t) * [1 + b8*p/(b8*p+Kt)]^(1/b8) c Kt = b6 + b7 * t c alpha = b1 + 2*b2*t + b3/t - b4/t^2 + b5/t^(1/2) c [thermo(11-18,id)]. c the variables in cst44 must be set before a call to gcpd. c--------------------------------------------------------------------- implicit double precision (a-g,o-z),integer (h-n) parameter (k1=5000,k4=22,k5=12,m9=10,m6=3,m7=12,k9=25,m8=9, * h5=7,k10=200) common/ cst1 /thermo(k4,k10),vf(2,k1),vs(h5,k1),uf(2),us(h5) * / cst88 /vu(2,k1),jprct,jmct/ cst11 /f(2) * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn/ cst113 /ifproj * / cst5 /p,t,xco2,u1,u2,tr,pr,r,ps/ cst12 /cp(k5,k1) * / cst203 /therdi(m8,m9),therlm(m7,m6,k9) * / cst204 /ltyp(k1),lmda(k1),idis(k1)/ cst44 /dlogt,dsqrtt save rkt data rkt/0d0/ graw = 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-ndu term: if (thermo(18,id).eq.0d0) then c normal polynomial: vdp = p * (thermo(3,id) * + (thermo(17,id) * t + thermo(12,id)) * t * + (thermo(16,id) * p + thermo(14,id)) * p) c -ndu * - vu(1,id) * u1 - vu(2,id) * u2 else if (thermo(18,id).gt.0d0) then c murghnahan EoS: c vdp = V(1,T)*KT**(1/K')/(K'-1) c * [(KT+K'*p)**(1-1/K') - c (KT+K'*pr)**(1-1/K')] rkt = thermo(16,id) + thermo(17,id)*t vdp = (thermo(3,id) + (thermo(11,id) + thermo(12,id)*t)*t + * thermo(13,id)*dlogt + thermo(14,id)/t + * thermo(15,id)*dsqrtt) * rkt**thermo(21,id)/thermo(22,id) * *((rkt+thermo(18,id)*p)**thermo(19,id) * -(rkt+thermo(20,id))**thermo(19,id)) * - vu(1,id) * u1 - vu(2,id) * u2 else c gottschalk. vdp = thermo(11,id)*dexp(thermo(13,id)*t)* * (1d0 - dexp(thermo(18,id)*(p - pr))) * - vu(1,id) * u1 - vu(2,id) * u2 end if graw = graw + vdp c check for lambda transitions: if (ltyp(id).ne.0) then if (ltyp(id).lt.4) then c ubc-type transitions call lamubc (p,t,dg,lmda(id),ltyp(id)) graw = graw + dg else if (ltyp(id).lt.7) then c standard transitions call lamhel (p,t,graw,vdp,lmda(id),ltyp(id)) else if (ltyp(id).lt.10) then c supcrt q/coe lambda transition du = - vu(1,id) * u1 - vu(2,id) * u2 call lamqtz (p,t,graw,du,lmda(id),id) else if (ltyp(id).gt.9 .and. ltyp(id).le.12 ) then c holland and powell (landau model) call lamlan (dg,rkt,lmda(id),ltyp(id)) graw = graw + 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) graw = graw + gss end if c modification for fluid c projections: if (ifproj.ne.0.and.thermo(3,id).eq.0d0) then c if pr or tr = 0 then c this is a special calculation c and we can skip this part. if (pr.eq.0d0.or.tr.eq.0d0) goto 999 if (cp(1,id).eq.0d0.and.cp(2,id).eq.0d0) goto 999 if (ifproj.eq.1) then xco2 = 0d0 call cfluid (fo2,fs2) graw = graw + r*t*cp(1,id)*f(1) else if (ifproj.eq.2) then xco2 = 1d0 call cfluid (fo2,fs2) graw = graw + r*t*cp(1,id)*f(2) else xco2 = cp(2,id) call cfluid (fo2,fs2) graw = graw + r*t*(cp(1,id)*f(1) + cp(2,id)*f(2)) end if end if 999 gcpd = graw end subroutine rmodel (tname,istot) c--------------------------------------------------------------------- c rmodel - reads solution models from LUN n9. c--------------------------------------------------------------------- implicit double precision (a-g,o-z),integer (h-n) parameter (k1=5000,h9=30,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) * / cst8 /names(k1) * / 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-z),integer (h-n) parameter (k1=5000,k2=3000,h9=30) 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) * / cst8 /names(k1)/ cst60 /ipoint,imyn/ cst79 /isoct * / cst204 /ltyp(k1),lmda(k1),idis(k1)/ csta7 /fname(h9) * / cst61 /ikp(k1),ivarrx(k2),ivarip(k2),isudo,ivar 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 i = 1, 3 jst(i) = 1 do j = 1, 4 ispct(i,j) = 0 do k = 1, 3 imsol(j,i,k) = 0 end do insp(i,j) = 1 end do end do do 70 h = 1, ipoint do i = 1, isp(1) do 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 jstot = jstot + 1 imsol(i,j,k) = h c do i1 = 1, isite do j1 = 1, isp(i1) if (ijk(i1).eq.j1) ispct(i1,j1)=ispct(i1,j1)+1 end do end do c set kill flag: if (iend(i,j,k).eq.2) ikp(h) = -1 c found all endmembers: if (jstot.eq.istot) goto 99 goto 70 100 continue end do end do 70 continue c missing endmember warnings: if (jstot.lt.2) then im = im - 1 call warn (25,xfx(1),jstot,tname) jstot = 0 else imiss = 0 do i = 1, isp(1) do 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 end do end do 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-z),integer (h-n) parameter (k1=5000,k2=3000) 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(k2),isudo,ivar do i = 1, 3 iosp(i) = 0 do j = 1, 4 insp(i,j) = 0 end do end do 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 j = 1, isp(1) if (ispct(1,j).gt.0) insp(1,j) = 1 end do end if c count the number of endmembers c left: 600 do i = 1, 3 jsp(i) = 0 do j = 1, 4 jnsp(i,j) = 1 end do end do c and make a set of valid species do 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 end do jstot = 0 do i = 1, isp(1) do 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 end do end do kstot = 1 do i = 1, 3 if (jsp(i).eq.0) jsp(i) = 1 kstot = kstot * jsp(i) jsitp(i) = i end do 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 i = 1, isite do 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 end do end do insp(ik,jk) = 0 goto 600 end if c make a pointer for old site c index to new: jsite = 0 jmax = isite do 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 end do 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' * ,1x,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-z),integer (h-n) parameter (k1=5000,k2=3000,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(k2),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) = 0d0 smult(lsite) = smult(i) else if (l3.eq.0) then if (l4.eq.0) then lsp = lsp + 1 a0(lsite,lsp) = 0d0 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 = 0d0 c for each site do 10 i = 1, nsite zt = 0d0 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.0d0.and.z.lt.1d0) then dlnw = dlnw - smult(i) * z * dlog (z) end if 20 continue z = 1d0 - zt call zchk (z,tname) if (z.gt.0d0.and.z.lt.1d0) * dlnw = dlnw - smult(i) * z * dlog (z) 10 continue end subroutine subdiv (tname) c--------------------------------------------------------------------- implicit double precision (a-g,o-z),integer (h-n) character*10 tname parameter (k1=5000) 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-z),integer (h-n) parameter (k1=5000) 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(i), i = 1, 3)/ 0., 1., 0./ data (vco2(i), i = 1, 3)/ 1., 0., 0./ data (vco1(i), i = 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(i), i = 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(i), i = 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(i), i = 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(i), i = 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(i), i = 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(i), i = 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(i), i = 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.0d0) 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 if (xmax1.eq.0d0.or.dabs(xmax1).gt.1d0) then delx = 1d0 x0 = 0d0 else if (xmax1.gt.0d0) then delx = xmax1 x0 = 0d0 else if (xmax1.lt.0d0) then delx = dabs(xmax1) x0 = 1d0 - delx end if c symmetric stretching imod = 1: b = xmin1 if (b.lt.1d0) 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) = x0 + (0.5 * (bp1-bm1*t1)/(1.+t1))*delx 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) = 2d0*x0 + delx - xy(1,jpairs) 60 jpairs = jpairs - 1 else if (imod.eq.2) then c assymmetric stretching imod = 2: b = xmin1 if (b.lt.1d0) 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.1d0) 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-z),integer (h-n) character tname*10 if (xmax.gt.1d0) then call warn (21,xmax,1,tname) xmax = 1d0 end if if (xmin.lt.0d0) then call warn (22,xmin,1,tname) xmin = 0d0 end if if (xmax.lt.xmin) then call warn (23,xmax,1,tname) xmax = 1d0 xmin = 0d0 end if end subroutine zchk (z, tname) implicit double precision (a-g,o-z),integer (h-n) character tname*10 if (z.gt.0d0.and.z.le.1d-3) then z = 0d0 else if (z.ge.0.999.and.z.lt.1d0) then z = 1d0 else if (z.ge.0d0.and.z.le.1d0) then return else if (z.gt.1d0.and.z.lt.1.01) then z = 1d0 else if (z.gt.-.01) then z = 0d0 else call error (125,z,1,tname) end if 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-z),integer (h-n) parameter (l2=5,k6=1000) common/ cst32 /ptx(k6),ipt2 * / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5 * / cst5 /v(l2),tr,pr,r,ps/ cst63 /delv(l2) ipt2 = ipt2 + 2 if (ipt2.gt.k6) ipt2 = k6 ptx(ipt2-1) = v(iv1) ptx(ipt2) = v(iv2) 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-z),integer (h-n) parameter (l2=5) common/ cst62 /blim(l2),ulim(l2),dgr * / cst9 /vmax(l2),vmin(l2),dv(l2) c--------------------------------------------------------------------- do i = 1, 5 if (dv(i).lt.0d0) call error (34,dv(i),i,'CONCRT') ulim(i) = vmax(i)+dv(i) blim(i) = vmin(i)-dv(i) ddv = vmax(i)-vmin(i) if (ddv.lt.0d0) call error (35,ddv,i,'CONCRT') end do end subroutine conver (g,s,v,a,b,c,d,e,f,gg,b1,b2,b3,b4,b5,b6,b7,b8, * b9,b10,b11,b12,tr,pr) 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-z),integer (h-n) c b8 is K' in the murgnahan EoS c or -beta in exponential volumetric c functions. if (b8.eq.0d0) then c normal polynomial vdp term: c first compute G* g = g c add the integral of sdt at tr,pr: * + s * tr - a * tr - b * tr * tr / 2d0 + c / tr * - d * tr**3 / 3.d0 - 2d0 * e * dsqrt(tr) * - f * dlog(tr) + gg / tr / tr / 2d0 + f c add the constant components of the c integral -vdp at (t,pr): * - v * pr + b2 * tr * pr + b4 * pr * pr / 2d0 * - 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/2d0 * + d * tr * tr / 2d0 - 2d0 * e / dsqrt(tr) - f/tr * - gg / tr**3 / 3.d0 c V* (P) v = v - b2 * tr - b4 * pr - b7 * tr * tr c a* (-T log T) c a = a c b* (-T*T) b = b7 * pr + b / 2d0 c c* (-1/T) c = c / 2d0 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/2d0 c b5* (exp (P/ c2)) c b5 = c2 * b5 c b6* (P**3) b6 = b6/3.d0 c b7* (P*T**2) c b7 = b7 else c nonlinear EoS's: c thermal terms: g = g c add the integral of sdt at tr,pr: * + s * tr - a * tr - b * tr * tr / 2d0 + c / tr * - d * tr**3 / 3d0 - 2d0 * e * dsqrt(tr) * - f * dlog(tr) + gg / tr / tr / 2d0 + f c S* (T) s = a - s + a * dlog(tr) + b * tr - c / tr / tr /2d0 * + d * tr * tr / 2d0 - 2d0 * e / dsqrt(tr) - f / tr * - gg / tr**3 / 3d0 c these terms are as above b = b / 2d0 c = c / 2d0 d = d / 6.d0 e = 4d0 * e gg = gg/6.d0 c volume terms: if (b8.gt.0d0) then c mughnahan vdp b1 = v*b1 b2 = v*b2/2d0 b3 = v*b3 b4 = -v*b4 b5 = 2d0*v*b5 v = v - (b1-b2*tr)*tr - b3*dlog(tr) - b4/tr - b5*dsqrt(tr) b6 = b6 b7 = b7 c b8 = kappap c these are unnecessary, but may cut time? b9 = (1d0-1d0/b8) b10 = b8*pr b11 = 1d0/b8 b12 = b8 - 1d0 else c gottschalk, b3 = alpha, -b8 = beta: b1 = -v / b8 / dexp (b3*tr) end if end if 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-n) common/ cst5 /p,v(4),pr,z(3) g = 0d0 if (tk .lt. tr) return tt = tk if (tk.gt.tmax) tt = tmax dh = d0 * (tt - tr) * + 2d0 * d1 * (tt**0.5 - tr**0.5) * - d2 * (1d0 / tt - 1d0 / tr) * + d4 * dlog(tt/tr) * + d5 * (tt*tt - tr*tr) / 2d0 * + d6 * (tt**3 - tr**3) / 3d0 ds = d0 * dlog(tt/tr) * - 2d0 * d1 * (tt**(-0.5) - tr**(-0.5)) * - d2 * (1d0/tt/tt - 1d0/tr/tr) / 2d0 * - d4 * (1d0/tt - 1d0 / tr) * + d5 * (tt - tr) * + d6 * (tt*tt - tr*tr) / 2d0 gss = dh - (tk * ds) if (d3.ne.0d0) gss = gss + dh/d3 * (p - pr) end subroutine getphi (name,ibase,icmpn,*) c----------------------------------------------------------------------- implicit double precision (a-g,o-z),integer (h-n) parameter (k0=25,k4=22,k5=12,m6=3,m7=12,m8=9) character*8 name, oname, record*1 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,1010,end=90,err=98) record if (record.eq.' ') then c check for comments, i.e., data c with a blank 1st character goto 30 else backspace (n2) read (n2,1020,end=90,err=98) * name, ibase, ikind, ilam, idiso end if read (n2,*,err=98) (comp(i), i=1, icmpn), (therm(i),i=1,18) 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.0d0) 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 - 2) 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 1010 format (a1) 1020 format (a8,i2,i2,i2,i2) end subroutine chkphi (ichk,name,icout,ibase,icmpn,*) c----------------------------------------------------------------------- implicit double precision (a-g,o-z),integer (h-n) parameter (k0=25,k4=22,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 c reject the data if excluded in input: do i = 1, ixct if (name.eq.exname(i)) goto 90 end do 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 i= 1, icmpn if (icout(i).eq.0.and.comp(i).ne.0d0) goto 90 end do c reject phases which consist entirely c of constrained components: if (iprct.eq.0) goto 99 do i= 1, icp if (comp(ic(i)).ne.0d0) goto 99 end do 90 return 1 99 end subroutine lamlan (dg,rkt,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 p = pressure in bar c ld = pointer to phase in therlm c rkt = bulk modulus c returned - dg - difference in free energy due to the transition c--------------------------------------------------------------------- implicit double precision (a-g,o-z),integer (h-n) parameter (m9=10,k9=25,m6=3,m7=12,m8=9,k1=5000,k4=22, * h5=7,k10=200) common/ cst203 /therdi(m8,m9),therlm(m7,m6,k9) * / cst5 /p,t,xco2,u1,u2,tr,pr,r,ps * / cst1 /thermo(k4,k10),vf(2,k1),vs(h5,k1),uf(2),us(h5) lct = ltype - 9 c don't do this stuff until TJBH c has a consistent policy: c create a constant which contains c the 0 to Tr contribution to the c free energy: c q = (1.-(tr/tc))**0.25 c if (tr.gt.tc) q = 0. c maximum enthalpy of disorder: c therlm (3,j,lamin) = 2./3.*smax*tc c enthalpy of disorder at tr: c therlm (4,j,lamin) = 2.*smax*tc * c * (q**6/6. - q*q/2. + 1./3.) c entropy of disorder at tr: c therlm (5,j,lamin) = dg = 0d0 do 10 i = 1, lct tc = therlm(1,i,ld) + therlm(3,i,ld)*(p-pr) smax = therlm (2,i,ld) if (t.gt.tc) then c fully disordered: dg = dg + therlm(4,i,ld) - t*therlm(5,i,ld) else c partially disordered: q2 = dsqrt(1d0-t/tc) c bulk modulus: dg = dg c hr-t*sr * + therlm(4,i,ld) - t*therlm(5,i,ld) c G(landau) * + smax * (q2*(t-tc) + tc*q2**3/3.d0) end if c if pressure dependent c add the vdp term: c pressure dependent c landau transitions are c only in the mughnahan c versions of HP data. if (therlm(3,i,ld).ne.0d0) then dg = dg c vt*kt/3 * + (therlm(6,i,ld) + therlm(7,i,ld)*t * + therlm(8,i,ld)*dsqrt(t))*rkt/3.d0 c *[(1+4p/kt)^(3/4) - 1] c this approximation is only c good for pr near 0. * *((1d0 + 4.d0*p/rkt)**(0.75d0) - 1d0) end if 10 continue end function gtrans (id,j) c----------------------------------------------------------------------- c gtrans computes the reference pressure free energy of a compound c aboves its jth transition. c--------------------------------------------------------------------- implicit double precision (a-g,o-z),integer (h-n) parameter (m9=10,k9=25,m6=3,m7=12,m8=9) common/ cst203 /therdi(m8,m9),therlm(m7,m6,k9) * / cst5 /p,t,xco2,u1,u2,tr,pr,r,ps/ cst44 /dlogt,dsqrtt gtrans = therlm(12,j,id) c -sdt * + t * (therlm(3,j,id) - therlm(5,j,id) * dlogt * - t * (therlm(6,j,id) + therlm(8,j,id) * t)) * - (therlm(7,j,id) + therlm(11,j,id) / t) / t * + therlm(9,j,id) * dsqrtt + therlm(10,j,id)*dlogt end subroutine lamhel (p,t,g,vdp,ld,ltype) 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-n) parameter (m9=10,k9=25,m6=3,m7=12,m8=9) common/ cst203 /therdi(m8,m9),therlm(m7,m6,k9)/ cst5 /v(6),pr,u(2) lct = ltype - 3 c T