implicit real (a-g,o-y),integer (h-m,z) character*8 cnames(3), names(10) integer ids(8,20), ida(3,20),idpsf(5) real vs(3), cp(9,10) common / cst23 /a(10,10),b(11),ipvt(11) common / cst52 /hcp,h,id(8),ic(7)/ cst3 /x(9,10) * / cst55 /s1,s2,s3,s4,s5,s6,sign,lgo open (20,file='stoich.dat') read (20,*) icp, ics read (20,1000) (cnames(i), i = 1, ics) ictot = icp + ics iphct = icp + ics + 1 icp1 = icp + 1 isct = 0 do i = 1, ictot ic(i) = i end do do i = 1, iphct read (20,1000) names(i) read (20,*) (cp(j,i), j = 1, ictot) end do do i = 1, iphct tot = 0.d0 do j = 1, ictot tot = tot + cp(j,i) end do do j = 1, iphct x(j,i) = cp(j,i) / tot end do end do c make h2o x(4,7) = 1.0 c make co2 x(5,8) = 1.0 c irct = 0 do i = 1, ictot do j = i+1, iphct irct = irct + 1 ida(1,irct) = i ida(2,irct) = j ict = 0 do k = 1, iphct if (k.ne.i.and.k.ne.j) then ict = ict + 1 ids(ict,irct) = k end if end do end do end do write (*,*) 'there are a total of ',irct,' two-phase absent ', * 'reactions.' c now try side testing: do k = 1, irct do i = 1, 4 idpsf(i) = ids(i,k) end do idpsf(5) = ida(1,k) call slop5 (idpsf) write (*,2000) (names(idpsf(i)),i = 1, 4) * ,names(ida(1,k)),sid5(ida(1,k)) write (*,2000) (names(idpsf(i)),i = 1, 4) * ,names(ida(2,k)),sid5(ida(2,k)) write (*,2000) (names(idpsf(i)),i = 1, 4) * ,'h2o',sid5(7) write (*,2000) (names(idpsf(i)),i = 1, 4) * ,'co2',sid5(8) end do 2000 format (4(a8,1x),'vs ',a8,' has sidc=',g12.3) do k = 1, irct do i = 1, icp do j = 1, icp a(j,i) = cp(j,ids(i,k)) end do end do do j = 1, icp b(j) = cp(j,ids(icp1,k)) end do b(icp1) = -1.0 call factor (a,icp,ipvt,ier) goto (10), ier call subst (a,b,ipvt,icp) do j = 1, ics vs(j) = 0.0 end do do i = 1, icp1 do j = icp1, ictot vs(j-icp) = vs(j-icp) + b(i) * cp(j,ids(i,k)) end do end do if (vs(1)*vs(2).lt.0.0) then write (*,1030) names(ida(1,k)), names(ida(2,k)) goto 20 end if xc = vs(2) / (vs(1) + vs(2)) isct = isct + 1 write (*,1010) names(ida(1,k)), names(ida(2,k)), * cnames(ics), xc goto 20 10 write (*,1030) names(ida(1,k)), names(ida(2,k)) 20 write (*,1040) (b(i),i = 1, icp1) end do write (*,*) 'there are a total of ',isct,' two-phase absent ', * 'singular reactions.' 1000 format (8(a8,1x)) 1010 format ('the (',a2,',',a2,') singular composition is x(',a3, * ') = ', g10.4) 1020 format ('the (',a2,',',a2,') singular point does not exist', * ' or is degenerate') 1030 format ('the (',a2,',',a2,') singular point does not exist.') 1040 format (7(g9.3,1x)) end subroutine factor (a,n,ipvt,ier) implicit real (a-g,o-y),integer (h-m,z) dimension a(10,10),d(11),ipvt(11) ier=0 c initialize ipvt,d do 10 i=1,n ipvt(i)=i rmax=0.d0 do 20 j=1,n 20 rmax= max1(rmax, abs(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= abs(a(i,i))/d(i) istr=i do 40 j=ip1,n tmax= abs(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 c subroutine subst (a,b,ipvt,n) c implicit real (a-g,o-y),integer (h-m,z) c dimension a(10,10),ipvt(11),b(11),x(11) 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 slop5 (idpsf) implicit real (a-g,o-y),integer (h-n,z) integer idpsf(5) common/ cst3 /x(9,10)/ cst55 /s1,s2,s3,s4,s5,s6,sign,lgo * / cst23 /a(10,10),b(11),ipvt(11) * / cst52 /hcp,h,id(8),ic1,ic2,ic3,ic4,ic5,ic6,ic7 save ict data ict/4/ c load matrix and vectors lgo = 3 i5 = idpsf(5) 20 do 10 i = 1, 4 a(i,1) = x(1,idpsf(i)) a(i,2) = x(2,idpsf(i)) a(i,3) = x(lgo,idpsf(i)) a(i,4) = 1.d0 10 b(i) = x(5,idpsf(i)) c solve for the plane: call factor (a,ict,ipvt,ier) c if a is singular (ier = 1) then c switch the components and reload. goto (30),ier call subst (a,b,ipvt,ict) c save parameters: s1 = b(1) s2 = b(2) s3 = b(3) s4 = b(4) c compute sign: sign = s1*x(1,i5)+s2*x(2,i5)+s3*x(lgo,i5)+s4-x(5,i5) return 30 lgo = ic4 goto 20 end function sid5 (kd) implicit real (a-g,o-y),integer (h-n,z) common/ cst55 /s1,s2,s3,s4,s5,s6,sign,lgo/ cst3 /x(9,10) * / cst52 /hcp,h,id(8),ic1,ic2,ic3,ic4,ic5,ic6,ic7 c lgo is a flag which indicates c which components were used to c define the plane in the simplex. sid5 = (s1*x(1,kd)+s2*x(2,kd)+s3*x(lgo,kd)+s4-x(5,kd)) * /sign end