编辑代码

program main         !主程序开始
REAL t,p,k,eps,ych4,yh2s,yh2,ys,ycs2,x,G,L          !定义实型变量,单精度
external f
common p,L

data t/1000/      !使用data命令赋值
data eps/1e-6/     !精度确定
data G1/41.17/,G2/17.45/,G3/19.38/,R/8.314e-03/     !使用DATA命令批量输入数值
G=G2+G1-G3           !定义数值,反应吉布斯自由能
L=exp(-G/R/t)         !定义数值,反应标准平衡常数
p=152000            !赋值
a=0.0               !赋值
b=1.0               !赋值
eps=1.0e-6          !精度确定
call half(a,b,eps,f,x)           !调用二分法子程序
ych4=(1-x)/(3-x)         !计算数值,平衡组成甲烷含量
yh2s=(x)/(3-x)           !计算数值,平衡组成硫化氢含量
ys=(2-3*x)/(3-x)         !计算数值,平衡组成硫含量
ycs2=(x)/(3-x)           !计算数值,平衡组成二硫化碳含量
yh2=x/(3-x)              !计算数值,平衡组成氢气含量
  Write(*,*)'T=:'        
write(*,*)t              !输出
 Write(*,*)'P=:' 
write(*,*)P              !输出
 Write(*,*)'x=' 
write(*,*)x              !输出结果
 Write(*,*)'ych4=:' 
write(*,*)ych4           !输出结果
Write(*,*)'ys=:' 
write(*,*)ys             !输出结果
Write(*,*)'ycs2=:' 
write(*,*)ycs2           !输出结果
 Write(*,*)'yh2s=:' 
write(*,*)yh2s           !输出结果
 Write(*,*)'yh2=:' 
write(*,*)yh2            !输出结果

end

function f(x)            !定义函数
real x,p,L               !定义实型变量,单精度
common p,L   
f=x*x*x*(3-x)*(1.013e5)/(1-x)/(2-3*x)/(2-3*x)/(2-3*x)/p-L            !函数表达式
return
end

SUBROUTINE HALF(A,B,EPS,F,X)              !二分法子程序开始
REAL A,B,EPS                              !定义实型变量,单精度
Y0=F(A)                                   !变量赋值
A1=A
B1=B
K=0
10 X=(A1+B1)*0.5
Y=F(X)
IF(Y*Y0.GT.0.0)A1=X                       !判断是否异号
IF(Y*Y0.LE.0.0)B1=X                       !判断是否异号
IF(ABS(B1-A1)/X.GT.EPS)GOTO 10            !迭代次数控制
K=K+1
IF(K.GT.50) THEN 
WRITE(*,*)'NO RESOLUTION'
GOTO 30
END IF 
X=(A1+B1)*0.5                             !输出返回x的值
30 RETURN
END                                       !二分法子程序结束