33
11 В1(N)=B1(N)+2.*Y(I)*Y(J)(GIJМ(I,J)+1D0-GN(N))**GN(N)*
*(QI(I)*QI(J)+1D0-QN(N))**QN(N)*((FI(I)*FI(J))**.5+
1D0-FN(N))**FN(N)*EIJM(I,J)**UN(N)*(KI(I)*KI(J)**1.5
7 CONTINUE
DO 13 N=8,53
13 C1(N)=AN(N)*(GM+1D0-GN(N))**GN(N)*(QM**2+1D0-QN(N))**
*QN(N)*(FM+1D0-FN(N))**FN(N)*UM**UN(N)
RETURN
END
SUBROUTINE PHASE1
IMPLICIT REAL*8(A-H,O-Z)
COMMON/Z/Z/RM/RM/T/T/P/P/AI1/AO,A1/AN/AN(53)
*/COEF1/B1(13),C1(53)/COEF2/B,C(53)/UN/UN(53)
CALL PCONT1(P,T)
IF(Z.EQ.0D0) GO TO 134
B=0D0
DO 1 N=1,13
1 B=B+AN(N)/T**UN(N)*B1(N)
DO 3 N=8,53
3 C(N)=C1(N)/T**UN(N)
PR=P/5.
RO=9D3*P/(RM*T*(1.1*PR+0.7D0))
CALL FUN1(RO)
Z=1D0+AO
134 RETURN
END
С Подпрограмма, реализующая итерационный процесс определения
С плотности из уравнения состояния (метод Ньютона)
SUBROUTINE FUN1(X)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/P/P/RM/RM/T/T/AI1/AO,A1
ITER=1
1 CONTINUE
CALL COMPL1(X)
Z=1.D0+AO
FX= 1 .D6*(P-(1.D-3*RM*T*Z*X))
F= 1 .D3*RM*(1.D0+A1)
DR=FX/F
X=X+DR
IF(ITER.GT.10) GO TO 4
ITER=ITER+1
IF(DABS(DR/X).GT.1.D-6) GO TO 1
4 CALL COMPL1(X)
RETURN
END
SUBROUTINE PCONT1(P,T)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/Z/Z
Z=-1D0
IF(T.LT.250D0.OR.T.GT.340D0)Z=0D0
IF(P.LE.0D0.OR.P.GT.12D0) Z=0D0
RETURN