关于fortran 90里面的编程问题(急)

2025-02-27 06:18:48
推荐回答(1个)
回答1:

subroutine equil (aik,y,fe,elem,nele,ngas,ntot,ind)
c aik is elemental specie composition matrix nele x ntotal dimension
c y is no moles of each specie ntotal dimension
c fe is free energy of each specie ntotal dimension
c elem is no moles of each element nele dimension
c nele is no elements
c ngas is no gas species
c ntot is total no species
c ind is -7 if singular matrix and -1 if other than last solid disappeared
dimension aik(250),y(25),fe(25),elem(10),x(25),bmat(36),amat(1296)
c this common statement special to bkw code
common/subvar/amaxe,aminx,aminy,tx(10)
sst=0.
cpt=0.0
do 2 i=1,ntot
if(y(i).lt.aminy) y(i)=aminy
2 continue
35 nm1=ntot+nele+1
nn1=ntot-ngas
nm1sq=nm1*nm1
c zero amatrix
do 4 i=1,nm1sq
amat(i)=0.
4 continue
c form sums
sum=0.
do 5 i=1,ngas
sum=sum+y(i)
5 continue
bary=sum
rbary=(1.)/bary
c fill b matrix
i=1
do 6 j=1,ngas
bmat(i)=-(fe(i)+alog(rbary*y(i)))
i=i+1
6 continue
if(nn1.eq.0) go to 70
do 7 j=1,nn1
bmat(i)=-(fe(i))
i=i+1
7 continue
70 continue
do 8 j=1,nele
bmat(i)=elem(j)
i=i+1
8 continue
bmat(i)=0.
c fill in amatrix
i=1
do 9 j=1,ngas
amat(i)=(1.)/(y(j))
i=i+nm1+1
9 continue
i=ntot+1
j=1
l=1
12 do 10 k=1,nele
amat(i)=aik(j)
i=i+1
j=j+1
10 continue
if(l.gt.ngas)go to 11
amat(i)=1.0
11 l=l+1
if(l.gt.ntot)go to 13
i=i+1+ntot
go to 12
13 i=i+1
do 14 k=1,ngas
amat(i)=-rbary
i=i+1
14 continue
i=i+nn1+nele
amat(i)=-1.0
i=i+1
l=1
j=1
17 do 15 k=1,ntot
amat(i)=aik(j)
j=j+nele
i=i+1
15 continue
i=i+nele+1
l=l+1
if(l.gt.nele) go to 16
j=l
go to 17
16 call lss(nm1,1,nm1,amat,bmat,d,det,ind)
if(ind.lt.0) go to 40
c put answers in x
do 18 i=1,ntot
x(i)=bmat(i)
18 continue
c test for too small x
do 19 i=1,ngas
if(x(i).lt.aminx) x(i)=aminx
19 continue
c test to see if any solids disappeared
if(nn1.eq.0) go to 22
i=ngas+1
do 20 j=1,nn1
if(x(i).lt.0.) go to 50
i=i+1
20 continue
22 continue
c test to see if converged
80 sum=0.
do 21 i=1,ntot
sum=sum+abs(y(i)-x(i))
21 continue
if(sum.lt.amaxe) go to 60
c special go ahead with less strict error limit
cpt=cpt+1.0
if(cpt.gt.50.)write (3,991)sum
991 format(25h equil amaxe accepted was,1e18.11)
if(cpt.gt.50.)go to 60
c resolve with y now having last answers x
31 do 32 i=1,ntot
y(i)=x(i)
32 continue
go to 35
c have conve-ged
60 do 38 i=1,ntot
y(i)=x(i)
38 continue
j=ntot+1
if(sst.le.0.)go to 34
y(j)=0.
ntot=ntot+1
34 return
c solid has disappeared
50 sst=1.0
if(i.lt.ntot)go to 33
ntot=ntot-1
go to 35
c an error has occurred as other than last solid disappeared
c ind=-1 if error has occurred
33 ind=-1
return
c error has occurred matrix is singular
40 ind=-7
return
end