c Please do not distribute any part of this source. 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 VLIB - a library of subprograms called by VERTEX. subroutine abload (*) c--------------------------------------------------------------------- c abload assembles the matrix 'a' and vector 'b' for hcp component c sytems and then solves the equation ax = b, the vector x is c returned in 'b'. c referenced by: simpl2,simpl3,simpl4,simpl5 c references to: subst,factor c---------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) parameter (k1=1200,k6=7) common/ cst23 /a(10,10),b(11),ipvt(10),idv(10),iophi,idphi, * iiphi,iflg1 * / cst2 /g(k1)/ cst12 /cp(k6,k1)/ cst52 /hcp,h,id(8),ic(7) c----------------------------------------------------------------------- do 10 i = 1, hcp do 20 j = 1, hcp 20 a(i,j) = cp(ic(j),id(i)) 10 continue call factor (a,hcp,ipvt,ier) goto (99),ier do 30 i = 1, hcp 30 b(i) = g(id(i)) call subst (a,ipvt,hcp,b) return 99 return 1 end subroutine asschk c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) parameter (k1=1200,k6=7,l2=5) common/ cst23 /a(10,10),b(11),ipvt(10),idv(10),iophi,idphi, * iiphi,iflg1/ cst87 /delt(l2),dtol,utol,ptol * / cst6 /icomp,istct,iphct,icp/ cst62 /blim(l2),ulim(l2),dgr * / cst2 /g(k1)/ cst12 /cp(k6,k1)/ cst7 /iflag c----------------------------------------------------------------------- c determine chemical potentials iflag = 0 do 10 i = 1, icp 10 b(i) = g(idv(i)) c call subst (a,ipvt,icp,b) c test phases against the c assemblage idv do 20 i = istct, iphct gphi = 0. do 30 j = 1, icp 30 gphi = gphi + cp(j,i)*b(j) dg = g(i) - gphi if (dg.gt.dtol) goto 20 c check that a phase is not metastable c with respect to itself, this c could be remedied by changing the c value of dtol. do 25 j = 1, icp 25 if (idv(j).eq.i) goto 20 c assemblage is metastable with respect c to phase idphi. iflag = iflag+1 idphi = i if (dabs(dg).gt.utol) iflag = iflag + 1 goto (20,40,40),iflag 20 continue 40 end block data c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) parameter (h5=7,k1=1200,k4=16,k6=7,l5=400,l6=400,m1=21) common/ cst5 /p,t,xco2,u1,u2,tr,pr,r,ps/ cst85 /pp,tt,yy,rr * / cst13 /id1(k6,l5),id2(m1,l6)/ cst134 /iperm(6,7,7) * / cst19 /hi,iqth(4),ivchk(k1)/ cst32 /ptx(450),logx,ipt2 * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn * / cst1 /thermo(k4,k1),vf(2,k1),vs(h5,k1),uf(2),us(h5) * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst89 /i2c(2,21),i3c(3,35),i4c(4,35),i5c(5,21),i6c(6,7) * / cst90 /i23(3,35),i34(4,35),i45(5,21),i56(6,7) * / cst93 /i24a(420),i24(6,35)/ cst64 /i7c(7,8),i8c(8,9) * / cst94 /i25a(420),i25(10,21),i35a(630),i35(10,21) * / cst95 /i16a(6,7),i26a(5,6,7),i26(15,7),i36a(420), * i36(20,7),i46a(420),i46(15,7) * / cst101 /i27a(6,7),i37a(15,7),i47a(20,7),i57a(15,7) c----------------------------------------------------------------------- data iff/2*0/,ipt2,logx/0,1/ c data us, uf/ h5*0.d0, 2*0.d0/ c data r,rr/8.3144126d0,83.14d0/ c the following data statement assigns c logical unit numbers for i/o. data n1,n2,n3,n4,n5,n6,n7,n8,n9/21,22,23,24,25,26,21, 6,29/ c------------------------------------------------------------------- c SUBSYSTEM COMPONENT ARRAYS: c i2c(i,j) is an array identifying the c i components in the jth binary. data i2c/ 1,2, 1,3, 2,3, 1,4, 2,4, 3,4, 1,5, 2,5, 3,5, 4,5, * 1,6, 2,6, 3,6, 4,6, 5,6, 1, 7, 2,7, 3,7, 4,7, 5,7, * 6,7 /, c i3c(i,j) is an array identifying the c i components in the jth ternary. * i3c/ 1,2,3, 1,2,4, 1,3,4, 2,3,4, 1,2,5, 1,3,5, 2,3,5, * 1,4,5, 2,4,5, 3,4,5, 1,2,6, 1,3,6, 2,3,6, 1,4,6, * 2,4,6, 3,4,6, 1,5,6, 2,5,6, 3,5,6, 4,5,6, 1,2,7, * 1,3,7, 2,3,7, 1,4,7, 2,4,7, 3,4,7, 1,5,7, 2,5,7, * 3,5,7, 4,5,7, 1,6,7, 2,6,7, 3,6,7, 4,6,7, 5,6,7/ c i4c(i,j) is an array identifying the c i components in the jth quaternary. data i4c/1,2,3,4, 1,2,3,5, 1,2,4,5, 1,3,4,5, 2,3,4,5, * 1,2,3,6, 1,2,4,6, 1,3,4,6, 2,3,4,6, 1,2,5,6, * 1,3,5,6, 2,3,5,6, 1,4,5,6, 2,4,5,6, 3,4,5,6, * 1,2,3,7, 1,2,4,7, 1,3,4,7, 2,3,4,7, 1,2,5,7, * 1,3,5,7, 2,3,5,7, 1,4,5,7, 2,4,5,7, 3,4,5,7, * 1,2,6,7, 1,3,6,7, 2,3,6,7, 1,4,6,7, 2,4,6,7, * 3,4,6,7, 1,5,6,7, 2,5,6,7, 3,5,6,7, 4,5,6,7/ c i5c(i,j) is an array identifying the c i components in the jth quinary. data i5c/1,2,3,4,5, 1,2,3,4,6, 1,2,3,5,6, 1,2,4,5,6, * 1,3,4,5,6, 2,3,4,5,6, 1,2,3,4,7, 1,2,3,5,7, * 1,2,4,5,7, 1,3,4,5,7, 2,3,4,5,7, 1,2,3,6,7, * 1,2,4,6,7, 1,3,4,6,7, 2,3,4,6,7, 1,2,5,6,7, * 1,3,5,6,7, 2,3,5,6,7, 1,4,5,6,7, 2,4,5,6,7, * 3,4,5,6,7/ c i6c(i,j) is an array identifying the c i components in the jth hexary. data i6c/1,2,3,4,5,6, 1,2,3,4,5,7, 1,2,3,4,6,7, 1,2,3,5,6,7, * 1,2,4,5,6,7, 1,3,4,5,6,7, 2,3,4,5,6,7/ c------------------------------------------------------------------- c i7c and i8c are arrays used for enumerating univariant fields c around invariant points in 6 and 7 component systems, the arrays c i3c through i6c are used for this purpose in simpler systems. c data i7c/1,2,3,4,5,6,7, 1,2,3,4,5,6,8, 1,2,3,4,5,7,8, * 1,2,3,4,6,7,8, 1,2,3,5,6,7,8, 1,2,4,5,6,7,8, * 1,3,4,5,6,7,8, 2,3,4,5,6,7,8/ data i8c/1,2,3,4,5,6,7,8, 1,2,3,4,5,6,7,9, 1,2,3,4,5,6,8,9, * 1,2,3,4,5,7,8,9, 1,2,3,4,6,7,8,9, 1,2,3,5,6,7,8,9, * 1,2,4,5,6,7,8,9, 1,3,4,5,6,7,8,9, 2,3,4,5,6,7,8,9/ c------------------------------------------------------------------- c component permutations iperm(i,j,k) c iperm(i(i=1,k),j(j=1,k),k) contains c the jth permutation of the k-1 components c in a system of k components. data iperm/ * 0,0,0,0,0,0, 0,0,0,0,0,0, 0,0,0,0,0,0, 0,0,0,0,0,0, 0,0,0,0,0,0, * 0,0,0,0,0,0, 0,0,0,0,0,0, 2,0,0,0,0,0, 1,0,0,0,0,0, 0,0,0,0,0,0, * 0,0,0,0,0,0, 0,0,0,0,0,0, 0,0,0,0,0,0, 0,0,0,0,0,0, 2,3,0,0,0,0, * 1,3,0,0,0,0, 1,2,0,0,0,0, 0,0,0,0,0,0, 0,0,0,0,0,0, 0,0,0,0,0,0, * 0,0,0,0,0,0, 2,3,4,0,0,0, 1,3,4,0,0,0, 1,2,4,0,0,0, 1,2,3,0,0,0, * 0,0,0,0,0,0, 0,0,0,0,0,0, 0,0,0,0,0,0, 2,3,4,5,0,0, 1,3,4,5,0,0, * 1,2,4,5,0,0, 1,2,3,5,0,0, 1,2,3,4,0,0, 0,0,0,0,0,0, 0,0,0,0,0,0, * 2,3,4,5,6,0, 1,3,4,5,6,0, 1,2,4,5,6,0, 1,2,3,5,6,0, 1,2,3,4,6,0, * 1,2,3,4,5,0, 0,0,0,0,0,0, 2,3,4,5,6,7, 1,3,4,5,6,7, 1,2,4,5,6,7, * 1,2,3,5,6,7, 1,2,3,4,6,7, 1,2,3,4,5,7, 1,2,3,4,5,6/ c------------------------------------------------------------------- c SUBSYSTEM ARRAYS: c c i23(i,j) is an array identifying the c ith binary in the jth ternary. data i23/ * 1, 2, 3, 1, 4, 5, 2, 4, 6, 3, 5, 6, 1, 7, 8, 2, 7, 9, * 3, 8, 9, 4, 7,10, 5, 8,10, 6, 9,10, 1,11,12, 2,11,13, * 3,12,13, 4,11,14, 5,12,14, 6,13,14, 7,11,15, 8,12,15, * 9,13,15, 10,14,15, 1,16,17, 2,16,18, 3,17,18, 4,16,19, * 5,17,19, 6,18,19, 7,16,20, 8,17,20, 9,18,20, 10,19,20, * 11,16,21, 12,17,21, 13,18,21, 14,19,21, 15,20,21/ c i24(i,j) is an array identifying the i c binaries of the jth quaternary. data i24/ * 1, 2, 3, 4, 5, 6, 1, 2, 3, 7, 8, 9, 1, 4, 5, 7, 8,10, * 2, 4, 6, 7, 9,10, 3, 5, 6, 8, 9,10, 1, 2, 3,11,12,13, * 1, 4, 5,11,12,14, 2, 4, 6,11,13,14, 3, 5, 6,12,13,14, * 1, 7, 8,11,12,15, 2, 7, 9,11,13,15, 3, 8, 9,12,13,15, * 4, 7,10,11,14,15, 5, 8,10,12,14,15, 6, 9,10,13,14,15, * 1, 2, 3,16,17,18, 1, 4, 5,16,17,19, 2, 4, 6,16,18,19, * 3, 5, 6,17,18,19, 1, 7, 8,16,17,20, 2, 7, 9,16,18,20, * 3, 8, 9,17,18,20, 4, 7,10,16,19,20, 5, 8,10,17,19,20, * 6, 9,10,18,19,20, 1,11,12,16,17,21, 2,11,13,16,18,21, * 3,12,13,17,18,21, 4,11,14,16,19,21, 5,12,14,17,19,21, * 6,13,14,18,19,21, 7,11,15,16,20,21, 8,12,15,17,20,21, * 9,13,15,18,20,21, 10,14,15,19,20,21/ c i25(i,j) identifies the i binaries of c the jth quinary. data i25/ * 1, 2, 3, 4, 5, 6, 7, 8, 9,10, 1, 2, 3, 4, 5, 6,11,12,13,14, * 1, 2, 3, 7, 8, 9,11,12,13,15, 1, 4, 5, 7, 8,10,11,12,14,15, * 2, 4, 6, 7, 9,10,11,13,14,15, 3, 5, 6, 8, 9,10,12,13,14,15, * 1, 2, 3, 4, 5, 6,16,17,18,19, 1, 2, 3, 7, 8, 9,16,17,18,20, * 1, 4, 5, 7, 8,10,16,17,19,20, 2, 4, 6, 7, 9,10,16,18,19,20, * 3, 5, 6, 8, 9,10,17,18,19,20, 1, 2, 3,11,12,13,16,17,18,21, * 1, 4, 5,11,12,14,16,17,19,21, 2, 4, 6,11,13,14,16,18,19,21, * 3, 5, 6,12,13,14,17,18,19,21, 1, 7, 8,11,12,15,16,17,20,21, * 2, 7, 9,11,13,15,16,18,20,21, 3, 8, 9,12,13,15,17,18,20,21, * 4, 7,10,11,14,15,16,19,20,21, 5, 8,10,12,14,15,17,19,20,21, * 6, 9,10,13,14,15,18,19,20,21/ c c data i26/ * 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15, * 1, 2, 3, 4, 5, 6, 7, 8, 9,10,16,17,18,19,20, * 1, 2, 3, 4, 5, 6,11,12,13,14,16,17,18,19,21, * 1, 2, 3, 7, 8, 9,11,12,13,15,16,17,18,20,21, * 1, 4, 5, 7, 8,10,11,12,14,15,16,17,19,20,21, * 2, 4, 6, 7, 9,10,11,13,14,15,16,18,19,20,21, * 3, 5, 6, 8, 9,10,12,13,14,15,17,18,19,20,21/ c i34(i,j) is an array identifying the c ith ternary in the jth quaternary. data i34/ * 1, 2, 3, 4, 1, 5, 6, 7, 2, 5, 8, 9, 3, 6, 8,10, 4, 7, 9,10, * 1,11,12,13, 2,11,14,15, 3,12,14,16, 4,13,15,16, 5,11,17,18, * 6,12,17,19, 7,13,18,19, 8,14,17,20, 9,15,18,20, 10,16,19,20, * 1,21,22,23, 2,21,24,25, 3,22,24,26, 4,23,25,26, 5,21,27,28, * 6,22,27,29, 7,23,28,29, 8,24,27,30, 9,25,28,30, 10,26,29,30, * 11,21,31,32, 12,22,31,33, 13,23,32,33, 14,24,31,34, 15,25,32,34, * 16,26,33,34, 17,27,31,35, 18,28,32,35, 19,29,33,35, 20,30,34,35/ c i35(i,j) identifies the i ternaries of c the jth quinary. data i35/ * 1, 2, 3, 4, 5, 6, 7, 8, 9,10, 1, 2, 3, 4,11,12,13,14,15,16, * 1, 5, 6, 7,11,12,13,17,18,19, 2, 5, 8, 9,11,14,15,17,18,20, * 3, 6, 8,10,12,14,16,17,19,20, 4, 7, 9,10,13,15,16,18,19,20, * 1, 2, 3, 4,21,22,23,24,25,26, 1, 5, 6, 7,21,22,23,27,28,29, * 2, 5, 8, 9,21,24,25,27,28,30, 3, 6, 8,10,22,24,26,27,29,30, * 4, 7, 9,10,23,25,26,28,29,30, 1,11,12,13,21,22,23,31,32,33, * 2,11,14,15,21,24,25,31,32,34, 3,12,14,16,22,24,26,31,33,34, * 4,13,15,16,23,25,26,32,33,34, 5,11,17,18,21,27,28,31,32,35, * 6,12,17,19,22,27,29,31,33,35, 7,13,18,19,23,28,29,32,33,35, * 8,14,17,20,24,27,30,31,34,35, 9,15,18,20,25,28,30,32,34,35, * 10,16,19,20,26,29,30,33,34,35/ data i36/ * 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15,16,17,18,19,20, * 1, 2, 3, 4, 5, 6, 7, 8, 9,10,21,22,23,24,25,26,27,28,29,30, * 1, 2, 3, 4,11,12,13,14,15,16,21,22,23,24,25,26,31,32,33,34, * 1, 5, 6, 7,11,12,13,17,18,19,21,22,23,27,28,29,31,32,33,35, * 2, 5, 8, 9,11,14,15,17,18,20,21,24,25,27,28,30,31,32,34,35, * 3, 6, 8,10,12,14,16,17,19,20,22,24,26,27,29,30,31,33,34,35, * 4, 7, 9,10,13,15,16,18,19,20,23,25,26,28,29,30,32,33,34,35/ c i45(i,j) is an array identifying the c ith quaternary in the jth quinary. data i45/ * 1, 2, 3, 4, 5, 1, 6, 7, 8, 9, 2, 6,10,11,12, 3, 7,10,13,14, * 4, 8,11,13,15, 5, 9,12,14,15, 1,16,17,18,19, 2,16,20,21,22, * 3,17,20,23,24, 4,18,21,23,25, 5,19,22,24,25, 6,16,26,27,28, * 7,17,26,29,30, 8,18,27,29,31, 9,19,28,30,31, 10,20,26,32,33, * 11,21,27,32,34, 12,22,28,33,34, 13,23,29,32,35, 14,24,30,33,35, * 15,25,31,34,35/ c c data i46/ * 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15, * 1, 2, 3, 4, 5,16,17,18,19,20,21,22,23,24,25, * 1, 6, 7, 8, 9,16,17,18,19,26,27,28,29,30,31, * 2, 6,10,11,12,16,20,21,22,26,27,28,32,33,34, * 3, 7,10,13,14,17,20,23,24,26,29,30,32,33,35, * 4, 8,11,13,15,18,21,23,25,27,29,31,32,34,35, * 5, 9,12,14,15,19,22,24,25,28,30,31,33,34,35/ c i56(i,j) is an array identifying the c ith quinary in the jth hexary. data i56/ * 1, 2, 3, 4, 5, 6, 1, 7, 8, 9,10,11, 2, 7,12,13,14,15, * 3, 8,12,16,17,18, 4, 9,13,16,19,20, 5,10,14,17,19,21, * 6,11,15,18,20,21/ c------------------------------------------------------------------- c SUBSYSTEM ABSENT ARRAYS c data i16a/ * 6, 5, 4, 3, 2, 1, 7, 5, 4, 3, 2, 1, 7, 6, 4, 3, 2, 1, * 7, 6, 5, 3, 2, 1, 7, 6, 5, 4, 2, 1, 7, 6, 5, 4, 3, 1, * 7, 6, 5, 4, 3, 2/ c i24a(i,j,k) indexes the three binaries c in quaternary k that are absent from c the jth ternary of quaternary k. c Note that the elements of i24a are loaded c into a one dimensional array here, but c are used as a 3d array (3,4,35) data (i24a(zi),zi = 1,180)/ * 4, 5, 6, 2, 3, 6, 1, 3, 5, 1, 2, 4, 7, 8, 9, 2, 3, 9, * 1, 3, 8, 1, 2, 7, 7, 8,10, 4, 5,10, 1, 5, 8, 1, 4, 7, * 7, 9,10, 4, 6,10, 2, 6, 9, 2, 4, 7, 8, 9,10, 5, 6,10, * 3, 6, 9, 3, 5, 8, 11,12,13, 2, 3,13, 1, 3,12, 1, 2,11, * 11,12,14, 4, 5,14, 1, 5,12, 1, 4,11, 11,13,14, 4, 6,14, * 2, 6,13, 2, 4,11, 12,13,14, 5, 6,14, 3, 6,13, 3, 5,12, * 11,12,15, 7, 8,15, 1, 8,12, 1, 7,11, 11,13,15, 7, 9,15, * 2, 9,13, 2, 7,11, 12,13,15, 8, 9,15, 3, 9,13, 3, 8,12, * 11,14,15, 7,10,15, 4,10,14, 4, 7,11, 12,14,15, 8,10,15, * 5,10,14, 5, 8,12, 13,14,15, 9,10,15, 6,10,14, 6, 9,13/ data (i24a(zi),zi = 181,360)/ * 16,17,18, 2, 3,18, 1, 3,17, 1, 2,16, 16,17,19, 4, 5,19, * 1, 5,17, 1, 4,16, 16,18,19, 4, 6,19, 2, 6,18, 2, 4,16, * 17,18,19, 5, 6,19, 3, 6,18, 3, 5,17, 16,17,20, 7, 8,20, * 1, 8,17, 1, 7,16, 16,18,20, 7, 9,20, 2, 9,18, 2, 7,16, * 17,18,20, 8, 9,20, 3, 9,18, 3, 8,17, 16,19,20, 7,10,20, * 4,10,19, 4, 7,16, 17,19,20, 8,10,20, 5,10,19, 5, 8,17, * 18,19,20, 9,10,20, 6,10,19, 6, 9,18, 16,17,21, 11,12,21, * 1,12,17, 1,11,16, 16,18,21, 11,13,21, 2,13,18, 2,11,16, * 17,18,21, 12,13,21, 3,13,18, 3,12,17, 16,19,21, 11,14,21, * 4,14,19, 4,11,16, 17,19,21, 12,14,21, 5,14,19, 5,12,17/ data (i24a(zi),zi = 361,420)/ * 18,19,21, 13,14,21, 6,14,19, 6,13,18, 16,20,21, 11,15,21, * 7,15,20, 7,11,16, 17,20,21, 12,15,21, 8,15,20, 8,12,17, * 18,20,21, 13,15,21, 9,15,20, 9,13,18, 19,20,21, 14,15,21, * 10,15,20, 10,14,19/ c i25a(i,j,k) indexes the 4 binaries c absent from the jth quaternary of the c kth quinary, loaded as 1d array: data (i25a(zi),zi = 1,160)/ * 7, 8, 9,10, 4, 5, 6,10, 2, 3, 6, 9, 1, 3, 5, 8, * 1, 2, 4, 7, 11,12,13,14, 4, 5, 6,14, 2, 3, 6,13, * 1, 3, 5,12, 1, 2, 4,11, 11,12,13,15, 7, 8, 9,15, * 2, 3, 9,13, 1, 3, 8,12, 1, 2, 7,11, 11,12,14,15, * 7, 8,10,15, 4, 5,10,14, 1, 5, 8,12, 1, 4, 7,11, * 11,13,14,15, 7, 9,10,15, 4, 6,10,14, 2, 6, 9,13, * 2, 4, 7,11, 12,13,14,15, 8, 9,10,15, 5, 6,10,14, * 3, 6, 9,13, 3, 5, 8,12, 16,17,18,19, 4, 5, 6,19, * 2, 3, 6,18, 1, 3, 5,17, 1, 2, 4,16, 16,17,18,20, * 7, 8, 9,20, 2, 3, 9,18, 1, 3, 8,17, 1, 2, 7,16/ data (i25a(zi),zi = 161,320)/ * 16,17,19,20, 7, 8,10,20, 4, 5,10,19, 1, 5, 8,17, * 1, 4, 7,16, 16,18,19,20, 7, 9,10,20, 4, 6,10,19, * 2, 6, 9,18, 2, 4, 7,16, 17,18,19,20, 8, 9,10,20, * 5, 6,10,19, 3, 6, 9,18, 3, 5, 8,17, 16,17,18,21, * 11,12,13,21, 2, 3,13,18, 1, 3,12,17, 1, 2,11,16, * 16,17,19,21, 11,12,14,21, 4, 5,14,19, 1, 5,12,17, * 1, 4,11,16, 16,18,19,21, 11,13,14,21, 4, 6,14,19, * 2, 6,13,18, 2, 4,11,16, 17,18,19,21, 12,13,14,21, * 5, 6,14,19, 3, 6,13,18, 3, 5,12,17, 16,17,20,21, * 11,12,15,21, 7, 8,15,20, 1, 8,12,17, 1, 7,11,16/ data (i25a(zi),zi = 321,420)/ * 16,18,20,21, 11,13,15,21, 7, 9,15,20, 2, 9,13,18, * 2, 7,11,16, 17,18,20,21, 12,13,15,21, 8, 9,15,20, * 3, 9,13,18, 3, 8,12,17, 16,19,20,21, 11,14,15,21, * 7,10,15,20, 4,10,14,19, 4, 7,11,16, 17,19,20,21, * 12,14,15,21, 8,10,15,20, 5,10,14,19, 5, 8,12,17, * 18,19,20,21, 13,14,15,21, 9,10,15,20, 6,10,14,19, * 6, 9,13,18/ c c data i26a/ * 11,12,13,14,15, 7, 8, 9,10,15, 4, 5, 6,10,14, 2, 3, 6, 9,13, * 1, 3, 5, 8,12, 1, 2, 4, 7,11, 16,17,18,19,20, 7, 8, 9,10,20, * 4, 5, 6,10,19, 2, 3, 6, 9,18, 1, 3, 5, 8,17, 1, 2, 4, 7,16, * 16,17,18,19,21, 11,12,13,14,21, 4, 5, 6,14,19, 2, 3, 6,13,18, * 1, 3, 5,12,17, 1, 2, 4,11,16, 16,17,18,20,21, 11,12,13,15,21, * 7, 8, 9,15,20, 2, 3, 9,13,18, 1, 3, 8,12,17, 1, 2, 7,11,16, * 16,17,19,20,21, 11,12,14,15,21, 7, 8,10,15,20, 4, 5,10,14,19, * 1, 5, 8,12,17, 1, 4, 7,11,16, 16,18,19,20,21, 11,13,14,15,21, * 7, 9,10,15,20, 4, 6,10,14,19, 2, 6, 9,13,18, 2, 4, 7,11,16, * 17,18,19,20,21, 12,13,14,15,21, 8, 9,10,15,20, 5, 6,10,14,19, * 3, 6, 9,13,18, 3, 5, 8,12,17/ c c data i27a/ * 16,17,18,19,20,21, 11,12,13,14,15,21, 7, 8, 9,10,15,20, * 4, 5, 6,10,14,19, 2, 3, 6, 9,13,18, 1, 3, 5, 8,12,17, * 1, 2, 4, 7,11,16/ c i35a(i,j,k) indexes the 6 ternaries c absent from the jth quaternary of the c kth quinary, loaded as 1d array: data (i35a(zi),zi = 1,180)/ * 5, 6, 7, 8, 9,10, 2, 3, 4, 8, 9,10, 1, 3, 4, 6, 7,10, * 1, 2, 4, 5, 7, 9, 1, 2, 3, 5, 6, 8, 11,12,13,14,15,16, * 2, 3, 4,14,15,16, 1, 3, 4,12,13,16, 1, 2, 4,11,13,15, * 1, 2, 3,11,12,14, 11,12,13,17,18,19, 5, 6, 7,17,18,19, * 1, 6, 7,12,13,19, 1, 5, 7,11,13,18, 1, 5, 6,11,12,17, * 11,14,15,17,18,20, 5, 8, 9,17,18,20, 2, 8, 9,14,15,20, * 2, 5, 9,11,15,18, 2, 5, 8,11,14,17, 12,14,16,17,19,20, * 6, 8,10,17,19,20, 3, 8,10,14,16,20, 3, 6,10,12,16,19, * 3, 6, 8,12,14,17, 13,15,16,18,19,20, 7, 9,10,18,19,20, * 4, 9,10,15,16,20, 4, 7,10,13,16,19, 4, 7, 9,13,15,18/ data (i35a(zi),zi = 181,360)/ * 21,22,23,24,25,26, 2, 3, 4,24,25,26, 1, 3, 4,22,23,26, * 1, 2, 4,21,23,25, 1, 2, 3,21,22,24, 21,22,23,27,28,29, * 5, 6, 7,27,28,29, 1, 6, 7,22,23,29, 1, 5, 7,21,23,28, * 1, 5, 6,21,22,27, 21,24,25,27,28,30, 5, 8, 9,27,28,30, * 2, 8, 9,24,25,30, 2, 5, 9,21,25,28, 2, 5, 8,21,24,27, * 22,24,26,27,29,30, 6, 8,10,27,29,30, 3, 8,10,24,26,30, * 3, 6,10,22,26,29, 3, 6, 8,22,24,27, 23,25,26,28,29,30, * 7, 9,10,28,29,30, 4, 9,10,25,26,30, 4, 7,10,23,26,29, * 4, 7, 9,23,25,28, 21,22,23,31,32,33, 11,12,13,31,32,33, * 1,12,13,22,23,33, 1,11,13,21,23,32, 1,11,12,21,22,31/ data (i35a(zi),zi = 361,540)/ * 21,24,25,31,32,34, 11,14,15,31,32,34, 2,14,15,24,25,34, * 2,11,15,21,25,32, 2,11,14,21,24,31, 22,24,26,31,33,34, * 12,14,16,31,33,34, 3,14,16,24,26,34, 3,12,16,22,26,33, * 3,12,14,22,24,31, 23,25,26,32,33,34, 13,15,16,32,33,34, * 4,15,16,25,26,34, 4,13,16,23,26,33, 4,13,15,23,25,32, * 21,27,28,31,32,35, 11,17,18,31,32,35, 5,17,18,27,28,35, * 5,11,18,21,28,32, 5,11,17,21,27,31, 22,27,29,31,33,35, * 12,17,19,31,33,35, 6,17,19,27,29,35, 6,12,19,22,29,33, * 6,12,17,22,27,31, 23,28,29,32,33,35, 13,18,19,32,33,35, * 7,18,19,28,29,35, 7,13,19,23,29,33, 7,13,18,23,28,32/ data (i35a(zi),zi = 541,630)/ * 24,27,30,31,34,35, 14,17,20,31,34,35, 8,17,20,27,30,35, * 8,14,20,24,30,34, 8,14,17,24,27,31, 25,28,30,32,34,35, * 15,18,20,32,34,35, 9,18,20,28,30,35, 9,15,20,25,30,34, * 9,15,18,25,28,32, 26,29,30,33,34,35, 16,19,20,33,34,35, * 10,19,20,29,30,35, 10,16,20,26,30,34, 10,16,19,26,29,33/ c c data (i36a(zi),zi = 1,200)/ * 11,12,13,14,15,16,17,18,19,20, 5, 6, 7, 8, 9,10,17,18,19,20, * 2, 3, 4, 8, 9,10,14,15,16,20, 1, 3, 4, 6, 7,10,12,13,16,19, * 1, 2, 4, 5, 7, 9,11,13,15,18, 1, 2, 3, 5, 6, 8,11,12,14,17, * 21,22,23,24,25,26,27,28,29,30, 5, 6, 7, 8, 9,10,27,28,29,30, * 2, 3, 4, 8, 9,10,24,25,26,30, 1, 3, 4, 6, 7,10,22,23,26,29, * 1, 2, 4, 5, 7, 9,21,23,25,28, 1, 2, 3, 5, 6, 8,21,22,24,27, * 21,22,23,24,25,26,31,32,33,34, 11,12,13,14,15,16,31,32,33,34, * 2, 3, 4,14,15,16,24,25,26,34, 1, 3, 4,12,13,16,22,23,26,33, * 1, 2, 4,11,13,15,21,23,25,32, 1, 2, 3,11,12,14,21,22,24,31, * 21,22,23,27,28,29,31,32,33,35, 11,12,13,17,18,19,31,32,33,35/ data (i36a(zi),zi = 201,420)/ * 5, 6, 7,17,18,19,27,28,29,35, 1, 6, 7,12,13,19,22,23,29,33, * 1, 5, 7,11,13,18,21,23,28,32, 1, 5, 6,11,12,17,21,22,27,31, * 21,24,25,27,28,30,31,32,34,35, 11,14,15,17,18,20,31,32,34,35, * 5, 8, 9,17,18,20,27,28,30,35, 2, 8, 9,14,15,20,24,25,30,34, * 2, 5, 9,11,15,18,21,25,28,32, 2, 5, 8,11,14,17,21,24,27,31, * 22,24,26,27,29,30,31,33,34,35, 12,14,16,17,19,20,31,33,34,35, * 6, 8,10,17,19,20,27,29,30,35, 3, 8,10,14,16,20,24,26,30,34, * 3, 6,10,12,16,19,22,26,29,33, 3, 6, 8,12,14,17,22,24,27,31, * 23,25,26,28,29,30,32,33,34,35, 13,15,16,18,19,20,32,33,34,35, * 7, 9,10,18,19,20,28,29,30,35, 4, 9,10,15,16,20,25,26,30,34, * 4, 7,10,13,16,19,23,26,29,33, 4, 7, 9,13,15,18,23,25,28,32/ c c data i37a/ * 21,22,23,24,25,26,27,28,29,30,31,32,33,34,35, * 11,12,13,14,15,16,17,18,19,20,31,32,33,34,35, * 5, 6, 7, 8, 9,10,17,18,19,20,27,28,29,30,35, * 2, 3, 4, 8, 9,10,14,15,16,20,24,25,26,30,34, * 1, 3, 4, 6, 7,10,12,13,16,19,22,23,26,29,33, * 1, 2, 4, 5, 7, 9,11,13,15,18,21,23,25,28,32, * 1, 2, 3, 5, 6, 8,11,12,14,17,21,22,24,27,31/ c c data (i46a(zi),zi = 1,220)/ * 6, 7, 8, 9,10,11,12,13,14,15, 2, 3, 4, 5,10,11,12,13,14,15, * 1, 3, 4, 5, 7, 8, 9,13,14,15, 1, 2, 4, 5, 6, 8, 9,11,12,15, * 1, 2, 3, 5, 6, 7, 9,10,12,14, 1, 2, 3, 4, 6, 7, 8,10,11,13, * 16,17,18,19,20,21,22,23,24,25, 2, 3, 4, 5,20,21,22,23,24,25, * 1, 3, 4, 5,17,18,19,23,24,25, 1, 2, 4, 5,16,18,19,21,22,25, * 1, 2, 3, 5,16,17,19,20,22,24, 1, 2, 3, 4,16,17,18,20,21,23, * 16,17,18,19,26,27,28,29,30,31, 6, 7, 8, 9,26,27,28,29,30,31, * 1, 7, 8, 9,17,18,19,29,30,31, 1, 6, 8, 9,16,18,19,27,28,31, * 1, 6, 7, 9,16,17,19,26,28,30, 1, 6, 7, 8,16,17,18,26,27,29, * 16,20,21,22,26,27,28,32,33,34, 6,10,11,12,26,27,28,32,33,34, * 2,10,11,12,20,21,22,32,33,34, 2, 6,11,12,16,21,22,27,28,34/ data (i46a(zi),zi = 221,420)/ * 2, 6,10,12,16,20,22,26,28,33, 2, 6,10,11,16,20,21,26,27,32, * 17,20,23,24,26,29,30,32,33,35, 7,10,13,14,26,29,30,32,33,35, * 3,10,13,14,20,23,24,32,33,35, 3, 7,13,14,17,23,24,29,30,35, * 3, 7,10,14,17,20,24,26,30,33, 3, 7,10,13,17,20,23,26,29,32, * 18,21,23,25,27,29,31,32,34,35, 8,11,13,15,27,29,31,32,34,35, * 4,11,13,15,21,23,25,32,34,35, 4, 8,13,15,18,23,25,29,31,35, * 4, 8,11,15,18,21,25,27,31,34, 4, 8,11,13,18,21,23,27,29,32, * 19,22,24,25,28,30,31,33,34,35, 9,12,14,15,28,30,31,33,34,35, * 5,12,14,15,22,24,25,33,34,35, 5, 9,14,15,19,24,25,30,31,35, * 5, 9,12,15,19,22,25,28,31,34, 5, 9,12,14,19,22,24,28,30,33/ c data i47a/ * 16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35, * 6, 7, 8, 9,10,11,12,13,14,15,26,27,28,29,30,31,32,33,34,35, * 2, 3, 4, 5,10,11,12,13,14,15,20,21,22,23,24,25,32,33,34,35, * 1, 3, 4, 5, 7, 8, 9,13,14,15,17,18,19,23,24,25,29,30,31,35, * 1, 2, 4, 5, 6, 8, 9,11,12,15,16,18,19,21,22,25,27,28,31,34, * 1, 2, 3, 5, 6, 7, 9,10,12,14,16,17,19,20,22,24,26,28,30,33, * 1, 2, 3, 4, 6, 7, 8,10,11,13,16,17,18,20,21,23,26,27,29,32/ c c data i57a/ * 7, 8, 9,10,11,12,13,14,15,16,17,18,19,20,21, * 2, 3, 4, 5, 6,12,13,14,15,16,17,18,19,20,21, * 1, 3, 4, 5, 6, 8, 9,10,11,16,17,18,19,20,21, * 1, 2, 4, 5, 6, 7, 9,10,11,13,14,15,19,20,21, * 1, 2, 3, 5, 6, 7, 8,10,11,12,14,15,17,18,21, * 1, 2, 3, 4, 6, 7, 8, 9,11,12,13,15,16,18,20, * 1, 2, 3, 4, 5, 7, 8, 9,10,12,13,14,16,17,19/ end subroutine assip (ip) c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k3=1200,k8=9,l2=5) character*8 names,ipnms common/ cst23 /a(10,10),b(11),ipvt(10),idv(10), * iophi,idphi,iiphi,iflg1 * / cst5 /v(l2),tr,pr,r,ps * / cst7 /iflag/ cst29 /vip(l2,k3),ipid(k3,k8),ipct * / cst9 /vmax(l2),vmin(l2),ddv(l2) * / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5 * / cst6 /icomp,istct,iphct,icp/ cst8 /names(k1) * / cst34 /ip1(k3),iipct,loop/ csta3 /ipnms(k3,k8) * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 c----------------------------------------------------------------------- icp1=icp+1 icp2 =icp+2 idv(icp1) =iophi idv(icp2) =idphi c determine if a new ip was found if (ipct.eq.0) goto 40 jps = 1 5 ips = jps do 10 i = ips, ipct do 20 j = 1,icp2 do 30 k= 1,icp2 30 if (idv(k).eq.ipid(i,j)) goto 20 c check if invariant point is equivalent c to the ith ip already found goto 10 20 continue c the new ip assemblage is equivalent to one found c earlier, check if the conditions match: do 25 j = 1, 5 if (dabs(vip(j,i)-v(j)).gt.ddv(j).and. * ddv(j).ne.0.) then if (i.lt.ipct) then jps = i + 1 goto 5 end if goto 40 end if 25 continue ip=i iflag= 0 goto 99 10 continue c assign the new elements of ipid 40 ipct=ipct + 1 if (ipct.gt.k3) call error (182,r,k3,'ASSIP') c check if the ip was encountered in c earlier runs or tiers if (io5.eq.1) goto (80),loop if (iipct.eq.0) goto 80 do 50 i = 1,iipct do 60 j = 1,icp2 do 70 k= 1,icp2 70 if (names(idv(k)).eq.ipnms(i,j)) goto 60 goto 50 60 continue c a match was found with ip i ip1(ipct) =i goto 100 50 continue c no match found 80 iipct=iipct + 1 if (iipct.gt.k3) call error (182,r,k3,'ASSIP') do 110 i = 1, icp2 110 ipnms(iipct,i) =names(idv(i)) ip1(ipct) =iipct 100 do 90 j = 1, icp2 90 ipid(ipct,j) =idv(j) vip(1,ipct) = v(1) vip(2,ipct) = v(2) vip(3,ipct) = v(3) vip(4,ipct) = v(4) vip(5,ipct) = v(5) ip=ipct 99 end c subroutine assir c---------------------------------------------------- c iirct - cumulative counter of univariant rxns. c irct - execution counter " " ". c ir1(irct) - index (=iirct) of the irct(th) rxn identified in an c execution. c irv(iirct) - number of phases in the iirct(th)(=ir1(irct)) rxn c (set in balanc as ivct cst25). c ivarrx(iirct) - variance flag of the iirct(th) rxn (set in balanc as c ivar in cst61. c jrchk(iirct) - a flag indicating if the iirct(th) rxn has been c encountered during an execution (initialized in input1 c and set in assir) c irchk(iirct) - a flag indicating if the iirct(th) rxn occurs on the c edge of a diagram (initialized in newhld and set in c coface or sfol1). c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k2=1200,k3=1200,k8=9,l2=5) c character*8 ipnms,irnms,names c common/ cst5 /v(l2),tr,pr,r,ps/ cst80 /jrchk(k2),irchk(k2) * / cst31 /vn(k2,11),ir(k2,11),irct,ird/ cst8 /names(k1) * / cst25 /vnu(11),idr(11),ivct * / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5 * / cst34 /ip1(k3),iipct,loop/ csta3 /ipnms(k3,k8) * / cst35 /ir1(k2),iirct,irpct,irv(k2)/ csta1 /irnms(k2,11) * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar c----------------------------------------------------------------------- c check if the equilibrium is stored in c the permanent rxn list 'irnms',if io5 = 0 if (iirct.eq.0) goto 40 do 70 i = 1, iirct if (ivct.ne.irv(i)) goto 70 do 80 j = 1, ivct do 90 k= 1, ivct c it is necessary to use names here c for testing against old reaction lists c cause the pointers are problem dependent 90 if (irnms(i,k).eq.names(idr(j))) goto 80 c no match with irnms(i,k) try i + 1. goto 70 80 continue c match found with irnms(i,k) ird=i if (jrchk(i).eq.1) then c whats this??? c if (irchk(i).eq.0) then c irct = irct + 1 c if (irct.gt.k2) call error (181,r,k2,'ASSIR') c ir1(irct) = i c end if goto 99 end if goto 60 70 continue c the reaction is new: 40 iirct = iirct + 1 if (iirct.gt.k2) call error (181,r,k2,'ASSIP') ird = iirct c the rxn was in the permanent list but c is new for the current execution. 60 if (jrchk(ird).eq.0) irct = irct + 1 if (irct.gt.k2) call error (181,r,k2,'ASSIR') ir1(irct) = ird irv(ird) = ivct ivarrx(ird) = ivar do 110 i = 1, ivct vn(ird,i) = vnu(i) ir(ird,i) = idr(i) 110 irnms(ird,i) = names(idr(i)) c match 99 end subroutine assri (ier) c---------------------------------------------------- c iirct - cumulative counter of univariant rxns. c irct - execution counter " " ". c ir1(irct) - index (=iirct) of the irct(th) rxn identified in an c execution. c irv(iirct) - number of phases in the iirct(th)(=ir1(irct)) rxn c (set in balanc as ivct cst25). c ivarrx(iirct) - variance flag of the iirct(th) rxn (set in balanc as c ivar in cst61. c jrchk(iirct) - a flag indicating if the iirct(th) rxn has been c encountered during an execution (initialized in input1 c and set in assir) c irchk(iirct) - a flag indicating if the iirct(th) rxn occurs on the c edge of a diagram (initialized in newhld and set in c coface or sfol1). c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k2=1200,k3=1200,k6=7,k8=9,l2=5) character*8 ipnms,irnms,names common/ cst5 /v(l2),tr,pr,r,ps/ cst80 /jrchk(k2),irchk(k2) * / cst31 /vn(k2,11),ir(k2,11),irct,ird/ cst8 /names(k1) * / cst25 /vnu(11),idr(11),ivct * / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5 * / cst34 /ip1(k3),iipct,loop/ csta3 /ipnms(k3,k8) * / cst35 /ir1(k2),iirct,irpct,irv(k2)/ csta1 /irnms(k2,11) * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar * / cst29 /vip(l2,k3),ipid(k3,k6+2),ipct * / cst9 /vmax(l2),vmin(l2),dv(l2) c----------------------------------------------------------------------- c check if the equilibrium is stored in c the permanent rxn list 'irnms',if io5 = 0 if (iirct.eq.0) goto 40 do 70 i = 1, iirct if (ivct.ne.irv(i)) goto 70 do 80 j = 1, ivct do 90 k= 1, ivct c it is necessary to use names here c for testing against old reaction lists c cause the pointers are problem dependent 90 if (irnms(i,k).eq.names(idr(j))) goto 80 c no match with irnms(i,k) try i + 1. goto 70 80 continue c match found with irnms(i,k), now check c that the conditions are different, if c ier = 1, univeq failed, skip the check goto (75), ier if (dabs(vip(iv1,i)-v(iv1)).gt.0.05*dv(iv1)) goto 40 75 ird = i ier = 1 if (jrchk(i).eq.1) goto 99 goto 60 70 continue c the reaction is new: 40 iirct = iirct + 1 ier = 0 if (iirct.gt.k2) call error (181,r,k2,'ASSRI') ird = iirct c the rxn was in the permanent list but c is new for the current execution. 60 if (jrchk(ird).eq.0) irct = irct + 1 if (irct.gt.k2) call error (181,r,k2,'ASSIP') ir1(irct) = ird irv(ird) = ivct ivarrx(ird) = ivar do 110 i = 1, ivct vn(ird,i) = vnu(i) ir(ird,i) = idr(i) 110 irnms(ird,i) = names(idr(i)) c match 99 end subroutine balanc (b,idv,idphi,ier) c----------------------------------------------------------------------- c balanc balances reactions between the phases in idv and idphi c if the system is saturated with any phases these are taken c into account. once the reaction is determined balanc determines c the thermodynamic parameters needed for subroutine univeq and c saves them in the array 'rxn'. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k2=1200,k3=1200,k4=16,k6=7, * h5=7,h6=500) dimension a(10,10), ipvt(10), idv(10), b(11) common/ cst1 /thermo(k4,k1),vf(2,k1),vs(h5,k1),uf(2),us(h5) * / cst6 /icomp,istct,iphct,icp/ cst25 /vnu(11),idr(11),ivct * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn/ cst20 /comps(k1,12) * / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2 * / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar * / cst88 /vu(2,k1),jprct,jmct,jmyn/ cst12 /cp(k6,k1) * / cst201 /vuf(2),vus(h5),iffr,isr c----------------------------------------------------------------------- c assemble transpose of the composition c matrix for an assemblage of icp phases ivar =1 do 10 i = 1, icp do 20 j = 1, icp 20 a(j,i) = cp(j,idv(i)) 10 idr(i) =idv(i) c assemble the composition vectore c for the icp+1 th phase tangent to c the free energy plane. ip=icp+1 do 70 i = 1, icp 70 b(i) =comps(idphi,i) idr(ip) =idphi b(ip) = -1.d0 c balance the equilibrium call factor (a,icp,ipvt,ier) goto (999),ier call subst (a,ipvt,icp,b) c eliminate phases with vnu= 0 ivct= 0 do 80 i = 1, ip if (dabs(b(i)).lt.1.d-05) goto 80 ivct=ivct + 1 vnu(ivct) =b(i) idr(ivct) =idr(i) 80 continue c test to determine if the reaction c is pseudo-univariant, and set flag c ivar =2 if so, and =1 otherwise. iv1=ivct-1 do 140 i = 1, iv1 if (ikp(idr(i)).eq.0) goto 140 i1=i + 1 do 150 j =i1,ivct if (ikp(idr(i)).ne.ikp(idr(j))) goto 150 ivar =2 c pseudo-univariant reject rxn if isudo=1 goto (999),isudo goto 170 c note that this test will fail to c distinguish between critical equilibria c and pseudo-univariant equilibria. 150 continue 140 continue 170 goto (40),isyn c determine stoichiometric coefficients c of saturated components: c (these aren't used for anything real) isr = 1 do 30 j = 1, isat vus(j) = 0.d0 do 120 i = 1, ivct 120 vus(j) = vus(j) + vnu(i) * comps(idr(i),icp+j) c flag isr is not used. if (vus(j).ne.0.d0) isr = 0 30 continue c ok, now vus has the total deltas c determine stoichiometric coefficients c of saturated phase components: iffr = 1 40 do 50 j = 1,2 vuf(j) = 0.d0 do 130 i = 1, ivct 130 vuf(j) = vuf(j) + vnu(i) * vf(j,idr(i)) c this is only the total of saturated c phase components for the reaction c (doesn't include saturated phases). c flag iffr is used 50 if (vuf(j).ne.0.d0) iffr = 0 c et finis 999 end subroutine chkass (jdv,ivd,ivi,iste) c----------------------------------------------------------------------- c chkass determines if a divariant assemblage (jdv) generated by c newass has already been identified. there are three possible cases: c in the reaction as determined by routine balanc. c (1) the assemblage is unique: the indices of the phases, the initial c conditions of stability, and the coordinate frame edge are recorde c (2) the assemblage has already been defined and it's stability field c determined: quit c (3) the assemblage has already been defined but it's stability field c has not been determined: reset the initial conditions of stability c to shorten traverses by the routine search, if possible. c input: jdv(icp) - a vector containing the indices of phases in a c divariant assemblage. c ipart - 1 if testing original divariant assemblages, 2 if c testing new divariant assemblages. c index - index of the divariant assemblage in array id7f c (ipart=1) or idns (ipart=2) currently being tested. c referenced by: newass c references to: no external references c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k2=1200,k3=1200,k6=7,l2=5) parameter (j9=1200,h7=6000,nb=20,k5=12) character*8 names integer jdv(k6), icrap common/ cst6 /icomp,istct,iphct,icp/ cst5 /v(l2),tr,pr,r,ps * / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar * / cst8 /names(k1)/ cst96 /id7f(k6,j9),i7fct * / cst65 /vt(j9),vvt(h7),jok(j9),jnk(h7),idns(k6,h7),idnct * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst99 /ipart,index,kok(j9),knk(h7)/ cst302 /jasmbl(j9) * / cst300 /cblk(nb,k5),atw(k1),vol(k1),iblk,jcby,jcth * / cst303 /vti(j9),vvti(h7)/ cst9 /vmax(l2),vmin(l2),dv(l2) save icrap data icrap /0/ c----------------------------------------------------------------------- c check for a match with an original c assemblage: do 10 i = 1, i7fct do 20 j = 1, icp do 30 k = 1, icp 30 if (jdv(k).eq.id7f(j,i)) goto 20 c no match with the ith assemblage goto 10 20 continue c jdv matches the ith assemblage c check if the bulk composition is c visible: if (jasmbl(i).eq.0.and.jcby.eq.1) goto 99 c check if within known stability c field of the assemblage if (iste.lt.kok(i)) then goto 99 else if (iste.eq.kok(i)) then if (vt(i).ge.v(ivd)) goto 99 c it extends known stability field: else if (ipart.eq.1.and.i.gt.index) then c if the assemblage hasn't been tested: call sreset (jok(i),iste,vt(i),v(ivd),vti(i),v(ivi)) goto 99 end if c the assemblage has already been tested c so count it as a new assemblage: goto 70 c 10 continue c no match with an original assemblage c try new assemblages: 70 if (idnct.eq.0) goto 100 do 40 i = 1, idnct do 50 j = 1, icp do 60 k = 1, icp 60 if (jdv(k).eq.idns(j,i)) goto 50 c no match with the ith assemblage goto 40 50 continue c jdv matches the ith assemblage c first check if within known stability c field of the assemblage if (iste.lt.knk(i)) then goto 99 else if (iste.eq.knk(i)) then if (iste.lt.3) then if (vvt(i).ge.v(ivd)-dv(ivd)) goto 99 else if (vvt(i).le.v(ivd)+dv(ivd)) goto 99 end if c it extends known stability field: c if the assemblage hasn't been tested c then reset start conditions else if (ipart.eq.1) then call sreset (jnk(i),iste,vvt(i),v(ivd),vvti(i),v(ivi)) goto 99 else if (i.gt.index) then call sreset (jnk(i),iste,vvt(i),v(ivd),vvti(i),v(ivi)) goto 99 end if c otherwise save as a new assemblage goto 100 40 continue c the assemblage is unique, store c parameters: 100 goto (99),icrap c do bulk composition check if (iblk.ne.0) then call bulkck (jdv,jbfg,ier) if (jcby.eq.1.and.jbfg.eq.0) goto 99 end if idnct = idnct + 1 if (idnct.gt.h7) then call warn (205,r,h7,'CHKASS') idnct = h7 icrap = 1 end if do 110 i = 1, icp 110 idns(i,idnct) = jdv(i) vvt(idnct) = v(ivd) vvti(idnct) = v(ivi) jnk(idnct) = iste knk(idnct) = iste 99 end subroutine coface (iovd,iovi,odiv,iste,irend) c---------------------------------------------------------------------- c once nohold has established that a stable equilibrium occurs on c the boundary of a diagram, coface traces the equilibrium condiitons c within the diagram. if the equilibrium terminates at an c invariant point coface traces the equilibria which emanate c from the invariant point and any invariant points encountered c subsequently. coface does not detect indifferent points. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k2=1200,k3=1200,k8=9,l2=5) character*8 names,irnms common/ cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar * / cst29 /vip(l2,k3),ipid(k3,k8),ipct * / cst7 /iflag/ cst80 /jrchk(k2),irchk(k2) * / cst23 /a(10,10),b(11),ipvt(10),idv(10), * iophi,idphi,iiphi,iflg1 * / cst31 /vn(k2,11),ir(k2,11),irct,ird * / cst25 /vnu(11),idr(11),ivct/ cst8 /names(k1) * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst35 /ir1(k2),iirct,irpct,irv(k2)/ csta1 /irnms(k2,11) * / cst5 /v(l2),tr,pr,r,ps/ cst32 /ptx(450),logx,ipt2 * / cst105 /ibug(k2)/ cst103 /isec,icopt,ifull,imsg,io3p c----------------------------------------------------------------------- ivi =iovi ivd=iovd jflg= 1 call balanc (b,idv,idphi,ier) c singular concentration matrix if (ier.eq.1) call error (68,v(1),ier,'COFACE') c call newass to define new (pseudo-) c divariant assemblages generated by c the reaction defined by balanc: call newass (b,idv,idphi,ivd,ivi,iste) c if the rxn is univariant (ivar =1) c continue. goto (80),ivar c the rxn is pseudo-univariant continue c only if isudo= 0 (input). goto (9999),isudo c call assir to determine if the c reaction has already been identified c and if so quit, there is a small c possibility that as a result of c quitting at this stage some equilibria c will not be identified, if this is a c significant concern the following c computed goto should be deleted. 80 goto (90), icopt c allow duplicated reactions in c mixed variable diagrams: call univeq (ivd,ier) c assign reaction: call assri (ier) goto (9999), ier call delrxn vip(1,ird) = v(1) vip(2,ird) = v(2) vip(3,ird) = v(3) vip(4,ird) = v(4) vip(5,ird) = v(5) call outrxn (0,0) goto 9999 90 call assir c this section replaced by next: c skip equilibria that have already c been defined on an edge, there is some c possibility that this may result in an c an incomplete diagram, but it saves c some time and redundancy in output. c check if the univariant rxn has alread c been found on the edge of a diagram. c set flags irchk and ier: c if (irchk(ird).eq.1) goto 9999 c call svrend to see if the end point c has already been found, if it has ier c is 1, and traverse can be skipped. call svrend (ird,irend,ier) goto (9999),ier irchk(ird) = 1 call univeq (ivd,ier) c get deltas for dependent extensities: call delrxn c if univeq fails on a bounding edge c write error message and return: goto (9000,9000),ier c save the index of the c+1th phase c as iophi. iophi =idphi c factor the concentration matrix: call pivots (ier) if (ier.eq.1) call error (69,v(1),ier,'COFACE') lphi = 0 c set the increment for the iv div = odiv c initialize counters ipt2 = 0 icter = 0 c assign the 1st point call assptx c follow the equilibrium 60 ikwk = 0 jflg = 0 call sfol1 (ivd,ivi,ier,div,lphi,ikwk,irend) c sfol1 returns ier = 1 goto (70,70),ier ivi = iovi ivd = iovd goto (10),iflag goto 9999 70 call switch (div,ivi,ivd,jer) goto (75),jer icter = icter + 1 if (icter.lt.4) goto 60 75 call warn (10,v(1),ier,'COFACE') call outrxn (ipct,2) ibug(iirct) = 1 ivi =iovi ivd=iovd c return on error goto 9999 c a new invariant point has been c encountered, generate the new c univariant fields. c iopct=the total number of ip's encountered c before each call to sfol2. c inpct=the total number after each call to c sfol2. 10 iopct=ipct call sfol2 (iovi,iovd,iopct,irend) if (iopct.eq.ipct) goto 9999 iopct=iopct + 1 inpct=ipct c loop for new invariant points 30 do 40 i =iopct,inpct 40 call sfol2 (iovi,iovd,i,irend) if (ipct.eq.inpct) goto 9999 iopct=inpct + 1 inpct=ipct goto 30 c error in univeq: 9000 call warn (79,v(1),ird,'COFACE') ipt2 = 0 call outrxn (ipct,2) ibug(iirct) = 1 9999 iflg1= 0 goto (999),jflg goto (999),io3p write (n3,1020) 1020 format (' Network traced, resuming boundary search.',/) 999 end subroutine delvar (dv,iflag,iflg1) c----------------------------------------------------------------------- c delvar determines if a variable increment 'dv' should c be incremented or decremented on a traverse to locate c to locate conditions where an equilibrium is metastable with respect c to exactly one phase (iflag = 1). the flag iflg1 indicates if the c equilibrium was stable (iflg1 = 0) or metastable (iflg1 = 1) at the c last conditions tested. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) c the equilibrium is stable (iflag = 0): goto (20,20),iflag c and was stable before (iflg1 = 0) then goto (10),iflg1 c continue looking in the same direction: return c and was metastable before (iflg1 = 1) then c switch sign and halve dv: 10 dv = -dv/2.d0 iflg1 = 0 return c the equilibrium is metastable (iflag = 2): 20 goto (30),iflg1 c and was stable before (iflg1 = 0) then c switch sign and halve dv: dv = -dv/2.d0 iflg1 = 1 c and was metastable before (iflg1 = 1) then c continue looking in the same direction: 30 end subroutine factor (a,n,ipvt,ier) c----------------------------------------------------------------------- c factor is a subroutine which calculates the triangular c decompositions of the matrix 'a'. factor is modified from c the subroutine of the same name given by conte and de boor c in 'elementary numerical analysis', mcgraw-hill, 1980. c factor uses scaled partial pivoting. c c input a- an n by n array containing the elements of matrix a. c n- the dimension of the matrix a. c output a- an n by n array containing the upper, u, and lower, l, c triangular decompositions of input matrix a. c ipvt- a vector indicating that row ipvt(k) was used to c eliminate the a(n,k). c ier- a flag, zero if a is of rank = n, and 1 if a is of c lower rank. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) c dimension a(10,10),d(10),ipvt(10) c----------------------------------------------------------------------- ier = 0 c initialize ipvt,d do 10 i = 1,n ipvt(i) = i rmax = 0.d0 do 20 j = 1,n 20 rmax = dmax1(rmax,dabs(a(i,j))) c ax = b is singular if rmax = 0 if (rmax.eq.0.d0) goto 9000 10 d(i) = rmax c begin decomposition: nm1 = n-1 do 30 i = 1,nm1 ip1 = i+1 c determine pivot row (istr). rmax = dabs(a(i,i))/d(i) istr = i do 40 j = ip1,n tmax = dabs(a(j,i))/d(j) if (tmax.le.rmax) goto 40 rmax = tmax istr = j 40 continue if (rmax.eq.0.d0) goto 9000 c if istr gt i, make i the pivot row c by interchanging it with row istr. if (istr.le.i) goto 50 j = ipvt(istr) ipvt(istr) = ipvt(i) ipvt(i) = j temp = d(istr) d(istr) = d(i) d(i) = temp do 60 j = 1,n temp = a(istr,j) a(istr,j) = a(i,j) 60 a(i,j) = temp c eliminate x(k) from rows k+1,...,n. 50 do 70 j = ip1,n a(j,i) = a(j,i)/a(i,i) ratio = a(j,i) do 80 k = ip1,n 80 a(j,k) = a(j,k)-ratio*a(i,k) 70 continue 30 continue if (a(n,n).eq.0.d0) ier = 1 return c algoritmic singularity. 9000 ier = 1 end subroutine gall (istart,iend) c----------------------------------------------------------------------- c subroutine gall computes molar free energies of the phases with c indices between and equal to istart and iend. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) parameter (k1=1200, k4=16, h5=7, h6=500) common/ cst1 /thermo(k4,k1),vf(2,k1),vs(h5,k1),uf(2),us(h5) * / cst2 /g(k1)/ cst6 /icomp,istct,iphct,icp * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn * / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2 c----------------------------------------------------------------------- c compute the chemical potential c of the projected components. call uproj do 10 id = istart, iend g(id) = gphase(id) - uf(1)*vf(1,id) - uf(2)*vf(2,id) goto (10),isyn do 20 j = 1, isat 20 g(id) = g(id) - us(j) * vs(j,id) 10 continue end function gproj (id) c----------------------------------------------------------------------- c gproj computes the projected molar free energy of the c phase with index id. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) parameter (k1=1200, k4=16, h5=7, h6=500) common/ cst1 /thermo(k4,k1),vf(2,k1),vs(h5,k1),uf(2),us(h5) * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn * / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2 c----------------------------------------------------------------------- gproj = gphase(id)-uf(1)*vf(1,id)-uf(2)*vf(2,id) goto (99),isyn do 10 j = 1, isat 10 gproj = gproj - vs(j,id) * us(j) 99 end subroutine flipit (ddv,vst,ivd,ist,inow,jer) c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,l2=5) character*8 names common/ cst23 /a(10,10),b(11),ipvt(10),idv(10), * iophi,idphi,iiphi,iflg1 * / cst5 /v(l2),x(4)/ cst7 /iflag/ cst6 /icomp,istct,iphct,icp * / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5/ cst38 /iedge * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst9 /vmax(l2),vmin(l2),dv(l2)/ cst8 /names(k1) c----------------------------------------------------------------------- c initialization call gall (istct,iphct) call asschk c the goto 20 is a dummy to avoid c compiler errors, change 90 to 20 to c get backup search, ha ha c goto (10,90,20),iflag if (iflag.ne.0.and.v(ivd).eq.vst.and.ist.eq.inow) then c write (*,1000) goto 90 end if goto (10,20,20),iflag c the assemblage is stable, return jer = 0 return c the assemblage is metastable with c respect to only one phase: 10 jer = 1 return c the assemblage is metastable with c respect to >1 phase, therefore c it must be stable at some condition c between the present condition and the c conditions of the univariant c equilibrium that generated the c assemblage. flipit does a reverse c traverse to find this condition. c initialize: 20 iflg1 = 1 ddv = -ddv c begin search: 100 v(ivd) = v(ivd)+ddv c out of range?: goto (120,120,110,110),ist c 110 if (v(ivd).gt.vmax(ivd)) v(ivd) = vmax(ivd) if (v(ivd).gt.vst) goto 130 ddv = dabs(ddv)/2.d0 v(ivd) = vst iflg1= 0 goto 100 c 120 if (v(ivd).lt.vmin(ivd)) v(ivd) = vmin(ivd) if (v(ivd).lt.vst) goto 130 ddv = -dabs(ddv)/2.d0 v(ivd) = vst iflg1 = 0 goto 100 c calculate phase energies: 130 call gall (istct,iphct) c test stability of the assemblage: call asschk c check if search is in range: goto (150,150,140,140),ist 140 if ((v(ivd).ge.vmax(ivd)).and.(iflag.gt.0)) goto 60 goto 160 150 if ((v(ivd).le.vmin(ivd)).and.(iflag.gt.0)) goto 60 c iflag=1, found a stable equilibrium c in bounds: if (iflag.eq.1) then write (n8,*) 'flipit worked please tell me, if youre not me' end if c 160 goto (10),iflag c refine increment if necessary: call delvar (ddv,iflag,iflg1) c check if increment has been refined c beyond the tolerance which is c arbitrarily set at 1.d-8. if (dabs(ddv).lt.1.d-8) goto 9000 c increment conditions again: goto 100 c next traverse. 60 call warn (76,vst,ist,'FLIPIT') goto 90 c write warning message: 9000 write (n8,1010) ivd,dv(ivd),(names(idv(j)),j = 1, icp) 90 jer = 2 c done: 1000 format (/,' **warning ver046** programming error, an assemblage', * ' is metastable at vst and ist in FLIPIT.',/) 1010 format (/,' **warning ver045** FLIPIT: > 1 equilibrium', * ' occurs within the',/,' minimum search increment for', * ' variable: ',i1,', this often occurs as YCO2 => 1', * ' or => 0, you may be able to correct this',/, * ' by reducing the default increment for this variable', * ' (',g12.3,') in the computational option file.',/, * ' Equilibria involving the following assemblage may', * ' not be delineated:',/,1x,7(1x,a8)) end subroutine outdel c----------------------------------------------------------------------- c output deltas of dependent extensities for a reaction c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k6=7,l2=5,h5=7,h6=500) character*5 cname,vname*8,gname*78,xname*18,exten(2)*7,names*8 common/ csta2 /xname(k6),vname(l2),gname(l2)/ csta4 /cname(12) * / cst21 /du(2),dv(2),jds(h5),ifr * / cst88 /vu(2,k1),jprct,jmct,jmyn * / cst6 /icomp,istct,iphct,icp * / cst201 /vuf(2),vus(h5),iffr,isr * / cst25 /vnu(11),idr(11),ivct/ cst8 /names(k1) * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn * / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2 * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 save exten c this is a bullshit trick and c will cause errors of someone c uses a function other than G. data exten/'-V(j/b)','S(j/k)'/ c composant stoichiometry: if (isyn.eq.0) then do 1 i = 1, isat 1 write (n3,1000) cname(icp+i),vus(i),names(jds(i)) end if c fluid stoichiometry: if (ifyn.eq.0) then do 2 i = 1, 2 2 if (iff(i).ne.0) write (n3,1010) names(i),vuf(i) end if c mobile components: if (jmyn.eq.0) then do 3 i = 1, jmct 3 write (n3,1020) cname(jprct+i),du(i),vname(3+i) end if c dependent extensities 1 and 2, c normally entropy and negative c volume: do 4 i = 1, 2 4 write (n3,1030) exten(i), dv(i),vname(i) 1000 format (10x,'Delta(',a7,') =',g9.3,1x, * '(saturated composant=',a8,')') 1010 format (10x,'Delta(',2x,a5,') =',g9.3,1x, * '(saturated phase component)') 1020 format (10x,'Delta(',a7,') =',g9.3,1x, * '(dependent conjugate of ',a8,')') 1030 format (10x,'Delta(',a7,') =',g9.3,1x, * '(dependent conjugate of ',a8,')') end subroutine delrxn c----------------------------------------------------------------------- c delrxn computes the change in the dependent extensities for a reaction c these are returned in cst21, delrxn is always called after UNIVEQ, so c the free energy at the present conditions is zero, if called c independently, then grxn must be called first. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k2=1200,k4=16,l2=5,h5=7,h6=500) common/ cst5 /v(5),tr,pr,r,ps/ cst21 /du(2),dv(2),jds(h5),ifr * / cst87 /delt(l2),dtol,utol,ptol/ cst102 /ivfl * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn * / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2 * / cst201 /vuf(2),vus(h5),iffr,isr/ cst208 /ifct,idfl * / cst25 /vnu(11),idr(11),ivct/ cst24 /ipot,jv(l2),iv(l2) * / cst31 /vn(k2,11),ir(k2,11),irct,ird * / cst88 /vu(2,k1),jprct,jmct,jmyn * / cst1 /thermo(k4,k1),vf(2,k1),vs(h5,k1),uf(2),us(h5) c----------------------------------------------------------------------- c save dependent extensities c it is necessary to calculate c the initial G cause iptol may c be large relative to the c change in G due to delt og = grxn(ird) do 1 i = 1, 2 v(i) = v(i) + delt(i) dv(i) = (og - grxn(i)) / delt(i) 1 v(i) = v(i) - delt(i) c get deltas on mobile comps if (jmct.gt.0) then do 40 i = 1, jmct 40 du(i) = 0. do 9 i = 1, ivct do 8 j = 1, jmct du(j) = du(j) + vnu(i) * vu(j,idr(i)) 8 continue 9 continue end if c do inverse projection to c get stoichiometries of c constrained components: if (isat.gt.0) then if (isat.gt.1) then do 3 j = 1, isat - 1 do 22 i = j + 1, isat 22 vus(j) = vus(j) - vus(i) * vs(j,idss(i)) 3 continue end if c now do the effects of c the inverse projection on c saturated phase components: if (ifct.ne.0) then do 5 i = 1, isat do 4 j = 1, ifct 4 vuf(j) = vuf(j) - vus(i) * vf(j,idss(i)) 5 continue end if c now do the effects of c the inverse projection on c mobile components: if (jmct.ne.0) then do 6 i = 1, isat do 7 j = 1, jmct 7 du(j) = du(j) - vus(i) * vu(j,idss(i)) 6 continue end if end if c see if coefficients should c change sign: invert = 1 if (ivfl.eq.1.and.dv(1).gt.0.) then invert = 0 else if (ivfl.eq.2.and.dv(2).gt.0.) then invert = 0 else if (ivfl.gt.2) then if (du(ivfl-3).gt.0.d0) then c diagram must have a mu variable invert = 0 else if (du(ivfl-3).eq.0.d0) then if (iv(2).eq.3) then if (vuf(2).lt.0.d0) invert = 0 else if (iv(1).eq.3) then if (vuf(2).lt.0.d0) invert = 0 else if (du(iv(1)-3).gt.0.d0) then invert = 0 end if end if end if if (invert.eq.1) then do 10 i = 1, ivct vnu(i) = -vnu(i) c why a different signs? 10 vn(ird,i) = vnu(i) do 20 i = 1, isat 20 vus(i) = -vus(i) do 25 i = 1, ifct 25 vuf(i) = -vuf(i) do 30 i = 1, jmct 30 du(i) = -du(i) dv(1) = -dv(1) dv(2) = -dv(2) end if c save id's of saturated comps if (isyn.eq.0) then do 2 i = 1, isat 2 jds(i) = idss(i) end if end function grxn (i) c----------------------------------------------------------------------- c grxn computes the free energy of univariant equilibria c defined by the data in commonn block cst21 which is initialized c in the subprogram balanc. grxn is partially redundant with c the function gphase but because of the frequency that these c these routines are used a significant increase in efficiency is c gained by maintaining separate functions. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (h5=7) common/ cst5 /p,t,xco2,u1,u2,tr,pr,r,ps * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn/ cst11 /fh2o,fco2 * / cst201 /vuf(2),vus(h5),iffr,isr * / cst6 /icomp,istct,iphct,icp/ cst25 /vnu(11),idr(11),ivct c----------------------------------------------------------------------- c compute potentials of saturated phases c and components, note that in this c version of vertex the stoichiometry of c such components may vary. c no saturated phase components and no c saturated components: if (iffr.eq.1) goto (10),isyn c note that this call to uproj makes a c subsequent call in gall redundant if c sfol1 is used to trace a univariant c curve. call uproj c compute free energy change of the rxn 10 grxn= 0.d0 do 20 j = 1, ivct 20 grxn = grxn + vnu(j) * gproj(idr(j)) end subroutine header c----------------------------------------------------------------------- c header writes the graphics file header for schreinemakers diagrams c to lun n4 c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k2=1200,k3=1200,k6=7,l2=5) parameter (j9=1200,h5=7,h6=500,h9=18) character*8 names,vname,gname*78,xname*18,title*72,fname*10 character*5 cname common/ cst8 /names(k1)/ cst6 /icomp,istct,iphct,icp/ cst79 /iso * / csta2 /xname(k6),vname(l2),gname(l2)/ cst83 /ig(l2) * / cst60 /iasmbl(j9),ipoint,imyn/ cst34 /ip1(k3),iipct,loop * / cst9 /vmax(l2),vmin(l2),dv(l2)/ csta7 /fname(h9) * / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5 / csta8 /title * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar * / cst102 /ivfl/ cst40 /ids(h5,h6),isct(h5),icp1,isat,io2 * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn/ csta4 /cname(12) c----------------------------------------------------------------------- c value to be read as icopt write (n4,1020) 1 c write phase and solution count: write (n4,1020) iphct, iso c write phase names: write (n4,1010) (names(i),i = 1, iphct) c write solution phase flags: write (n4,1060) (ikp(i),i = 1, iphct) c write solution names: if (iso.ne.0) write (n4,1040) (fname(i),i = 1, iso) c for ens graphics: c write (n4,*) iset c write (n4,1220) loop,iv1,iv2,iv3,vmax(iv1),vmax(iv2),vmax(iv3) c * ,vmin(iv1),vmin(iv2),vmin(iv3) c write (n4,1000) vname(iv1),vname(iv2),vname(iv3) c for di3000 graphics write (n4,1150) title c write saturated and buffered c component names: if (isyn.eq.0) then write (n4,1070) (cname(z+icp),z= 1, isat) else write (n4,1150) ' ' end if write (n4,1080) vname(ivfl) write (n4,1020) ipot write (n4,1030) (vmax(jv(z)),vmin(jv(z)),z= 1, ipot) write (n4,1000) (ig(jv(z)),gname(jv(z)),z= 1, ipot) 1000 format (i2,a78) 1010 format (10a8) 1020 format (20(i3,1x)) 1040 format (8a10) 1030 format (6(d11.5,1x)) 1060 format (26(i2,1x)) 1070 format ('Component saturation hierarchy: ',7(a5,1x)) 1080 format ('Reaction equations are written such that the high ', * a8,/,'assemblage is on the right of the = sign') 1150 format (a72) end function lchk (lphi) c----------------------------------------------------------------------- c a function subprogram to determine if a phase, lphi, lies c below a g-x plane defined by the vertices idv. c lchk has a value of 0 if the plane is stable, and 1 if not. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k6=7,l2=5) common/ cst23 /a(10,10),b(11),ipvt(10),idv(10), * iophi,idphi,iiphi,iflg1 * / cst62 /blim(l2),ulim(l2),dg/ cst12 /cp(k6,k1) * / cst6 /icomp,istct,iphct,icp/ cst2 /g(k1) * / cst87 /delt(l2),dtol,utol,ptol c----------------------------------------------------------------------- call uproj do 10 i = 1, icp 10 b(i) =gproj(idv(i)) g(lphi) =gproj(lphi) lchk= 0 call subst (a,ipvt,icp,b) gphi = 0.d0 do 30 j = 1, icp 30 gphi = gphi + cp(j,lphi) * b(j) if (gphi.lt.g(lphi)) goto 99 lchk= 1 99 end subroutine newass (b,idv,idphi,ivd,ivi,iste) c----------------------------------------------------------------------- c newass determines the (pseudo-) divariant assemblages generated by c a reaction based on the stoichiometric coefficients of the phases c in the reaction as determined by routine balanc. c input: idv(icp) - a vector containing the indices of phases in a c divariant assemblage. c idphi - the index of the phase which is involved with the c the divariant assemblage (idv) in a univariant equilibrium. c b(icp) - the stoichiometric coefficient of the phases inedxed c by idv in the reaction of these phases to form the phase c idphi. by convention the stoichiometric coeffeicient of c idphi is -1.0. c output: jdv(icp) - a vector containing the indices of phases in a c divariant assemblage generated by the reaction between the c phases idphi and idv this vector is used by the routine c chkass. this assemblage is stable on the opposite 'side' of c the univariant equilibrium from the assemblage idv. c referenced by: coface c references to: chkass c algorithm: phases with null or negative coefficients are always c compatible with eachother on the 'side' of a univariant reaction c which idv is metastable. the k phases with positive coefficients c define a chemographic simplex surrounding phases with null or c negative coefficients. consequently a univariant reaction will c generate k new divariant assemblages consisting of k-1 of the phases c in the original assemblage (idv) and all the phases with null or c negative coefficients. if only one phase has a positive coefficient c this phase is unstable and only one new assemblage will be generated c which consists of all the other phases (including idphi). c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k2=1200,k3=1200,k6=7) character*8 names dimension idv(10),b(11),idpos(7),jdv(k6) common/ cst6 /icomp,istct,iphct,icp * / cst8 /names(k1) * / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar common/ cst89 /i2c(2,21),i3c(3,35),i4c(4,35),i5c(5,21),i6c(6,7) * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 c----------------------------------------------------------------------- c initialize counter: ipos = 0 do 10 i = 1, icp c is the b(i) coefficient zero/negative? if (b(i).lt.1.d-8) goto 20 c no, save in the array idpos and count. ipos=ipos+1 idpos(ipos) =idv(i) goto 10 c yes, load into the array jdv 20 jdv(i-ipos) =idv(i) 10 continue c load idphi ineg=icp-ipos+1 jdv(ineg) =idphi ineg=ineg+1 c generate the divariant assemblages do 30 j = 1, ipos goto (30,40,50,60,70,80,90),ipos c ipos=2, 2 assemblages generated: 40 jdv(ineg) =idpos(j) goto 30 c ipos=3, 3 assemblages generated: 50 jdv(ineg) =idpos(i2c(1,j)) jdv(ineg+1) =idpos(i2c(2,j)) goto 30 c ipos=4, 4 assemblages generated: 60 jdv(ineg) =idpos(i3c(1,j)) jdv(ineg+1) =idpos(i3c(2,j)) jdv(ineg+2) =idpos(i3c(3,j)) goto 30 c ipos=5, 5 assemblages generated: 70 jdv(ineg) =idpos(i4c(1,j)) jdv(ineg+1) =idpos(i4c(2,j)) jdv(ineg+2) =idpos(i4c(3,j)) jdv(ineg+3) =idpos(i4c(4,j)) goto 30 c ipos=6, ... 80 itic = 0 do 1 k = ineg, ineg + 4 itic = itic + 1 1 jdv(ineg) = idpos(i5c(itic,j)) goto 30 c ipos=7, ... 90 itic = 0 do 2 k = ineg, ineg + 5 itic = itic + 1 2 jdv(ineg) = idpos(i6c(itic,j)) c send to the routine chkass to store th c the assemblage if unique. 30 call chkass (jdv,ivd,ivi,iste) c done end subroutine onedim c---------------------------------------------------------------------- c for phase diagrams with one extrinsic variable (iv1) onedim sorts c equilibria, identified by newhld and located by coface, by veq(iv1). c the sorted equilibria are then output to units n3 and n4 as requested. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k2=1200,k3=1200,k8=9,l2=5) character*8 names,irnms common/ cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar * / cst29 /vip(l2,k3),ipid(k3,k8),ipct * / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5 * / cst31 /vn(k2,11),ir(k2,11),irct,ird/ cst8 /names(k1) * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst35 /ir1(k2),iirct,irpct,irv(k2)/ csta1 /irnms(k2,11) c----------------------------------------------------------------------- c any rxns to be sorted? if (irct.eq.1) goto 30 c initialize: ist = 2 c begin sort loop: 20 ist1 = ist - 1 itemp = ir1(ist1) vst = vip(iv1,itemp) do 10 i = ist, irct if (vip(iv1,ir1(i)).gt.vst) goto 10 ir1(ist1) = ir1(i) ir1(i) = itemp itemp = ir1(ist1) vst = vip(iv1,itemp) 10 continue ist = ist + 1 if (ist.gt.irct) goto 30 goto 20 c end of sort loop 30 if (irct.eq.0) goto 99 c output invariant reactions and conditions: call outirn 99 end subroutine outier c---------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k2=1200,k3=1200,k6=7,k8=9,l2=5) character*8 ipnms, irnms, vname*8, gname*78, xname*18, * text(256)*1, names common/ cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar * / cst29 /vip(l2,k3),ipid(k3,k8),ipct * / cst24 /ipot,jv(l2),iv(l2)/ cst8 /names(k1) * / csta2 /xname(k6),vname(l2),gname(l2)/ cst83 /ig(l2) * / cst34 /ip1(k3),iipct,loop/ csta3 /ipnms(k3,k8) * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst81 /ipchk(k3),ip2(k3),icp2,ifake * / cst103 /isec,icopt,ifull,imsg,io3p * / cst35 /ir1(k2),iirct,irpct,irv(k2)/ csta1 /irnms(k2,11) c----------------------------------------------------------------------- c output for each tier goto (65),io4 c write end of ptx data for univariant c equilibria marker to graphics file. write (n4,1070) 1,1,1,1,1,1. c write ip names and locations to c units n6 and n3. write (n6,1070) ipct,icp2 65 if (ipct.eq.0) goto 99 goto (55),io3 write (n3,1050) write (n3,1030) 55 do 30 j = 1, ipct goto (75), io3 goto (75), io3p c n3: call iptext (text,iend,j,icp2) write (n3,1080) ip1(j),ivarip(j),(text(z),z= 1, iend) write (n3,1020) write (n3,1040) (vname(jv(z)),vip(jv(z),j),z= 1, ipot) 75 goto (85),io4 c n4: write (n6,1000) ip1(j),ivarip(j),(ipid(ip1(j),i),i = 1, icp2) write (n6,1060) (vip(jv(z),j),z= 1, ipot) c check if the ip was encountered in c an earlier tier if not save id in ip2 c for the cumulative list 85 if (ipchk(ip1(j)).eq.1) goto 30 ipchk(ip1(j)) =1 irpct=irpct + 1 ip2(irpct) =ip1(j) 30 continue 1000 format (i5,1x,i1,9(1x,i3)) 1020 format (15x,'occurs at:') 1030 format (' (pseudo-) invariant points are summarized below:') 1040 format (25x,a8,'=',g12.6) 1050 format (/,'------------------------------------------------------' * ,'-----------------',/) 1060 format (5(1x,g11.5)) 1070 format (i5,1x,4(i3,1x),f3.1) 1080 format (/,' (',i4,'-',i1,') ',200a1) 99 end subroutine outirn c----------------------------------------------------------------------- c outirn writes the identity and equilibrium conditions of invariant c equilibria to units n3 and n4 after they have been sorted by onedim. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k2=1200,k3=1200,k6=7,k8=9,l2=5) parameter (m1=21,j1=50,h9=18) character*8 irnms, vname, gname*78, xname*18, rxnstr*132, * fname*10 common/ cst31 /vn(k2,11),ir(k2,11),irct,ird/ cst3 /x(k1,k6) * / cst9 /vmax(l2),vmin(l2),dv(l2)/ csta7 /fname(h9) * / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5/ cst79 /iso * / cst6 /icomp,istct,iphct,icp/ cst104 /rxnstr(k2) * / cst15 /idc(k6),idbv(m1,j1),ibvct(m1) * / cst29 /vip(l2,k3),ipid(k3,k8),ipct/ cst5 /v(l2),q(4) * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst35 /ir1(k2),iirct,irpct,irv(k2)/ csta1 /irnms(k2,11) * / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar * / csta2 /xname(k6),vname(l2),gname(l2)/ cst83 /ig(l2) c----------------------------------------------------------------------- c write header for graphics file c currently setup only for binary systems: goto (5), io4 c write header info: call outgrf c write graphics code names for x variables: write (n4,1010) (xname(i), i = 2, icp) c write solution phase flags: write (n4,1020) (ikp(i), i = 1, iphct) if (iso.ne.0) write (n4,1080) (fname(i),i = 1, iso) c write min and max values for the c independent potentials: write (n4,*) vmin(iv1),vmax(iv1) c write invariant reactions: 5 do 50 i = 1, irct j = ir1(i) c the index on these two guys used to be j c must be wrong? if so why did it work c before? changed to i. you know sometimes c i'm so stupid, i just can't believe it. k = irv(j) l = ivarrx(j) goto (10),io3 c output to print file: if (l.eq.1) then write (n3,1030) j,l,rxnstr(j) else write (n3,1040) j,l,rxnstr(j) end if write (n3,1070) vname(iv1),vip(iv1,j),vname(iv2),vip(iv2,j) 10 write (n3,1160) c output to graphics file goto (50), io4 write (n4,1000) j,k,l,vip(iv1,j),(ir(j,zl),zl= 1,k) write (n4,1050) (vn(j,zl), zl = 1, k) 50 continue 1000 format (i5,1x,2(i3,1x),g12.5,1x,6(i3,1x)) 1010 format (5(a18)) 1020 format (20(i3,1x)) 1030 format (' The equilibrium of the invariant reaction:',//, * ' (',i4,'-',i1,') ',a123) 1040 format (' The equilibrium of the pseudoinvariant reaction:',//, * ' (',i4,'-',i1,') ',a123) 1050 format (12x,6(g9.2,1x)) 1070 format (/,' occurs at ',a8,'=',g12.6,' and ',a8,'=',g12.6) 1160 format (' ----------------------------------------') 1080 format (8a10) end subroutine outlst c---------------------------------------------------------------------- c outlst writes cumulative lists of the univariant and invariant c equilibria identified during a single execution of the program c to the graphics and print files. outlst also calls the subroutine c out5 to write a permanent list (on n5) of equilibria for multiple c executions. c---------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k2=1200,k3=1200,k8=9,l2=5) character*8 names,ipnms,irnms,cname*5,rxnstr*132 integer io(2), n(2) common/ cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar * / cst81 /ipchk(k3),ip2(k3),icp2,ifake * / cst31 /vn(k2,11),ir(k2,11),irct,ird/ cst105 /ibug(k2) * / cst8 /names(k1)/ cst29 /vip(l2,k3),ipid(k3,k8),ipct * / cst34 /ip1(k3),iipct,loop/ csta3 /ipnms(k3,k8) * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst35 /ir1(k2),iirct,irpct,irv(k2)/ csta1 /irnms(k2,11) * / cst37 /iprct,ixct,iexyn,istbyn/ csta4 /cname(12) * / cst80 /jrchk(k2),irchk(k2)/ cst104 /rxnstr(k2) c------------------------------------------------------------------------ c printing to n4 is suppresses for the c erlanger system io(1) =io3 io(2) =io4 n(1) =n3 n(2) =n4 jbug = 0 c invariant point lists: if (irpct.eq.0) goto 80 goto (20),io4 c output graphics file heading: c write (n4,*) irpct+2 c output print file heading: c output ip lists: 20 do 60 j =1,1 if (io(j).eq.1) goto 60 nn=n(j) write (nn,1030) do 70 i = 1, irpct c call iptext (text,iend,i,icp2) c write (nn,1000) ip2(i),ivarip(i),(text(k),k= 1, iend) 70 write (nn,1150) ip2(i),ivarip(i),(ipnms(ip2(i),zj),zj = 1, icp2) 60 continue c write end of invariant point list c marker to graphics file. 80 goto (110),io3 write (n3,1160) itic= 0 c 110 do 140 j = 1, 1 if (io(j).eq.1) goto 140 nn=n(j) c write list of stable univariant c equilibria to units n3 and n4. write (nn,1070) do 120 i = 1, iirct if (jrchk(i).eq.0) goto 120 if (ibug(i).ne.0) jbug = 1 itic=itic+1 120 write (nn,1000) i,ivarrx(i),rxnstr(i) c if itic= 0 no equilibria were identified c goto 9000 to write no eq mess. if (itic.eq.0) goto 9000 write (nn,1160) 140 continue c jbug ne 0, one or more equilibria c may not have completely defined, write c list of possible failures: if (jbug.ne.0) then write (n3,1160) write (*,1160) write (n3,1190) write (*,1190) do 141 i = 1, iirct if (ibug(i).eq.0) goto 141 write (*,1000) i,ivarrx(i),rxnstr(i) write (n3,1000) i,ivarrx(i),rxnstr(i) 141 continue write (n3,1160) write (*,1160) end if c write permanent lists if requested c---------------------------------------------------------------------- c write to unit n5. this file contains a permanent list c of univariant and invariant equilibria. see comments for input5. c----------------------------------------------------------------------- goto (99),io5 rewind n5 c write cumulative list of invariant eq: write (n5,*) iipct if (iipct.eq.0) goto 150 write (n5,1060) ((ipnms(zi,zj),zj = 1, icp2),zi = 1, iipct) 150 write (n5,*) iirct if (iirct.eq.0) goto 99 c write list of univariant reactions do 160 i = 1, iirct ip=irv(i) write (n5,*) ip 160 write (n5,1060) (irnms(i,zj),zj = 1, ip) goto 99 c write no equilibrium messages: 9000 goto (100),io3 write (n3,1180) 100 goto (99),io4 write (n4,1180) 1000 format (' (',i4,'-',i1,') ',a72) 1030 format (' (pseudo-) invariant points are summarized below:',/) 1060 format (9(a8,1x)) 1070 format (' (pseudo-) univariant equilibria are summarized ', * 'below:',/) 1150 format (' (',i4,'-',i1,') ',9(a8,1x)) 1160 format (//,'-------------------------------------------------' * ,'---------------',//) 1180 format (' no univariant or invariant equilibria occur.') 1190 format (' WARNING!! The stability fields of the following', * ' equilibria may',/,' have been entirely or', * ' partially skipped in the calculation: ',/) 99 end subroutine iptext (text,iend,jp,icp2) c----------------------------------------------------------------------- c subprogram to write a text label for an invariant point c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k2=1200,k3=1200,k8=9,l2=5,h9=18) character*8 names, text(256)*1, fname*10 common/ cst8 /names(k1)/ cst34 /ip1(k3),iipct,loop * / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar * / cst29 /vip(l2,k3),ipid(k3,k8),ipct/ csta7 /fname(h9) c dump the names into the array text ist = 1 do 20 i = 1, icp2 id = ipid(jp,i) if (ikp(id).eq.0) then c simple compound: read (names(id),1000) (text(j), j = ist, ist + 7) ist = ist + 8 else c solution phases: read (fname(ikp(id)),1000) (text(j), j = ist, ist + 9) text(ist + 10) = '(' read (names(id),1000) (text(j), j = ist + 11, ist + 18) text(ist + 19) = ')' ist = ist + 20 end if text(ist) = ' ' 20 ist = ist + 1 c now filter out blanks iend = 1 do 40 i = 2, ist - 1 if (text(i).eq.' '.and.text(i + 1).eq.' ') then goto 40 else if (text(i).eq.' '.and.text(i + 1).eq.')') then goto 40 else if (text(i).eq.' '.and.text(i + 1).eq.'(') then goto 40 end if iend = iend + 1 text(iend) = text(i) 40 continue 1000 format (20a1) end subroutine rxntxt (iend,jend,text,alpha) c----------------------------------------------------------------------- c subprogram to write a text label for a reaction c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k2=1200,k3=1200,h9=18) character*8 names, text(256)*1, fname*10, * alpha(90)*1, string*132, rxnstr*132 integer iplus(8), iminus(8), jplus(8), jminus(8), * ione(8), jone(8) common/ csta7 /fname(h9)/ cst8 /names(k1) * / cst25 /vnu(11),idr(11),ivct * / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar * / cst31 /vn(k2,11),ir(k2,11),irct,ird * / cst80 /jrchk(k2),irchk(k2)/ cst104 /rxnstr(k2) ip = 0 im = 0 jend = 1 c load id's of phases with + or - coeffs do 10 i = 1, ivct if (vnu(i).gt.0.0) then ip = ip + 1 iplus(ip) = idr(i) jplus(ip) = i else im = im + 1 iminus(im) = idr(i) jminus(im) = i end if 10 continue do 60 i = 1, im jone(i) = jminus(i) 60 ione(i) = iminus(i) do 70 i = 1, ip j = i + im jone(j) = jplus(i) 70 ione(j) = iplus(i) c dump coefficients into string: write (string,1010) (vnu(jone(i)),i = 1, ivct) c read alpha from string: jst = 7 + ivct * 10 read (string,1000) (alpha(i), i = 1, jst) c filter blanks and add right ) do 50 i = 1, jst if (alpha(i).eq.' '.and.alpha(i + 1).eq.' ') then goto 50 else if (alpha(i).eq.' '.and.alpha(i + 1).eq.',') then goto 50 end if alpha(jend) = alpha(i) jend = jend + 1 50 continue if (alpha(jend-1).eq.' '.or.alpha(jend-1).eq.',') * jend = jend - 1 alpha(jend) = ')' c now dump the names into the array text ist = 1 is = 1 80 do 20 i = is, im id = ione(i) if (ikp(id).eq.0) then c simple compound: read (names(id),1000) (text(j), j = ist, ist + 7) text(ist+8) = ' ' ist = ist + 9 else c solution phases: read (fname(ikp(id)),1000) (text(j), j = ist, ist + 9) text(ist + 10) = '(' read (names(id),1000) (text(j), j = ist + 11, ist + 18) text(ist + 19) = ')' text(ist+20) = ' ' ist = ist + 21 end if 20 continue if (is.eq.1) then text(ist) = '=' text(ist + 1) = ' ' ist = ist + 2 is = im + 1 im = ivct goto 80 else text(ist) = ' ' end if c now filter out blanks iend = 1 do 40 i = 2, ist if (text(i).eq.' '.and.text(i + 1).eq.' ') then goto 40 else if (text(i).eq.' '.and.text(i + 1).eq.')') then goto 40 else if (text(i).eq.' '.and.text(i + 1).eq.'(') then goto 40 end if iend = iend + 1 text(iend) = text(i) 40 continue c this is completely wrong c if old lists are being used, c i.e. iirct ne irct. if (jrchk(ird).eq.0) then jrchk(ird) = 1 if (iend.gt.72) then write (rxnstr(ird),1000) (text(i),i=1,72) else write (rxnstr(ird),1000) (text(i),i=1,iend) end if end if if (iend.gt.200) iend = 200 1000 format (132a1) 1010 format ('Alpha(',8(g9.3,', ')) end subroutine outrxn (ip,ier) c----------------------------------------------------------------------- c outrxn writes the identity and v1-v2 loci of univariant equilibria c to units n3 and n4. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k2=1200,k3=1200,h9=18) character*8 names, text(256)*1, alpha(90)*1, fname*10 common/ cst31 /vn(k2,11),ir(k2,11),irct,ird/ cst7 /iflag * / cst25 /vnu(11),idr(11),ivct/ cst8 /names(k1) * / cst32 /ptx(450),logx,ipt2/ csta7 /fname(h9) * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst35 /ir1(k2),iirct,irpct,irv(k2) * / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar * / cst105 /ibug(k2)/ cst103 /isec,icopt,ifull,imsg,io3p c set bug flag to indicate that c the equilibrium was at least c partially located. if (ier.lt.2) ibug(iirct) = 0 if (ifull.eq.0) then call rxntxt (iend,jend,text,alpha) else call fultxt (iend,text) end if if (icopt.eq.3) goto 99 if (imsg.eq.0) write (*,1110) ird,(text(i),i = 1, iend) goto (10), io3p write (n3,1120) ird,ivarrx(ird),(text(i),i = 1, iend) if (ifull.ne.0) goto 80 write (n3,1130) (alpha(i),i = 1,jend) if (ipt2.lt.3) then write (n3,*) goto 99 end if call outdel 80 write (n3,*) write (n3,1000) (ptx(i),i = 1, ipt2) write (n3,*) if (ier.ne.0) goto 10 goto (40),iflag goto 10 40 write (n3,1070) ip write (n3,1040) 10 goto (99), io4 write (n4,1080) ipt2,ird,ivar,ivct,(idr(i),i = 1, ivct) write (n4,1090) (vnu(i),i = 1, ivct) write (n4,1000) (ptx(i),i = 1, ipt2) 1000 format (1x,g12.6,1x,g12.6,2x,g12.6,1x,g12.6,2x,g12.6,1x,g12.6) 1040 format (/) 1070 format (' the equilibrium extends to invariant point (',i4,')') 1080 format (16(i4,1x)) 1090 format (6(g12.4,1x)) 1110 format (' finished with equilibrium (',i4,') ',200a1) 1120 format (' (',i4,'-',i1,') ',200a1) 1130 format (/,10x,90a1) 99 end subroutine pchk (igo) c----------------------------------------------------------------------- c a subprogram which returns the difference in free energy c beteen a g-x plane, defined by the assemblage idv, and a phase lphi. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k6=7,l2=5) common/ cst23 /a(10,10),b(11),ipvt(10),idv(10), * iophi,idphi,iiphi,iflg1 * / cst6 /icomp,istct,iphct,icp/ cst2 /g(k1)/ cst7 /iflag * / cst12 /cp(k6,k1)/ cst87 /delt(l2),dtol,utol,ptol c----------------------------------------------------------------------- c compute energies: igo = 0 call uproj do 10 i = 1, icp 10 b(i) =gproj(idv(i)) dg=gproj(idphi) c solve for chemical potentials: call subst (a,ipvt,icp,b) c compute energy difference: do 30 j = 1, icp 30 dg = dg - cp(j,idphi)*b(j) if (dabs(dg).lt.ptol) then igo= 1 iflag= 1 else if (dg.gt.dtol) then iflag= 0 else iflag= 1 end if end if end subroutine pivots (ier) c-------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k6=7) common/ cst23 /a(10,10),b(11),ipvt(10),idv(10), * iophi,idphi,iiphi,iflg1 * / cst6 /icomp,istct,iphct,icp/ cst12 /cp(k6,k1) c-------------------------------------------------------------------- c calculate pivots for the matrix 'a' do 1 k = 1, icp do 2 j = 1, icp 2 a(k,j) = cp(j,idv(k)) 1 continue call factor (a,icp,ipvt,ier) end subroutine schk (lphi) c----------------------------------------------------------------------- c schk my ___ c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k6=7,l2=5) common/ cst23 /a(10,10),u(11),ipvt(10),idv(10), * iophi,idphi,iiphi,iflg1 * / cst6 /icomp,istct,iphct,icp/ cst62 /blim(l2),ulim(l2),dg * / cst87 /delt(l2),dtol,utol,ptol * / cst2 /g(k1)/ cst12 /cp(k6,k1)/ cst7 /iflag c----------------------------------------------------------------------- c determine chemical potentials iflag= 0 do 10 i = 1, icp 10 u(i) = g(idv(i)) call subst (a,ipvt,icp,u) c test phases, not=iophi, against the c assemblage idv(i),i = 1, icp do 20 i = istct, iphct if (i.eq.iophi.or.i.eq.lphi) goto 20 dg = g(i) do 30 j = 1, icp 30 dg = dg - cp(j,i) * u(j) if (dg.gt.dtol) goto 20 c check that a phase is not metastable c with respect to itself, this c could be remedied by changing the c value of dtol. do 25 j = 1, icp 25 if (idv(j).eq.i) goto 20 iflag = iflag + 1 idphi = i goto (20,40,40), iflag 20 continue 40 end subroutine nschk c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k6=7,l2=5) common/ cst23 /a(10,10),u(11),ipvt(10),idv(10), * iophi,idphi,iiphi,iflg1 * / cst6 /icomp,istct,iphct,icp/ cst62 /blim(l2),ulim(l2),dg * / cst87 /delt(l2),dtol,utol,ptol * / cst2 /g(k1)/ cst12 /cp(k6,k1)/ cst7 /iflag c----------------------------------------------------------------------- c determine chemical potentials iflag= 0 do 10 i = 1, icp 10 u(i) = g(idv(i)) call subst (a,ipvt,icp,u) c test phases, not=iophi, against the c assemblage idv(i),i = 1, icp do 20 i = istct, iphct if (i.eq.iophi) goto 20 dg = g(i) do 30 j = 1, icp 30 dg = dg - cp(j,i) * u(j) if (dg.gt.dtol) goto 20 c check that a phase is not metastable c with respect to itself, this c could be remedied by changing the c value of dtol. do 25 j = 1, icp 25 if (idv(j).eq.i) goto 20 iflag = iflag + 1 idphi = i goto (20,40,40), iflag 20 continue 40 end subroutine search (vst,vsti,ist,iste,ivi,ivd,div,ier) c----------------------------------------------------------------------- c search looks for an equilibrium around the periphery of the coordinate c frame specified for a spd calculation. if a stable equilibrium is c found ier = 0, if an error occurs ier =2, and if no equilibrium is c found ier =1. c---------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,l2=5) character*8 names common/ cst23 /a(10,10),b(11),ipvt(10),idv(10), * iophi,idphi,iiphi,iflg1 * / cst5 /v(l2),x(4)/ cst7 /iflag * / cst6 /icomp,istct,iphct,icp * / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5/ cst38 /iedge * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst9 /vmax(l2),vmin(l2),dv(l2)/ cst8 /names(k1) c----------------------------------------------------------------------- c initialization iflg1 = 0 ier = 0 iugh = 0 v(iv1) = vst v(iv2) = vst do 60 i = ist, iedge iste = i c set default dependent and independent c variables and increments for each edge c to be searched. goto (10,20,30,40),i c traverse 1. 10 ivi = iv2 ivd = iv1 ddv = dv(ivd) div = dv(ivi) v(ivi) = vmin(ivi) goto 80 c traverse 2. 20 ivi = iv1 ivd = iv2 ddv = dv(ivd) div = -dv(ivi) v(ivi) = vmax(ivi) goto 80 c traverse 3. 30 ivi =iv2 ivd=iv1 ddv = -dv(ivd) div = -dv(ivi) v(ivi) = vmax(ivi) goto 80 c traverse 4. 40 ivi =iv1 ivd=iv2 ddv = -dv(ivd) div = dv(ivi) v(ivi) = vmin(ivi) c in some cases newass may generate c a metastable assemblage call flipit c to check. if the assemblage is c stable flipit will return jer = 0, c if it is metastable flipit tries a c reverse traverse and returns jer =1 c if it is succesful, and jer =2 if not. 80 if (i.eq.ist) then v(ivd) = vst v(ivi) = vsti end if call flipit (ddv,vst,ivd,ist,i,jer) if (jer.eq.2) then write (*,1000) (names(idv(zk)),zk= 1, icp) write (*,1010) v ier = 1 goto 99 end if vlast = v(ivd) goto (99),jer c begin search: 100 v(ivd) = v(ivd) + ddv c out of range?: goto (110,110,120,120),i 110 if (v(ivd).gt.vmax(ivd)) then v(ivd) = vmax(ivd) else if (i.eq.ist) then if (v(ivd).gt.vst) goto 130 ddv = dabs(ddv)/2.d0 v(ivd) = vst iflg1= 0 goto 100 end if goto 130 120 if (v(ivd).lt.vmin(ivd)) then v(ivd) = vmin(ivd) else if (i.eq.ist) then if (v(ivd).lt.vst) goto 130 ddv = -dabs(ddv)/2.d0 v(ivd) = vst iflg1= 0 goto 100 end if c calculate phase energies: 130 call gall (istct,iphct) c test stability of the assemblage: call asschk if (iflag.eq.0) vlast = v(ivd) c iflag= 0 still stable c iflag=1 metastable with respect to one c phase indexed by idphi. c iflag=2 metastable with respect to c multiple phases refine the c search increment. goto (99),iflag c check if search is in range: goto (140,140,150,150),i 140 if (v(ivd).ge.vmax(ivd).and.iflag.eq.0) goto 60 goto 160 150 if (v(ivd).le.vmin(ivd).and.iflag.eq.0) goto 60 c refine increment if necessary: 160 call delvar (ddv,iflag,iflg1) c check if increment has been refined c beyond the tolerance which is c arbitrarily set at 1.d-8. if (dabs(ddv).lt.1.d-8) then iugh = iugh + 1 if (i.lt.3) then c traverses 1,2. ddv = dv(ivd) else c traverses 3,4. ddv = -dv(ivd) end if v(ivd) = vlast ddv = ddv/10. iflg1 = 0 if (iugh.gt.3) then write (n8,1020) ivd,dv(ivd),(names(idv(j)),j = 1, icp) write (n8,*) ' ' ier = 2 goto 99 end if goto 100 end if c increment conditions again: goto 100 c next traverse. 60 continue ier = 1 1000 format (' the assemblage is:',4x,7(1x,a8)) 1010 format (' v =',5(g12.6,1x)) 1020 format (/,' **warning ver047** SEARCH: > 1 equilibrium', * ' occurs within the minimum search',/,' increment for', * ' variable: ',i1,'. This often occurs as YCO2 =>', * ' 1 or => 0.',/,' You may be able to correct this by', * ' assigning slightly different limits',/,' for YCO2 or', * ' by reducing the default increment for this variable', * ' (',g12.3,')',/,' in the computational option file.',/, * ' Equilibria involving the following assemblage may', * ' not be delineated:',/,7(1x,a8)) 99 end subroutine sfol1 (ivd,ivi,ier,dv,lphi,ikwk,irend) c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k2=1200,k3=1200,k8=9,l2=5) common/ cst23 /a(10,10),b(11),ipvt(10),idv(10), * iophi,idphi,iiphi,iflg1 * / cst5 /v(l2),tr,pr,r,ps/ cst80 /jrchk(k2),irchk(k2) * / cst7 /iflag/ cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5 * / cst9 /vmax(l2),vmin(l2),ddv(l2)/ cst63 /delv(l2) * / cst6 /icomp,istct,iphct,icp/ cst32 /ptx(450),logx,ipt2 * / cst29 /vip(l2,k3),ipid(k3,k8),ipct * / cst31 /vn(k2,11),ir(k2,11),irct,ird * / cst87 /delt(l2),dtol,utol,ptol * / cst62 /blim(l2),ulim(l2),dg * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 c----------------------------------------------------------------------- if (ikwk.eq.1) then irchk(ird) = 1 goto 9000 end if c begin traverse: 10 v(ivi) = v(ivi) + dv c is search in range? if (v(ivi).gt.vmax(ivi)) then v(ivi) = vmax(ivi) else if (v(ivi).lt.vmin(ivi)) then v(ivi) = vmin(ivi) end if c solve for the equilibrium conditions: call univeq (ivd,ier) c on error return: c calling routine will switch variables. goto (999,999),ier c test if the equilibrium is stable: call gall (istct,iphct) call nschk c iflag= 0 stable. c iflag=1 metastable w/respect a phase. c iflag=2 metastable w/respect to more c than one phase. goto (20,70),iflag c iflag= 0: if (ipt2.gt.449) goto 9000 c dependent v in range? if c greater than the maximum value for v c or less than the minimum value for v c reset conditions, and c switch independent/dependent variables if (v(ivd).gt.vmax(ivd)) then v(ivd) = vmax(ivd) else if (v(ivd).lt.vmin(ivd)) then v(ivd) = vmin(ivd) else call assptx c if independent variable has reached a c limit of the diagram output traverse c otherwise increment the variable and c continue the traverse if ((v(ivi).eq.vmax(ivi)).or.(v(ivi).eq.vmin(ivi))) then c flag irchk indicates that the c traverse ended on an edge. irchk(ird) = 1 goto 9000 end if c try for next point goto 10 c no route out of this construction end if c solve for the equilibrium with c the switched variables: call univeq (ivi,ier) c flag irchk indicates that the c traverse ended on an edge. irchk(ird) = 1 c on error output traverse anyway goto (9000,9000),ier c assign final value call assptx goto 9000 c iflag=2 reset to last stable condition c and redefine increment dv 70 v(ivi) = v(ivi) - dv v(ivd) = v(ivd) dv = dv/2.d0 if (dabs(dv).gt.delt(ivi)) goto 10 call warn (73,delt(ivi),ivi,'SFOL2 ') iflag = 0 goto 9000 c iflag=1 (an invariant equilirium): c iterate independent variable to refine c the conditions of the invariant point: 20 v(ivi) = v(ivi) - dv v(ivd) = v(ivd) 80 dv = dv/2.d0 v(ivi) = v(ivi) + dv if (dabs(dv).gt.delt(ivi)) goto 90 call warn (73,delt(ivi),ivi,'SFOL2 ') iflag = 0 goto 9000 90 call univeq (ivd,ier) if ((ier.eq.2).or.(ier.eq.1)) goto 9000 call pchk (igo) goto (50),igo goto (20),iflag c check if in range (a solution may c go out of range because of a c singularity in an equation of c state). if (((v(ivi).gt.vmax(ivi)).or.(v(ivi).lt.vmin(ivi))).or. * ((v(ivd).gt.vmax(ivd)).or.(v(ivd).lt.vmin(ivd)))) then irchk(ird) = 1 goto 9000 end if call assptx c these conditions may appear c redundant in formatted output c because of roundoff, in most cases c it may be sufficient to record c the values of v1 and v2 as simple c variables rather than in the out- c put array ptx, in which case the c apparent redundancy will be removed. goto 80 c assign the invariant assemblage: 50 call assip (ip) c save the conditions: call assptx c output the traversed equilibrium: 9000 ier = 0 call outrxn (ip,ier) c save conditions of the endpoint call svrend (ird,irend,jer) 999 end subroutine sfol2 (iovi,iovd,ip,irend) c---------------------------------------------------------------------- c sfol2 generates and computes the loci of stable equilibria c emanating from an ip indexed by the array ipid and the c arguement (i). since one of these equilibria is calculated c by sfol1 or by previous calls to sfol2 only c+1 (at most) c curves are calibrated. c c the phases identified by idv define a plane in g-x space, idphi c is the phase which is tangent to this plane, and lphi is the c ,or one of the, phases which is present at the ip but absent c along the univariant point. c---------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k2=1200,k3=1200,k8=9,l2=5,k5=12,nb=20) parameter (k6=7) character*8 names, ipnms, irnms, text(256)*1, alpha(90)*1 integer jdv(k6) common/ cst8 /names(k1)/ cst6 /icomp,istct,iphct,icp * / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar * / cst23 /a(10,10),b(11),ipvt(10),idv(10), * iophi,idphi,iiphi,iflg1 * / cst5 /v(l2),tr,pr,r,ps/ cst32 /ptx(450),logx,ipt2 * / cst29 /vip(l2,k3),ipid(k3,k8),ipct * / cst35 /ir1(k2),iirct,irpct,irv(k2)/ csta1 /irnms(k2,11) * / cst31 /vn(k2,11),ir(k2,11),irct,ird/ cst105 /ibug(k2) common/ cst89 /i2c(2,21),i3c(3,35),i4c(4,35),i5c(5,21),i6c(6,7) * / cst64 /i7c(7,8),i8c(8,9)/ cst25 /vnu(11),idr(11),ivct * / cst34 /ip1(k3),iipct,loop/ csta3 /ipnms(k3,k8) * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5 * / cst103 /isec,icopt,ifull,imsg,io3p * / cst300 /cblk(nb,k5),atw(k1),vol(k1),iblk,jcby,jcth c----------------------------------------------------------------------- jchk = 0 icp2 = icp + 2 ivi =iovi ivd=iovd ivarip(ip) = 1 goto (150), io3 goto (150), io3p call iptext (text,iend,ip,icp2) write (n3,1040) write (n3,1020) ip1(ip), (text(i),i = 1, iend) write (n3,1050) 150 do 10 j = 2, icp2 goto (20,30,40,50,51,52),icp c septary systems do 1 i = 1, icp 1 idv(i) = ipid(ip,i8c(i,j)) idphi = ipid(ip,i8c(8,j)) lphi = ipid(ip,icp2+1-j) goto 55 c hexary systems 52 do 2 i = 1, icp 2 idv(i) = ipid(ip,i7c(i,j)) idphi = ipid(ip,i7c(7,j)) lphi = ipid(ip,icp2+1-j) goto 55 c quinary systems 51 idv(1) =ipid(ip,i6c(1,j)) idv(2) =ipid(ip,i6c(2,j)) idv(3) =ipid(ip,i6c(3,j)) idv(4) =ipid(ip,i6c(4,j)) idv(5) =ipid(ip,i6c(5,j)) idphi =ipid(ip,i6c(6,j)) lphi =ipid(ip,icp2+1-j) goto 55 c quaternary systems 50 idv(1) =ipid(ip,i5c(1,j)) idv(2) =ipid(ip,i5c(2,j)) idv(3) =ipid(ip,i5c(3,j)) idv(4) =ipid(ip,i5c(4,j)) idphi =ipid(ip,i5c(5,j)) lphi =ipid(ip,icp2+1-j) goto 55 c ternary systems 40 idv(1) =ipid(ip,i4c(1,j)) idv(2) =ipid(ip,i4c(2,j)) idv(3) =ipid(ip,i4c(3,j)) idphi =ipid(ip,i4c(4,j)) lphi =ipid(ip,icp2+1-j) goto 55 c binary systems 30 idv(1) =ipid(ip,i3c(1,j)) idv(2) =ipid(ip,i3c(2,j)) idphi =ipid(ip,i3c(3,j)) lphi =ipid(ip,icp2+1-j) goto 55 c unary systems 20 idv(1) =ipid(ip,i2c(1,j)) idphi =ipid(ip,i2c(2,j)) lphi =ipid(ip,icp2+1-j) c balance the equilibrium, ier = 1 c if the plane idv is degenerate 55 itic= 0 c set invariant point conditions v(ivd) = vip(ivd,ip) v(ivi) = vip(ivi,ip) c check if assemblage sees bulk comp if (iblk.ne.0) then do 56 i = 1, icp 56 jdv(i) = idv(i) call bulkun (jdv,lphi,jbfg) if (jbfg.eq.0.and.jcby.eq.1) goto 10 end if 110 call balanc (b,idv,idphi,ier) goto (80),ier goto (130),ivar ivarip(ip) =2 goto (10),isudo c call pivots for lchk. 130 call pivots (ier) if (ier.eq.0) goto 90 c in some cases singularity of the c transpose of the concentration matrix c may not be detected by balanc, but, c will be detected for the concentration c matrix by pivots, in which case the c columns are permuted and the stoichiometry c of the reaction is redetermined. 80 if (itic.eq.icp-1) goto 100 icter =idphi c if a degenerate facet has been c generated by permutation permute c idphi with the vertices idv until c balanc succeeds. for a c component c system no more than c-1 permutations c will be necessary. idphi =idv(icp-itic) idv(icp-itic) =icter itic=itic+1 goto 110 100 goto 10 90 call assir c TESTING c test endpoint. call svrend (ird,irend,ier) goto (10),ier c calculate the loci of the equilibia, ipt2 = 0 call delrxn call assptx c determine direction of traverse and c dependent and independent variables call wway (div,lphi,ivi,ivd,ikwk,ier) c begin traverse if (ier.ne.0) then c call rxntext to save the equilibrium c in the list in case it is not found c again: call rxntxt (iend,jend,text,alpha) c set the bug flag: ibug(iirct) = 1 goto 10 end if icter = 0 60 call sfol1 (ivd,ivi,ier,div,lphi,ikwk,irend) jchk = jchk + 1 goto (70,70),ier ivd = iovd ivi = iovi goto 10 70 call switch (div,ivi,ivd,jer) goto (75),jer icter = icter + 1 if (icter.lt.4) goto 60 75 call warn (20,v(ivi),ivi,'SFOL2 ') call outrxn (ip,1) ivi =iovi ivd=iovd 10 continue if (jchk.eq.0) call warn (74,v(ivi),ivi,'SFOL2 ') if (io3p.eq.0) write (n3,1040) 1020 format (' equilibria about invariant point (',i3,'):',//, * 3x,200a1) 1040 format (' ------') 1050 format (/,' are listed below:',/) end subroutine sreset (jok,iste,vt,v,vti,vi) c---------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) c if the assemblages conditions c of stability are already known (jok=5) c or are known to extend beyond the c current condition (jok>iste) quit: if (jok.gt.iste) goto 99 c if the known conditions of stability c have been extened, i.e., either: c jok=iste or v>vt if (jok.ne.iste) then jok=iste else if (iste.lt.3) then if (vt.gt.v) goto 99 else if (vt.lt.v) goto 99 end if vt = v vti = vi 99 end subroutine svrend (ird,irend,ier) c---------------------------------------------------------------------- c save end points of univariant equilibria and c check if the intersection has already been recorded (ier =1). c two cases: c 1) reaction occurs on an edge, if the intersection has already c occurred then ignore future occurrences. c 2) reaction occurs at an IP only ignore occurences being traced in c the same direction if the reaction is degenerate. c---------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (l2=5,k2=1200,h5=7,h6=500) dimension rends(2,8,k2),irends(k2) common/ cst24 /ipot,jv(l2),iv(l2)/ cst32 /ptx(450),logx,ipt2 * / cst35 /ir1(k2),iirct,irpct,irv(k2) * / cst103 /isec,icopt,ifull,imsg,io3p * / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2 * / cst5 /v(l2),tr,pr,r,ps/ cst9 /vmax(l2),vmin(l2),ddv(l2) save rends, irends c----------------------------------------------------------------------- ier = 0 c if isec = 1, don't do any checks goto (99) isec dx = ddv(iv(1)) dy = ddv(iv(2)) x = v(iv(1)) y = v(iv(2)) if (x.ne.vmin(iv(1)).and.x.ne.vmax(iv(1)) .and. * y.ne.vmin(iv(2)).and.y.ne.vmax(iv(2))) then c the endpoint is not on an edge: if (isec.eq.2) then c if isec = 2, don't check unless c curve intersects at an edge. goto 99 else if (isec.eq.3) then c if isec = 3, don't do a check for c degenerate interior endpoints: if (irv(ird).lt.icp1) goto 99 end if else c the endpoint is on an edge: c if isec = 4, don't do a check for c degenerate equilibria. if (isec.eq.4.and.irv(ird).lt.icp1) goto 99 end if c isec = 5, do check on all ends if (ird.le.irend) then c reaction has been found already c look to see if its at the same c endpoint (tolerance +/- dv). jct = irends(ird) do 1 j = 1, jct xx = rends(1,j,ird) yy = rends(2,j,ird) if ((xx.gt.x-dx).and.(xx.lt.x+dx).and. * (yy.gt.y-dy).and.(yy.lt.y+dy)) then c condition match has occurred ier = 1 goto 99 end if 1 continue c no condition match has occurred jct = jct + 1 if (jct.gt.8) then call warn (999,x,ier,'SVREND') jct = 8 end if else c reaction has not been recorded jct = 1 irend = irend + 1 end if if (ird.gt.k2) call error (206,x,k2,'SVREND') irends(ird) = jct rends(1,jct,ird) = x rends(2,jct,ird) = y 99 end subroutine ufluid (fo2) c---------------------------------------------------------------------- c subroutine ufluid computes the potential of the components c of a saturated fluid phase. if the mole fraction of a component is les c less than 1.d-38 the chemical potential is set to -9.9d09. c ufluid may call one of three molecular fluid equations of state, or c alternatively users may supply their own routines, however, c the routines currently in use return the log of a components fugacity c which is then added to the reference state potential computed by the c function gphase. c---------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) parameter (k1=1200,k4=16,h5=7) dimension xf(2) common/ cst1 /thermo(k4,k1),vf(2,k1),vs(h5,k1),uf(2),us(h5) * / cst5 /p,t,xco2,u1,u2,tr,pr,r,ps/ cst11 /f(2) * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn/ cst208 /ifct,idfl c----------------------------------------------------------------------- c compute the chemical potentials of c fluid components in fluid saturated c systems. call cfluid (fo2,fs2) if (idfl.ne.0) then uf(idfl) = gphase(idfl) + r * t * f(idfl) else xf(1) = 1.d0 - xco2 xf(2) = xco2 do 10 i = 1, 2 if (iff(i).eq.0) goto 10 if (xf(i).gt.1.d-38) goto 70 uf(i) = -9.9d09 goto 10 70 uf(i) = gphase(i) + r * t * f(i) 10 continue end if end subroutine uproj c---------------------------------------------------------------------- c subroutine uproj computes the potential of mobile c or saturated components. the potentials of c saturated non-volatile components are projected through c saturated volatile components. c c uproj is always called before the routines gphase or gproj, in c therefore in order to increase program efficiency the values c log(t),dsqrt(t),(p-pr)*dexp(-t/300.d0), and dexp(-p/35000.d0) c are maintained in common block 'cst44'. c---------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) parameter (k1=1200,k4=16,k5=12,h5=7,h6=500) double precision uss(h6) common/ cst1 /thermo(k4,k1),vf(2,k1),vs(h5,k1),uf(2),us(h5) * / cst5 /p,t,xco2,u1,u2,tr,pr,r,ps * / cst6 /icomp,istct,iphct,icp/ cst20 /comps(k1,k5) * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn * / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2 * / cst44 /dlogt,dsqrtt,dexpt,dexpp c----------------------------------------------------------------------- c compute isobaric isothermal constants. if (t.lt.0.0.or.tr.eq.0.0) then dlogt = 0.0 dsqrtt = 0.0 dexpt = 0.0 dexpp = 0.0 else dexpt = (p-pr)*dexp(t/(-300.d0)) dexpp = dexp(p/(-35000.d0)) dlogt = dlog(t) dsqrtt = dsqrt(t) end if fo2 = 0.0 goto (50),ifyn c compute the chemical potentials of c fluid components in fluid saturated c systems. call ufluid (fo2) c determine stable saturated composants c and the corresponding chemical potentials c (these are projected through the fluid c components). 50 goto (99),isyn do 60 i = 1, isat ict = isct(i) ll = icp+i do 20 j = 1, ict k = ids(i,j) uss(j) = gphase(k)-vf(1,k)*uf(1)-vf(2,k)*uf(2) c if multiple component saturation constraints c project by saturation hierarchy: goto (20),i i1 = i-1 do 40 l = 1, i1 40 uss(j) = uss(j)-vs(l,k)*us(l) 20 uss(j) = uss(j)/comps(k,ll) c if O2, check if fo2 has been c determined by a fluid phase routine, c if so do the projection: if (io2.eq.i) then do 100 j = 1, ict 100 uss(j) = uss(j) + r * t * fo2 end if c now find stable "composant": u = uss(1) id = 1 goto (80),ict do 30 j = 2, ict if (uss(j).gt.u) goto 30 id = j u = uss(j) 30 continue c save the id of the stable composant. 80 idss(i) = ids(i,id) c and its chemical potential. 60 us(i) = u 99 end subroutine wway (div,lphi,ivi,ivd,ikwk,ier) c----------------------------------------------------------------------- c wway determines from which side of an invariant point a specific c univariant equilibrium emanates. the appropriate sign for the c increment of the independent variable (dv) is also determined. c where necessary wway also refines dv and/or may switch the c independent and independent intensive variables (ivi,ivd). vo1 c and vo2 are the values of the intensive variables iv1 and iv2 at c the invariant point. lphi is the id of the absent phase. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k2=1200,l2=5) character*8 names,irnms dimension vo(l2),vvi(2),vvd(2),ddv(2) common/ cst25 /vnu(11),idr(11),ivct/ cst8 /names(k1) * / cst35 /ir1(k2),iirct,irpct,irv(k2)/ csta1 /irnms(k2,11) * / cst31 /vn(k2,11),ir(k2,11),irct,ird * / cst6 /icomp,istct,iphct,icp * / cst23 /a(10,10),b(11),ipvt(10),idv(10), * iophi,idphi,iiphi,iflg1 * / cst5 /v(l2),tr,pr,r,ps/ cst9 /vmax(l2),vmin(l2),dv(l2) * / cst24 /ipot,jv(l2),iv1,iv2,iv3,iv4,iv5/ cst7 /iflag * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 c----------------------------------------------------------------------- iophi =idphi ikwk= 0 icter = 0 ier = 0 vo(iv1) = v(iv1) vo(iv2) = v(iv2) c first try iv1 as the iv ivi =iv2 ivd=iv1 ifail= 0 5 call univeq (ivd,ier) if (ier.eq.0) then c is the equilibrium condition c for this reaction consistent c with the invariant point? c if ( dabs( (v(ivd)-vo(ivd))/dv(ivd) ).gt.0.1d0) goto 9700 if ( dabs( (v(ivd)-vo(ivd))/dv(ivd) ).gt.0.1d0) then vo(iv1) = v(iv1) vo(iv2) = v(iv2) call warn (18,v(ivd),ird,'WWAY') end if else goto (9800),ifail ivi =iv1 ivd=iv2 v(iv1) = vo(iv1) v(iv2) = vo(iv2) ifail=ifail+1 goto 5 end if c test that the ip isn't on an edge if ((v(iv1).eq.vmin(iv1)).or.(v(iv1).eq.vmax(iv1)).or. * (v(iv2).eq.vmin(iv2)).or.(v(iv2).eq.vmax(iv2))) goto 9600 c order = 1.d0 c set iv increment 60 div = dv(ivi)/order ifail= 0 c why a do loop? do 10 i = 1, 2 c because of numerical slop c an equilibrium may appear stable c on its metastable extension close to c an ip, because of this wway always c tests for stability on both extensions c if wway cannot solve for equilibrium c conditions on both sides of an ip c the increment is decreased and if this c fails the independent and dependent c variables are switched. v(ivi) = vo(ivi) + div c check if the increment is too large if (v(ivi).gt.vmax(ivi)) v(ivi) = vmax(ivi) if (v(ivi).lt.vmin(ivi)) v(ivi) = vmin(ivi) c set the dependent variable to ip: v(ivd) = vo(ivd) c solve for the equilibrium call univeq (ivd,ier) if ((v(ivd).gt.vmax(ivd)).or.(v(ivd).lt.vmin(ivd))) goto 20 c on error switch dependent and c independent variables goto (20,20),ier vvi(i) = v(ivi) vvd(i) = v(ivd) ddv(i) = div c determine if the equilibrium is c metastable with respect to lphi c (lchk=1): if (lchk(lphi).eq.0) then ifail = ifail + 1 isign = i end if c end of looop 10 div = -div c at this stage there are three c possibilities indicated by ifail goto (50,20),ifail c ifail= 0 the equilibrium was stable c on both sides, reduce the increment c by an order of magnitude: order = order / 2. if (order.lt.0.001d0) goto 9500 ivi =iv1 ivd=iv2 goto 60 c---------------------------------------------------------------------- c ifail=2 the equilibrium was metastable c on c both sides if this has occurred twice c (icter =2) write error message, c otherwise switch ivs: c20 icter =icter+1 c if (icter.gt.3) goto 9900 c switch ivs c iv =ivi c ivi =ivd c ivd=iv c and re-enter the loop c---------------------------------------------------------------------- c test segment to replace above c ifail = 2, expand the search c increment 20 order = order / 2. if (order.lt.1.d-6) goto 9500 ivi = iv1 ivd = iv2 goto 60 c ifail=1, c the stable extension was identified c determine the sign of the increment c and its magnitude: 50 if (isign.eq.1) then div = dv(ivi) del= vmax(ivi)-vo(ivi) if (del.lt.div) div = del/3.d0 else div = -dv(ivi) del= vo(ivi)-vmin(ivi) if (del.lt.-div) div = -del/3.d0 end if c reset variable values and save v(ivi) = vvi(isign) v(ivd) = vvd(isign) odiv = ddv(isign) icter = 0 c now check the stability with respect c to all phases order = 1.d0 70 call gall (istct,iphct) call schk (lphi) c stable if iflag= 0 from schk c now check if dependent variable is in c range. if (iflag.eq.0) goto 9999 c metastable try refining the increment order = order * 5.d0 div = div / order v(ivi) = vo(ivi) + div v(ivd) = vo(ivd) 100 call univeq (ivd,ier) if (ier.eq.0) goto 90 c error from univeq, more than once quit goto (9400), icter div =odiv icter = icter + 1 goto 100 c no error from univeq 90 if (order.gt.10.d6) goto 9500 goto 70 c all systems go 9999 call assptx if (dabs(div).lt.dv(ivi)) then div = dv(ivi)*div/dabs(div) end if if (v(ivi).eq.vmax(ivi).or.v(ivi).eq.vmin(ivi)) ikwk= 1 return 9400 call warn (87,dv(1),ivi,'WWAY') ier = 1 return 9500 call warn (24,dv(ivi)/order,ivi,'WWAY') ier = 1 return 9600 call warn (63,dv(1),ivi,'WWAY') ier = 1 return c9700 write (n3,1040) v(iv2),v(iv1) c goto 9810 9800 call warn (58,dv(1),ier,'WWAY') c goto 9810 c9900 call warn (108,dv(1),ier,'WWAY') c9810 if (ivct.gt.4) goto 9820 write (n3,1050) ir1(ird),(vnu(zl),names(idr(zl)),zl= 1, ivct) goto 9830 9820 write (n3,1050) ir1(ird),(vnu(zl),names(idr(zl)),zl= 1,4) write (n3,1060) (vnu(zl),names(idr(zl)),zl=5,ivct) 9830 write (n3,*) ier = 1 1050 format (1x,'(',i3,')',4(1x,g9.3,1x,a8)) 1060 format (6x,4(1x,g9.3,1x,a8),/,6x,4(1x,g9.3,1x,a8)) c1040 format (/,' **warning ver059** wway, the equilibrium conditions', c * '(',2(1x,g10.4),')',/,' for the following reaction are' c * ,' inconsistent with the invariant conditions',/) end subroutine fultxt (nchar,rtxt) c----------------------------------------------------------------------- c subprogram to write text for complete reaction reactions c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k1=1200,k2=1200,k3=1200,h5=7,h6=500) character*8 names, text(256,2)*1, rtxt(256)*1, * ppart(8)*30, rxnstr*132, char8, mpart(8)*30, * cname*5, exten(2)*8, ptext*256, mtext*256 integer ichar(2) common/ csta4 /cname(12)/ cst21 /du(2),dv(2),jds(h5),ifr * / cst88 /vu(2,k1),jprct,jmct,jmyn * / cst201 /vuf(2),vus(h5),iffr,isr * / cst25 /vnu(11),idr(11),ivct/ cst8 /names(k1) * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn * / cst40 /ids(h5,h6),isct(h5),icp1,isat,io2 * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst61 /ikp(k1),ivarrx(k2),ivarip(k3),isudo,ivar * / cst31 /vn(k2,11),ir(k2,11),irct,ird * / cst80 /jrchk(k2),irchk(k2)/ cst104 /rxnstr(k2) * / cst103 /isec,icopt,ifull,imsg,io3p save exten c this is a bullshit trick and c will cause errors if someone c uses a function other than G. data exten/'-V(j/b) ','S(j/k) '/ ip = 0 im = 0 do 1 i = 1, 8 mpart(i) = ' ' 1 ppart(i) = ' ' c load id's of phases with + or - coeffs do 10 i = 1, ivct c id = idr(i) if (vnu(i).lt.0.) then im = im + 1 call wrpart (-vnu(i),ikp(id),names(id),mpart(im)) else ip = ip + 1 call wrpart (vnu(i),ikp(id),names(id),ppart(ip)) end if 10 continue c composant stoichiometry: if (isyn.eq.0) then do 20 i = 1, isat if (vus(i).gt.0.) then im = im + 1 call wrpart (vus(i),0,names(jds(i)),mpart(im)) else if (vus(i).lt.0) then ip = ip + 1 call wrpart (-vus(i),0,names(jds(i)),ppart(ip)) end if 20 continue end if c fluid stoichiometry: if (ifyn.eq.0) then do 30 i = 1, 2 if (iff(i).ne.0) then if (vuf(i).gt.0.) then im = im + 1 call wrpart (vuf(i),0,names(i),mpart(im)) else if (vuf(i).lt.0.) then ip = ip + 1 call wrpart (-vuf(i),0,names(i),ppart(ip)) end if end if 30 continue end if c mobile components: if (jmyn.eq.0) then do 40 i = 1, jmct char8 = cname(jprct+i) if (du(i).gt.0.) then im = im + 1 call wrpart (du(i),0,char8,mpart(im)) else if (du(i).lt.0.) then ip = ip + 1 call wrpart (-du(i),0,char8,ppart(ip)) end if 40 continue end if if (ifull.eq.1.or.ifull.eq.3) goto 55 do 50 i = 1, 2 if (dv(i).gt.0.) then im = im + 1 call wrpart (dv(i),0,exten(i),mpart(im)) else if (dv(i).lt.0.) then ip = ip + 1 call wrpart (-dv(i),0,exten(i),ppart(ip)) end if 50 continue 55 if (ip.gt.8.or.im.gt.8) call error (997,du(1),ip,'FULTXT') write (ptext,1000) (ppart(i), i = 1, ip) write (mtext,1000) (mpart(i), i = 1, im) ichar(1) = im * 31 read (mtext,1010) (text(i,1), i = 1, ichar(1)) ichar(2) = ip * 31 read (ptext,1010) (text(i,2), i = 1, ichar(2)) c filter out blanks: do 60 i = 1, 2 nchar = 1 do 70 j = 2, ichar(i) if (text(j-1,i).eq.' '.and.text(j,i).eq.' ') goto 70 nchar = nchar + 1 text(nchar,i) = text(j,i) 70 continue c filter out blanks before parentheses: ichar(i) = nchar nchar = 0 do 80 j = 1, ichar(i)-1 if ((text(j,i).eq.' '.and.text(j+1,i).eq.')').or. * (text(j,i).eq.' '.and.text(j+1,i).eq.'(')) goto 80 nchar = nchar + 1 text(nchar,i) = text(j,i) 80 continue ichar(i) = nchar 60 continue c concatenate parts into rxnstr: if (ichar(1)+ichar(2).gt.253) call error (208,du(1),ip,'FULTXT') nchar = 0 do 90 i = 1, 2 do 100 j = 1, ichar(i) nchar = nchar + 1 100 rtxt(nchar) = text(j,i) if (i.eq.2) goto 90 nchar = nchar + 3 rtxt(nchar-2) = ' ' rtxt(nchar-1) = '=' rtxt(nchar) = ' ' 90 continue c this is completely wrong c if old lists are being used, c i.e. iirct ne irct. if (jrchk(ird).eq.0) then jrchk(ird) = 1 if (nchar.gt.132) then write (rxnstr(ird),1010) (rtxt(i),i=1,132) else write (rxnstr(ird),1010) (rtxt(i),i=1,nchar) end if end if 1000 format (8(a30,1x)) 1010 format (256a1) end subroutine wrpart (vnu,ikp,name,part) c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (h9=18) character*8 name, fname*10, part*30 common/ csta7 /fname(h9) * / cst103 /isec,icopt,ifull,imsg,io3p if (ikp.eq.0) then if (ifull.gt.2) then write (part,1000) vnu,name else write (part,2000) name end if else if (ifull.gt.2) then write (part,1010) vnu,fname(ikp),name else write (part,2010) fname(ikp),name end if end if 1000 format (g9.3,1x,a8) 1010 format (g9.3,1x,a10,'(',a8,')') 2000 format (a8) 2010 format (a10,'(',a8,')') end subroutine bulkun (jdv,lphi,jbfg) c---------------------------------------------------------------------- c check if any permutation of a c+1 phase assemblage sees a c specified bulk composition, if no jbfg = 1. c---------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (k6=7) integer jdv(k6),kdv(k6+1) common/ cst6 /icomp,istct,iphct,icp/ cst64 /i7c(7,8),i8c(8,9) * / cst89 /i2c(2,21),i3c(3,35),i4c(4,35),i5c(5,21),i6c(6,7) * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst23 /a(10,10),b(11),ipvt(10),idv(7),id8,id9,id10, * iophi,idphi,iiphi,iflg1 c----------------------------------------------------------------------- do 15 i = 1, icp 15 kdv(i) = jdv(i) kdv(icp+1) = lphi do 10 j = 1, icp+1 goto (20,30,40,50,51,52),icp c septary systems do 1 i = 1, icp 1 idv(i) = kdv(i8c(i,j)) goto 55 c hexary systems 52 do 2 i = 1, icp 2 idv(i) = kdv(i7c(i,j)) goto 55 c quinary systems 51 idv(1) = kdv(i6c(1,j)) idv(2) = kdv(i6c(2,j)) idv(3) = kdv(i6c(3,j)) idv(4) = kdv(i6c(4,j)) idv(5) = kdv(i6c(5,j)) goto 55 c quaternary systems 50 idv(1) = kdv(i5c(1,j)) idv(2) = kdv(i5c(2,j)) idv(3) = kdv(i5c(3,j)) idv(4) = kdv(i5c(4,j)) goto 55 c ternary systems 40 idv(1) = kdv(i4c(1,j)) idv(2) = kdv(i4c(2,j)) idv(3) = kdv(i4c(3,j)) goto 55 c binary systems 30 idv(1) = kdv(i3c(1,j)) idv(2) = kdv(i3c(2,j)) goto 55 c unary systems 20 idv(1) = kdv(i2c(1,j)) c check the assemblage 55 call bulkck (idv,jbfg,ier) c ier = 1, degenerate if (ier.eq.1) goto 10 c jbfg ne 0, bounds a comp. if (jbfg.ne.0) goto 99 10 continue 99 end subroutine bulkck (jdv,jbfg,ier) c---------------------------------------------------------------------- c bulkck determines phases identified by jdv bound c any of the specified bulk compositions. c---------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) character cname*5, names*8 parameter (k1=1200,nb=20,h5=7,k4=16,k0=25,k5=12,l2=5) integer jdv(7) double precision fbulk(k5) common/ cst6 /icomp,istct,iphct,icp/ cst20 /comps(k1,k5) * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn * / cst5 /v(l2),tr,pr,r,ps * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst43 /therm(k4),comp(k0),atwt(k0),idh2o,idco2,ikind * / cst300 /cblk(nb,k5),atw(k1),vol(k1),iblk,jcby,jcth * / cst103 /isec,icopt,ifull,imsg,io3p/ cst8 /names(k1) * / cst301 /a(k5,k5),b(k5),ipvt(k5)/ csta4 /cname(k5) c---------------------------------------------------------------------- jbfg = 0 c load the transpose of the c concentration matrix of the pseudo- c invariant assemblage. do 10 i = 1, jcth if (i.le.icp) then k = jdv(i) else k = idss(i-icp) end if do 20 j = 1, jcth 20 a(j,i) = comps(k,j) 10 continue c factor the matrix call factr1 (jcth,ier) if (ier.eq.1) goto 99 ier = 1 c test for bounded compositions do 50 j = 1, iblk do 60 i = 1, jcth c load composition vector, the alpha c vector is returned in the same array 60 b(i) = cblk(j,i) c solve for the alpha vector call subst1 (jcth) c a negative coefficient indicates c the composition is not bounded do 70 i = 1, icp 70 if (b(i).lt.-1.d-5) goto 50 c the composition is bounded jbfg = 1 ier = 0 tot = 0. atot = 0. vtot = 0. do 80 i = 1, jcth if (i.le.icp) then k = jdv(i) if (b(i).lt.0.) b(i) = 0. else if (b(i).lt.0.) jbfg = 2 k = idss(i-icp) end if vtot = vtot + b(i) * vol(k) atot = atot + b(i) * atw(k) 80 tot = tot + b(i) c normalize molar proportions: rho = atot / vtot tot = tot / 100. atot = atot / 100. vtot = vtot / 100. do 90 i = 1, icomp 90 fbulk(i) = 0. write (n3,1000) j,(names(jdv(i)), i = 1, icp) write (n3,1010) write (n3,1020) c get true bulk do 100 i = 1, jcth if (i.le.icp) then k = jdv(i) else k = idss(i-icp) end if write (n3,1030) names(k), b(i)/tot, b(i)*atw(k)/atot, * b(i)*vol(k)/vtot do 110 l = 1, icomp 110 fbulk(l) = fbulk(l) + b(i) * comps(k,l) 100 continue if (jbfg.eq.2) then write (n3,1050) goto 50 end if gtot = 0. atot = 0. do 120 i = 1, icomp atot = atot + fbulk(i)*atwt(i) 120 gtot = gtot + fbulk(i) write (n3,1040) do 130 i = 1, icomp 130 write (n3,1030) cname(i), fbulk(i)/gtot*1.d2, * fbulk(i)*atwt(i)/atot*1.d2 write (n3,1060) rho write (n3,1070) v 50 continue 1000 format (/,9x,'---------',//,' Composition ',i3,' is bounded' * ,' by the assemblage:',/,2x,7(a8,1x)) 1010 format (' in the thermodynamic composition space.') 1020 format (/,12x,' mol %',6x,'wt %',5x,'vol %',/) 1030 format (1x,a8,3x,3(f6.2,4x)) 1040 format (/,12x,' mol %',6x,'wt %',/) 1050 format (/,' WARNING - negative mass indicates the bulk' * ,' composition is',/,' inconsistent with the' * ,' saturation constraints',/) 1060 format (/,' Reference state density (g/cm3): ',f6.3,/) 1070 format (/,' (',5(g12.6,1x),')',//,9x,'---------',/) 99 end subroutine subst1 (n) c----------------------------------------------------------------------- c subst uses the lu decomposition of the matrix 'a' contained c in the array 'a' to solve ax = b for x. subst is modified from the c the subroutine of the same name listed by conte and de boor c in 'elementary numerical analysis', mcgraw-hill, 1980. c factor uses scaled partial pivoting. c input a- an n by n array containing the non-zero elements of c the u and l decompositions of a, as output by factor. c n- the dimension of the matrix a. c ipvt- a vector indicating that row ipvt(k) was used to c eliminate the coefficient a(n,k). c b- the vector b. c output b- the solution vector x. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) parameter (k5=12) double precision x(k5) common / cst301 /a(k5,k5),b(k5),ipvt(k5) c solve ly = b for y: ip = ipvt(1) x(1) = b(ip) do 10 i = 2, n sum = 0.d0 im1 = i-1 do 20 j = 1, im1 20 sum = a(i,j)*x(j)+sum ip = ipvt(i) 10 x(i) = b(ip)-sum c solve ux = y for x: x(n) = x(n)/a(n,n) nm1 = n-1 do 30 ii = 1, nm1 i = n-ii ip1 = i+1 sum = 0. do 40 j = ip1, n 40 sum = a(i,j)*x(j)+sum x(i) = (x(i)-sum)/a(i,i) 30 b(i) = x(i) b(n) = x(n) end subroutine factr1 (n,ier) c----------------------------------------------------------------------- c factr1 is a subroutine which calculates the triangular c decompositions of the matrix 'a'. factor is modified from c the subroutine of the same name given by conte and de boor c in 'elementary numerical analysis', mcgraw-hill, 1980. c factor uses scaled partial pivoting. c c input a- an n by n array containing the elements of matrix a. c n- the dimension of the matrix a. c output a- an n by n array containing the upper, u, and lower, l, c triangular decompositions of input matrix a. c ipvt- a vector indicating that row ipvt(k) was used to c eliminate the a(n,k). c ier- a flag, zero if a is of rank = n, and 1 if a is of c lower rank. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) parameter (k5=12) common/ cst301 /a(k5,k5),d(k5),ipvt(k5) c----------------------------------------------------------------------- ier = 0 c initialize ipvt,d do 10 i = 1, n ipvt(i) = i rmax = 0. do 20 j = 1,n 20 rmax = dmax1(rmax,dabs(a(i,j))) c ax = b is singular if rmax = 0 if (rmax.eq.0.) goto 9000 10 d(i) = rmax c begin decomposition: nm1 = n-1 do 30 i = 1,nm1 ip1 = i+1 c determine pivot row (istr). rmax = dabs(a(i,i))/d(i) istr = i do 40 j = ip1,n tmax = dabs(a(j,i))/d(j) if (tmax.le.rmax) goto 40 rmax = tmax istr = j 40 continue if (rmax.eq.0.) goto 9000 c if istr gt i, make i the pivot row c by interchanging it with row istr. if (istr.le.i) goto 50 j = ipvt(istr) ipvt(istr) = ipvt(i) ipvt(i) = j temp = d(istr) d(istr) = d(i) d(i) = temp do 60 j = 1,n temp = a(istr,j) a(istr,j) = a(i,j) 60 a(i,j) = temp c eliminate x(k) from rows k+1,...,n. 50 do 70 j = ip1,n a(j,i) = a(j,i)/a(i,i) ratio = a(j,i) do 80 k = ip1,n 80 a(j,k) = a(j,k)-ratio*a(i,k) 70 continue 30 continue if (a(n,n).eq.0.) ier = 1 return c algoritmic singularity. 9000 ier = 1 end subroutine money c open (97, file = 'money', iostat = ierr, status = 'old') if (ierr.ne.0) then c system could not find the file write (*,1000) stop 1000 format (/,' There is no money file, make one and', * ' deposit at least a dime.',/) end if read (97,*) imoney rewind (97) if (imoney.le.0) then write (*,1010) stop 1010 format (/,' You are out of money, deposit at least a dime'/) else imoney = imoney - 1 write (97,*) imoney end if close (97) end subroutine subst (a,ipvt,n,b) c----------------------------------------------------------------------- c subst uses the lu decomposition of the matrix 'a' contained c in the array 'a' to solve ax = b for x. subst is modified from the c the subroutine of the same name listed by conte and de boor c in 'elementary numerical analysis', mcgraw-hill, 1980. c factor uses scaled partial pivoting. c input a- an n by n array containing the non-zero elements of c the u and l decompositions of a, as output by factor. c n- the dimension of the matrix a. c ipvt- a vector indicating that row ipvt(k) was used to c eliminate the coefficient a(n,k). c b- the vector b. c output b- the solution vector x. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) dimension a(10,10),ipvt(10),b(11),x(10) c solve ly = b for y: ip = ipvt(1) x(1) = b(ip) do 10 i = 2,n sum = 0.d0 im1 = i-1 do 20 j = 1, im1 20 sum = a(i,j)*x(j)+sum ip = ipvt(i) 10 x(i) = b(ip)-sum c solve ux = y for x: x(n) = x(n)/a(n,n) nm1 = n-1 do 30 ii = 1, nm1 i = n-ii ip1 = i+1 sum = 0.d0 do 40 j = ip1, n 40 sum = a(i,j)*x(j)+sum x(i) = (x(i)-sum)/a(i,i) 30 b(i) = x(i) b(n) = x(n) end subroutine testit c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) parameter (k1=1200,k6=7,l2=5) character*8 names dimension inc(7) common/ cst23 /a(10,10),b(11),ipvt(10),idv(10),iophi,idphi, * iiphi,iflg1 * / cst8 /names(k1)/ cst2 /g(k1) * / cst6 /icomp,istct,iphct,icp/ cst12 /cp(k6,k1) * / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / cst52 /hcp,h,id(8),ic(7)/ cst87 /delt(l2),dtol,utol,ptol c optional routine for testing c for metastable g-x planes. igo = 0 istop = 1 if (hcp.eq.icp) igo = 1 call abload (*9030) goto(2050),igo iq = 0 do 20 j = 1, icp do 30 i = 1, hcp 30 if (ic(i).eq.j) goto 20 iq = iq+1 inc(iq) = j 20 continue 2050 do 2020 j = istct,iphct goto (2060),igo do 60 i = 1, iq 60 if (cp(inc(i),j).ne.0.0) goto 2020 2060 dg = g(j) do 10 i = 1, hcp 10 dg = dg - cp(ic(i),j) * b(i) if (dg.gt.dtol) goto 2020 c check that a phase is not metastable c with respect to itself, this c could be remedied by changing the c value of dtol. do 25 k = 1, hcp 25 if (id(k).eq.j) goto 2020 c write (n3,*) '**error ver004** testit metastable plane' write (n3,*) (names(id(i)),i = 1,hcp),'vs. ',names(j) write (n3,*) dg istop = 0 2020 continue goto (99),istop call outchm 9030 call error (14,dg,i,'TESTIT') 99 end function sidc (kd) c----------------------------------------------------------------------- c given the parameters of a plane spanning the join of a c-1 c dimensional simplex, e.g. id1,...,idc-1, sidc determines if c a point 'kd' lies on the same side of the join as the vertex c idc. parameters for the simplex id1,...,idc are initialized c by slopc. c sidc is 0 if kd lies in the plane, lt 0 if kd lies on the opposite c side of the join from idc, and gt 0 if kd is on the same side as c idc. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) parameter (k1=1200,k6=7) common/ cst55 /s(6),sign,ipc(6),jcp,kcp/ cst3 /x(k1,k6) sidc = s(jcp) - x(kd,ipc(jcp)) do 10 i = 1, kcp 10 sidc = sidc + s(i) * x(kd,ipc(i)) sidc = sidc / sign end subroutine slopc (j) c----------------------------------------------------------------------- c given the join of a simplex of c-1 dimensions, e.g. id1,...,idc-1, c slopc finds the equation of the plane spanning the join. the c displacement of the vertex not on the join, idc, is determined c relative to this plane (i.e. the parameter 'sign'). c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n,z) parameter (k1=1200,k6=7,j8=4800) common/ cst3 /x(k1,k6)/ cst134 /iperm(6,7,7) * / cst55 /s(6),sign,ipc(6),jcp,kcp * / cst52 /hcp,h,id(8),ic(7) * / cst23 /a(10,10),b(11),ipvt(10),idv(10),iophi,idphi, * iiphi,iflg1/ cst46 /ipsf,idpsf(j8,k6) lgo = 1 idtest = idpsf(j,hcp) c load matrix and vectors 60 do 70 i = 1, jcp 70 ipc(i) = ic(iperm(i,lgo,hcp)) do 10 i = 1, jcp do 20 k = 1, kcp 20 a(i,k) = x(idpsf(j,i),ipc(k)) a(i,jcp) = 1.d0 10 b(i) = x(idpsf(j,i),ipc(jcp)) c solve for the plane: call factor (a,jcp,ipvt,ier) c if a is singular (ier = 1) then c switch the components and reload. goto (30),ier call subst (a,ipvt,jcp,b) c c save parameters and compute sign: do 40 i = 1, jcp 40 s(i) = b(i) c compute sign: sign = s(jcp) - x(idtest,ipc(jcp)) do 50 i = 1, kcp 50 sign = sign + s(i) * x(idtest,ipc(i)) if (sign.eq.0.) then write (*,*) 'huffa puffa hubba' goto 30 end if sign = sign / dabs(sign) return 30 lgo = lgo + 1 if (lgo.gt.hcp) call error (998,sign,lgo,'SLOPC') goto 60 end