H=47.479*BMY+128.64D0
RETURN
3 Т=Т1
ТС=Т1-Т0
P=PRESS
IF(PRESS.LE.0D0.OR.PRESS.GT.12D0)Z=0D0
IF(T1.LT.250D0.OR.T1.GT.340D0)Z=0D0
IF(Z.EQ.0D0) GO TO 133
CALL B11BER(T,H,B11)
CALL BBER(T,B11,B,Z)
IF(Z.EQ.0D0) GO TO 133
CALL CBER(T,H,C,Z)
IF(Z.EQ.0D0) GO TO 133
CALL ITER2(P,T,B,C,Z)
133 RETURN
END
SUBROUTINE B11BER(T,H,B11)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/BBLOK/BR11H0(3),BR11H1(3),BR11H2(3),BR22(3),BR23(3),BR33(3)
T2=T*T
B11=BR11H0(1)+BR11H0(2)*T+BR11H0(3)*T2+
*(BR11H1(1)+BR11H1(2)*T+BR11H1(3)*T2)*H+
*(BR11H2(1)+BR11H2(2)*T+BR11H2(3)*T2)*H*H
END
SUBROUTINE BBER(T,B11,BEFF,Z)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/BBLOK/BR11H0(3),BR11H1(3),BR11H2(3),BR22(3),BR23(3),BR33(3)
COMMON/ZETA/Z12,Z13,Y12,Y13,Y123
COMMON/XBLOK/X1,X2,X3,X11,X12,X13,X22,X23,X33
T2=T*T
B22=BR22(1)+BR22(2)*T+BR22(3)*T2
B23=BR23(1)+BR23(2)*T+BR23(3)*T2
B33=BR33(1)+BR33(2)*T+BR33(3)*T2
BA13=B11*B33
IF(BA13.LT.0D0) THEN
Z=0D0
RETURN
ENDIF
ZZZ=Z12+(320D0-T)**2*1.875D-5
BEFF=X11*B11+X12*ZZZ*(B11+B22)+2.*X13*Z13*DSQRT(BA13)+
*X22*B22+2.*X23*B23+X33*B33
END
SUBROUTINE CBER(T,H,CEFF,Z)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/CBLOK/CR111H0(3),CR111H1(3),CR111H2(3),CR222(3),CR223(3),
*CR233(3),CR333(3)
COMMON/ZETA/Z12,Z13,Y12,Y13,Y123
COMMON/XBLOK/X1,X2,X3,X11,X12,X13,X22,X23,X33
T2=T*T
C111=CR111 H0(1)+CR111H0(2)*T+CR111H0(3)*T2+
*(CR111H1(1)+CR111H1(2)*T+CR111H1(3)*T2)*H+
*(CR111H2(1)+CR111H2(2)*T+CR111H2(3)*T2)H*H
C222=CR222(1)+CR222(2)*T+CR222(3)*T2
C223=CR223(1)+CR223(2)*T+CR223(3)*T2
C233=CR233(1)+CR233(2)*T+CR233(3)*T2
C333=CR333(1)+CR333(2)*T+CR333(3)*T2
CA112=C111*C111*C222
CA113=C111*C111*C333
CA122=C111*C222*C222
CA123=C111*C222*C333
CA133=C111<C333*C333
IF(CA112.LT.0D0.OR.CA113.LT.0D0.OR.CA122.LT.0DO0.OR.
*CA123.LT.0D0.OR.CA133.LT.0D0)THEN
Z=0D0
RETURN
ENDIF
D3REP=1D0/3D0
CEFF=X1*X11*C111+3D0*X11*X2*(CA112)**D3REP*(Y12+(T-270D0)*.0013D0)
*+3.*X11*X3*(CA113)**D3REP*Y13+
*3.*X1*X22*(CA122)**D3REP*(Y12+(T-270D0)*.0013D0)+
*6.*X1*X2*X3*(CA123)**D3REP*Y123+3.*X1*X33*(CA133)**D3REP*Y13+
*X22*X2*C222+3.*X22*X3*C223+3.*X2*X33*C233+X3*X33*C333
END
С Подпрограмма, реализующая схему Кардано для определения
С фактора сжимаемости из уравнения состояния
SUBROUTINE ITER2(P,T,Bm,Cm,Z)
IMPLICIT REAL*8(A-H,O-Z)
B1=1D3*P/2.7715/T
B0=Bl*Bm
C0=Bl**2*Cm
A1=1D0+B0
A0=1D0+1.5*(B0+C0)
A01=A0**2-A1**3
IF(A01.LE.0D0) THEN
Z=0D0
RETURN
ENDIF
A=A0-A01**0.5
A2=DABS(A)**(1D0/3D0)
IF(A-LT.0D0) A2=-A2
Z=(1D0+A2+A1/A2)/3.
END
BLOCK DATA BDGRG2
IMPLICIT REAL*8(A-H,O-Z)
COMMON/BBLOK/BR11H0(3),BR11H1(3),BR11H2(3),BR22(3),BR23(3),
*BR33(3)/CBLOK/CR111H0(3),CR111H1(3),CR111H2(3),CR222(3),
*CR223(3),CR233(3),CR333(3)
COMMON/ZETA/Z12,Z13,Y12,Y13,Y123
COMMON/MBLOK/GM2,GM3,FA,FB,TO,R
DATA BR11H0/-.425468D0,.2865D-2,-.462073D-5/,
* BR11H1/.877118D-3,-.556281D-5,.881514D-8/,
* BR11H2/-.824747D-6,.431436D-8,-.608319D-11/,
* BR22/-.1446D0,.74091D-3,-.91195D-6/,
* BR23/-.339693D0,.161176D-2,-.204429D-5/,
* BR33/-.86834D0,.40376D-2,-.51657D-5/
DATA CR111H0/-.302488D0,.195861D-2,-.316302D-5/,
* CR111 H1/.646422D-3,-.422876D-5,.688157D-8/,
* CR111H2/-.332805D-6,.22316D-8,-.367713D-11/,
* CR222/.78498D-2,-.39895D-4,.61187D-7/,
* CR223/.552066D-2,-.168609D-4,.157169D-7/,
* CR233/.358783D-2,.806674D-5,-.325798D-7/,
* CR333/.20513D-2,.34888D-4,-.83703D-7/
DATA Z12/.72D0/,Z13/-.865D0/,Y12/.92D0/,Y13/.92D0/,Y123/1.1D0/
DATA GM2/28.0135D0/,GM3/44.01D0/,
* FA/22.414097D0/,FB/22.710811D0/,
* TO/273.15D0/,R/.0831451D0/
END 46
C **********************************************************
С * *
С * Подпрограмма расчета коэффициента сжимаемости природного *
С * газа по уравнению состояния AGA8-92DC. *
C * *
C **********************************************************
SUBROUTINE AGA8DC(ICALC)
IMPLICIT REAL*8(A-H,O-Z)
REAL*8 KI,KIJ,KD
COMMON/RM/RM/Y1/Y(19)/NC1/NC/NI1/NI(19)/EFI/EI(19),KI(19),
*GI(19),QI(19),FI(19)
*/INTER1/EIJ(19,19),UIJ(19,19),KIJ(19,19),GIJ(19,19)
*/EFD/ED(19),KD(19),GD(19),QD(19),FD(19)/Z/Z
RM=8.31448D0
IF(ICALC.NE.l) GO TO 3
CALL COMPO1
IF(Z.EQ.0D0) GO TO 133
CALL PARIN1
DO 75 I=1,NC
EI(I)=ED(NI(I))
KI(I)=KD(NI(I))
GI(I)=GD(NI(I))
QI(I)=QD(NI(I))
FI(I)=FD(NI(I))
DO 123 J=1,NC
IF(I.GE.J) GO TO 123
EIJ(I,J)=EIJ(NI(I),NI(J))
UIJ(I,J)=UIJ(NI(I),NI(J))
KIJ(I,J)=KIJ(NI(I),NI(J))
GIJ(I,J)=GIJ(NI(I),NI(J))
123 CONTINUE
75 CONTINUE
CALL PARMI1
3 CALL PHASE1
133 RETURN
END
SUBROUTINE COMPO1
IMPLICIT REAL*8(A-h,O-Z)
DIMENSION ZNI(25),YI(25)
COMMON/YI/Y(19)/YI/YC(25)/NC1/NC/NT1/NI(19)/NPR/NPR
DATA ZNI/.9981D0,.992D0,.9834D0,.9682D0,.971D0,.9997D0,.9947D0,
*.99D0,.993D0,.994D0,985D0,.945D0,.953D0,1D0,.919D0,
*.936D0,.876D0,.892D0,3*1D0,1.0005D0,1.0006D0,.9996D0,.9993D0/
DO 100 I=1,25
100 YI(I)=YC(I)
YI(13)=YI(13)+YI(14)
YI(14)=0D0
IF(NPR.EQ.0D0) GO TO 5
YI(17)=YI(17)+YI(19)+YI(20)+YI(21)
YI(19)=0D0
YI(20)=0D0
YI(21)=0D0
SUM=0D0
DO 7 I=1,25
7 SUM=SUM+YI(I)/ZNI(I)
DO 9 I=1,25
9 YI(I)=YI(I)/ZNI(I)/SUM
5 YI(2)=YI(2)+YI(9)+YI(10)
YI(9)=0D0
YI(10)=O0D0
YI(3)=YI(3)+YI(11)
YI(11)=0D0
YI(15)=YI(15)+YI(16)
YI(16)=0D0
YI(17)=YI(17)+YI(18)
YI(18)=0D0
NC=0
ИС=0
YSUM=0D0
DO 11 1=1,25
IF((I.GE.9.AND.I.LE.11).OR.I.EQ.14.0R.I.EQ.16.0R.I.EQ.18)
*ИС=ИС+1
IF(YI(I).EQ.0D0) GO TO 11
NC=NC+1
NI(NC)=I-ИC
Y(NC)=YI(I)
YSUM=YSUM+Y(NC)
11 CONTINUE
CALL MOLDOl(YI)
DO 13 I=1,NC
13 Y(I)=Y(I)/YSUM
RETURN
END
SUBROUTINE MOLDO1(YI)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION YI(25)
COMMON/Z/Z
Z=-1D0
YS=0D0
DO 1 I=9,25
1 YS=YS+YI(I)
1F(YI(1).LT.0.65D0.OR.YI(2).GT.0.15D0.OR.YI(3).GT.0.035D0.OR.
*YI(4).GT.0.015D0.OR.YI(5).GT.0.015D0.OR.YS.GT.0.01D0) Z=0D0
IF(YI(6).GT.0.2D0.OR.YI(7).GT.0.15D0.OR.Y1(8).GT.5D-5) Z=O0D0
RETURN
END
SUBROUTINE PARIN1
IMPLICIT REAL*8(A-H,O-Z)
REAL*8 KIJ
COMMON/INTER1/EIJ(19,19),UIJ(19,19),KIJ(19,19),GIJ(19,19)
DO 1 I=1,19
DO 1 J=1,19
EIJ(I,J)=1D0
UIJ(I,J)=1D0
KIJ(I,J)=1D0
1 GIJ(I,J)=1D0
EIJ(1,6)=0.97164D0
UIJ(1,6)=0.886106D0
KIJ(1,6)=1.00363D0
EIJ(1,7)=0.960644D0
UIJ(1,7)=0.963827D0
KIJ(1,7)=0.995933D0
GIJ(1,7)=0.807653D0
EIJ(1,3)=0.99605D0
UIJ(1,3)=1.02396D0
EU(1,17)=1.17052D0
UIJ(1,17)=1.15639D0
KIJ(1,17)=1.02326D0
GIJ(1,17)=1.95731D0
EIJ(1,18)=0.990126D0
EIJ(1,5)=1.01953D0
EIJ(1,4)=0.995474D0
UIJ(1,4)=1.02128D0
EIJ(1,10)=1.00235D0
EIJ(1,9)=1.00305D0
EIJ(1,11)=1.01293D0
EIJ(1,12)=0.999758D0
EIJ(1,13)=0.988563D0
EIJ(6,7)=1.02274D0
UIJ(6,7)=0.835058D0
KIJ(6,7)=0.982361D0
GIJ(6,7)=0.982746D0
EIJ(2,6)=0.97012D0
UIJ(2,6)=0.816431D0
KIJ(2,6)=1.00796D0
EIJ(3,6)=0.945939D0
UIJ(3,6)=0.915502D0
EIJ(6,17)=1.08632D0
UIJ(6,17)=0.408838D0
KIJ(6,17)=1.03227D0
EIJ(6,18)=1.00571D0
EIJ(5,6)=0.946914D0
EIJ(4,6)=0.973384D0
UIJ(4,6)=0.993556D0
EIJ(6,10)=0.95934D0
EIJ(6,9)=0.94552D0
EIJ(6,11)=0.93788D0
EIJ(6,12)=0.935977D0
EIJ(6,13)=0.933269D0
EIJ(2,7)=0.925053D0
UIJ(2,7)=0.96987D0
KIJ(2,7)=1.00851D0
GIJ(2,7)=0.370296D0
EIJ(3,7)=0.960237D0
EIJ(7,17)=1.28179D0
EIJ(7,18)=1.5D0
UIJ(7,18)=0.9D0
EIJ(5,7)=0.906849D0
EIJ(4,7)=0.897362D0
EIJ(7,10)=0.726255D0
EIJ(7,9)=0.859764D0
EIJ(7,11)=0.766923D0
EIJ(7,12)=0.782718D0
EIJ(7,13)=0.805823D0
EIJ(2,3)=1.03502D0
UIJ(2,3)=1.0805D0
KIJ(2,3)=1.00046D0
EIJ(2,17)=1.16446D0
UIJ(2,17)=1.61666D0
KIJ(2,17)=1.02034D0
UIJ(2,5)=1.25D0
EIJ(2,4)=1.01306D0
UIJ(2,4)=1.25D0
UIJ(2,10)=1.25D0
EIJ(2,9)=1.00532D0
UIJ(2,9)=1.25D0
EIJ(3,17)=1.034787D0
EIJ(3,4)=1.0049D0
EIJ(17,18)=1.1D0
EIJ(5,17)=1.3D0
EIJ(4,17)=1.3D0
RETURN
END
SUBROUTINE PARMI1
IMPLICIT REAL*8(A-H,O-Z)
REAL*8 KI,KIJ,KM
INTEGER GN,QN,FN
DIMENSION EIJM(19,19),GIJM(19,19)
COMMON/Y1/Y(19)/NC1/NC/EFI/EI(19),KI(19),GI(19),QI(19),FI(19)
*/INTER1/EIJ(19,19),UIJ(19,19),KIJ(19,19),GIJ(19,19)
*/KM/KM/COEF1/B1(13),C1(53)/AN/AN(53)
*/GQFN/GN(53),QN(53),FN(53)/UN/UN(53)
DO 1 I=1,NC
EIJM(I,I)=EI(I)
GIJM(I,I)=GI(I)
DO 1 J=1,NC
IF(I.GE.J) GO TO 1
EIJM(I,J)=EIJ(I,J)*(EI(I)*EI(J))**.5
GIJM(I,J)=GIJ(I,J)*(GI(I)+GI(J))/2.
1 CONTINUE
KM=0D0
UM=0D0
KM=0D0
UM=0D0
GM=0D0
QM=0D0
FM=0D0
DO 3 I=1,NC
KM=KM+Y(I)*KI(I)**2.5
UM=UM+Y(I)*EI(I)**2.5
GM=GM+Y(I)*GI(I)
QM=QM+Y(I)*QI(I)
3 FM=FM+Y(I)**2*FI(I)
KM=KM*KM
UM=UM*UM
DO 5 I=1,NC-1
DO 5 J=I+1,NC
UM=UM+2.*Y(I)*Y(J)*(UIJ(I,J)**5-1D0)*(EI(I)*EI(J))**2.5
GM=GM+2.*Y(I)*Y(J)*(GIJ(I,J)-1D0)*(GI(I)+GI(J))
5 KM=KM+2.*Y(I)*Y(J)*(KIJ(I,J)**5-1D0)*(KI(I)*KI(J))**2.5
KM=KM**.6
UM=UM**.2
DO 7 N=1,13
B1(N)=0D0
DO 9 I=1,NC
9 В1(N)=B1(N)+Y(I)*Y(I)(GIJM(I,I)+ 1D0-GN(N))**GN(N)*
*(QI(I)*QI(I)+1D0-QN(N))**QN(N)*(FI(I)+1D0-FN(N))*FN(N)*
*EIJM(I,I)"UN(N)*KI(I)*KI(I)*KI(I)
DO 11 I=1,NC-1
DO 11 J=I+1,NC
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
END
SUBROUTINE COMPL1(RO)
IMPLICIT REAL*8(A-H,O-Z)
REAL*8 KM
INTEGER BN,CN
COMMON/KM/KM/COEF2/B,C(53)/BCKN/BN(53),CN(53),KN(53)/AI1/AO,A1
ROR=KM*RO
S1=0D0
S2=0D0
S3=0D0
DO 1 N=8,53
EXP=DEXP(-CN(N)*ROR**KN(N))
IF(N.LE.13) S1=S1+C(N)
S2=S2+C(N)*(BN(N)-CN(N)*KN(N)*ROR**KN(N))*ROR**BN(N)*EXP 1
S3=S3+C(N)*(-CN(N)*KN(N)**2*KM*ROR**(KN(N)-1)*ROR**BN(N)*
*EXP+(BN(N)-CN(N)*KN(N)*ROR**KN(N))*BN(N)*KM*ROR**(BN(N)-1)*
*EXP-(BN(N)-CN(N)*KN(N)*ROR**KN(N))*ROR**BN(N)*EXP*CN(N)*KN(N)*
*KM*ROR**(KN(N)-1))AO1=B*RO-ROR*S1
AO=AO1+S2
A1=AO+AO1+RO*S3
RETURN
END
BLOCK DATA DCAGA8
IMPLICIT REAL*8(A-H,O-Z)
REAL*8 KD
INTEGER BN,CN,GN,QN,FN
COMMON/EFD/ED(19),KD(19),GD(19),QD(19),FD(19)
*/BCKN/BN(53),CN(53),KN(53)/UN/UN(53)
*/AN/AN(53)/GQFN/GN(53),QN(53),FN(53)
DATA ED/151.3183D0,244.1667D0,298.1183D0,337.6389D0,324.0689D0,
*99.73778D0,241.9606D0,296.355D0,370.6823D0,365.5999D0,
*402.8429D0,427.5391D0,450.6472D0,472.1194D0:488.7633D0,
*2.610111D0,26.95794D0,105.5348D0,122.7667D0/
DATA KD/.4619255D0,.5279209D0,.583749D0,.6341423D0,.6406937D0,
*.4479153D0,.4557489D0,.4618263D0,.6798307D0,.6738577D0,
*.7139987D0,.7503628D0,.7851933D0,.8157596D0,.8389542D0,
*.3589888D0,.3514916D0,.4533894D0,.4186954D0/
DATA GD/0D0,.0793D0,.141239D0,.281835D0,.256692D0,.027815D0,
*.189065D0,.0885D0,.366911D0,.332267D0,.432254D0,.512507D0,
*.576242D0,.648601D0,.716574D0,0D0,.034369D0,.038953D0,.021D0/
DATA QD/6*0D0,.69D0,12*0D0/,FD/16*0D0,1D0,2*0D0/
DATA AN/.1538326D0,1.341953D0,-2.998583D0,-.04831228D0,
*.3757965D0,-1.589575D0,-.05358847D0,2.29129D-9,1576724D0,
*-.4363864D0,-.04408159D0,-.003433888D0,.03205905D0,.02487355D0, -
*.07332279D0,-.001600573D0,.6424706D0,-.4162601D0,-.06689957D0,
*.2791795D0,-.6966051D0,-.002860589D0,-.008098836D0,3.150547D0,
*.007224479D0,-.7057529D0,.5349792D0,-.07931491D0-1.418465D0,
*-5.99905D-17,.1058402D0,.03431729D0,-.007022847D0,.02495587D0,
*.04296818D0,.7465453D0,-.2919613D0,7.294616D0,-9.936757D0,
*-.005399808D0,-.2432567D0,.04987016D0,.003733797D0,1.874951D0,
*.002168144D0,-.6587164D0,.000205518D0,.009776195D0,-.02048708D0,
*.01557322D0,.006862415D0,-.001226752D0,.002850906D0/
DATA BN/13*1,9*2,10*3,7*4,5*5,2*6,2*7,3*8,2*9/
DATA CN/7*0,6*1,2*0,7*1,0,9*1,2*0,5*1,0,4*1,0,1,0,6*1/
DATA KN/7*0,3,3*2,2*4,2*0,3*2,4*4,0,2*1,2*2,2*3,3*4,2*0,3*2,
*2*4,0,2*2,2*4,0,2,0,2,1,4*2/
DATA UN/0D0,.5D0,1D0,3.5D0,-.5D0,4.5D0,.5D0,-6D0,2D0,3D0,2*2D0,
*11D0,-.5D0,.5D0,0D0,4D0,6D0,21D0,23D0,22D0,-1D0,-.5D0,7D0,-1D0,
*6D0,4D0,1D0,9D0,-13D0,21D0,8D0,-.5D0,0D0,2D0,7D0,9D0,22D0,23D0,
*1D0,9D0,3D0,8D0,23D0,1.5D0,5D0,-.5D0,4D0,7D0,3D0,0D0,1D0,0D0/
DATA GN/4*0,2*1,13*0,1,3*0,1,2*0,3*1,16*0,1,2*0,1,0,1,2*0/
DATA QN/6*0,1,3*0,1,9*0,1,0,1,8*0,1,4*0,1,4*0,1,0,1,2*0,1,5*0,1/
DATA FN/7*0,1,13*0,1,2*0,1,4*0,1,23*0/
END
C ***********************************************************
С * *
С * Подпрограмма расчета коэффициента сжимаемости природного *
С * газа по уравнению состояния ВНИЦ СМВ. *
С * *
C ***********************************************************
SUBROUTINE VNIC(ICALC)
IMPLICIT REAL*8(A-H,O-Z)
REAL*8 LIJ(8,8)
DIMENSION VC(8),TC(8),PII(8),DIJ(8,8)
COMMON/PARCD/VCD(8),TCD(8),PIID(8)/ABIJ/AIJ(10,8),BIJ(10,8)
*/B/B(10,8)/RM/RM/Y/Y(8)/BM/BM(8)/NI/NI(8)/NC/NC/Z/Z
RM=8.31451DO
IF(ICALC.NE.1) GO TO 1
CALL COMPON
IF(Z.EQ.0D0) GO TO 133
CALL DDIJ(DIJ,LIJ)
DO 75 I=1,NC
TC(I)=TCD(NI(I))
VC(I)=BM(I)/VCD(NI(I))
PII(I)=PIID(NI(I))
DO 123 J=1,NC
IF(I.GE.J) GO TO 123
DIJ(I,J)=DIJ(NI(I),NI(J))
LIJ(I,J)=LIJ(NI(I),NI(J))
123 CONTINUE
75 CONTINUE
CALL PARMIX(DIJ,LIJ,TC,VC,PII,PIM)
DO 27 I=1,10
DO 27 J=l,8
27 B(I,J)=AIJ(I,J)+BIJ(I,J)*PIM
1 CALL PHASE
133 RETURN
END
SUBROUTINE COMPON
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION BMI(25),ROI(8),GI(8),YI(25)
COMMON/Y/Y(8)/BMM/BMM/BM/BM(8)/YI/YC(25)/NI/NI(8)/NC/NC/NPR/NPR
DATA BMI/16.043D0,30.07D0,44.097D0,2*58.123D0,28.0135D0,
*44.01D0,34.082D0,26.038D0,28.054D0,42.081D0,3*72.15D0,
*86.177D0,78.114D0,100.204D0,92.141D0,114.231D0,128.259D0,
*142.286D0,4.0026D0,2.0159D0,28.01D0,31.9988D0/
DATAR0I/0.6682D0,1.2601D0,1.8641D0,2.4956D0,2.488D0,
*1.1649D0,1.8393D0,1.4311D0/
DO 100 I=1,25
100 YI(I)=YC(I)
IF(NPR.EQ.1) GO TO 333
BMM=0D0
DO 3333 I=1,25
3333 ВММ=ВММ+YI(I)*ВМI(I)
333 YS=0D0
DO 55 I=9,25
YS=YS+YI(I)
55 CONTINUE
YS1=0D0
DO 67 I=12,21
67 YS1=YS1+YI(I)
YS2=0D0
DO 69 1=22,25
69 YS2=YS2+YI(I)
YI(2)=YI(2)+YI(9)+YI(10)
YI(3)=YI(3)+YI(11)
YI(4)=YI(4)+YS1
YS3=YI(4)+YI(5)
IF(NPR.EQ.1.AND.YI(5).LT.0.01D0.AND.YS3.LT.0.03D0) YI(4)=YS3
IF(NPR.EQ.1.AND.YI(5).LT.0.01D0.AND.YS3.LT.0.03D0) YI(5)=0D0
IF(NPR.EQ.0.AND.Y1(5).LT.0.01D0.AND.YS3.LE.0.03D0) YI(4)=YS3
IF(NPR.EQ.0.AND.YI(5).LT.0.01D0.AND.YS3.LE.0.03D0)YI(5)=0D0
YI(6)=YI(6)+YS2
IF(NPR.EQ.0) GO TO 555
ROM=0D0
DO 7 I=1,8
7 ROM=ROM+YI(I)*ROI(I)
DO 9 I=1,8
9 GI(I)=YI(I)*ROI(I)/ROM
SUM=0D0
DO 11 1=1,8
11 SUM=SUM+GI(I)/BMI(I)
SUM=1./SUM
DO 13 I=1,8
13 YI(I)=GI(I)*SUM/BMI(I)
555 NC=0
YSUM=0D0
DO 155 I=1,8
IF(YI(I).EQ.0D0) GO TO 155
NC=NC+1
NI(NC)=I
Y(NC)=YI(I)
YSUM=YSUM+Y(NC)
BM(NC)=BMI(I)
155 CONTINUE
CALL MOLDOL(YI,YS)
DO 551 I=1,NC
551 Y(I)=Y(I)/YSUM
RETURN
END
SUBROUTINE MOLDOL(YI,YS)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION YI(25)
COMMON/Z/Z
Z=-1D0
IF(YI(1).LT.0.65D0.OR.YI(2).GT.0.15D0.OR.YI(3).GT.0.035D0.OR.
*YI(4).GT.0.015D0.OR.YI(5).GT.0.015D0.OR.YS.GT.0.01D0)Z=0D0
IF(YI(6).GT.0.2D0.OR.YI(7).GT.0.15D0.OR.YI(8).GT.0.3D0)Z=0D0
RETURN
END
SUBROUTINE DDIJ(DIJ,LIJ)
IMPLICIT REAL-8(A-H,O-Z)
REAL*8 LIJ(8,8)
DIMENSION DIJ(8,8)
DO 1 I=1,8
DO 1 J=l,8
LIJ(I,J)=0.D0
1 DIJ(I,J)=0.D0
DIJ(1,2)=0.036D0
DIJ(1,3)=0.076D0
DIJ(1,4)=0.121D0
DIJ(1,5)=0.129D0
DIJ(1,6)=0.06D0
DIJ(1,7)=0.074D0
DIJ(2,6)=0.106D0
DIJ(2,7)=0.093D0
DIJ(6,7)=0.022D0
DIJ(1,8)=0.089D0
DIJ(2,8)=0.079D0
DU(6,8)=0.211D0
DIJ(7,8)=0.089D0
LIJ(1,2)=-0.074D0
LIJ(1,3)=-0.146D0
LIJ(1,4)=-0.258D0
LIJ(1,5)=-0.222D0
LIJ(1,6)=-0.023D0
LIJ(1,7)=-0.086D0
LIJ(6,7)=-0.064D0
LIJ(7,8)=-0.062D0
RETURN
END
SUBROUTINE PARMIX(DIJ,LIJ,TC,VC,PII,PIM)
IMPLICIT REAL*8(A-H,O-Z)
REAL*8 LIJ(8,8)
DIMENSION Y(8),DIJ(8,8),VCIJ(8,8),TCIJ(8,8),V13(8),TC(8),VC(8),
*PII(8),PIIJ(8,8)
COMMON/PARCM/TCM,VCM/Y/Y/NC/NC/PCM/PCM
DO 1 I=1,NC
1 V13(I)=VC(I)**(1.DO/3.DO)
DO 3 I=1,NC
VCIJ(I,I)=VC(I)
PIIJ(I,I)=PII(I)
TCIJ(I,I)=TC(I)
DO 3 J=1,NC
IF(I.GE.J) GO TO 3
VCIJ(I,J)=(1.DO-LIJ(I,J))*((V13(I)+V13(J))/2.)**3
PIIJ(I,J)=(VC(I)*PII(I)+VC(J)*PII(J))/(VC(I)+VC(J))
TCU(I,J)=(1.D0-DIJ(I,J))*(TC(I)*TC(J))**0.5
VCIJ(J,I)=VCIJ(I,J)
PIIJ(J,I)=PIIJ(I,J)
TCIJ(J,I)=TCIJ(I,J)
3 CONTINUE
VCM=0.D0
PIM=0.D0
TCM=0.D0
DO 5 I=1,NC
DO 5 J=1,NC
VCM=VCM+Y(I)*Y(J)*VCIJ(I,J)
PIM=PIM+Y(I)*Y(J)*VCIJ(I,J)*PIIJ(I,J)
5 TCM=TCM+Y(I)*Y(J)*VCIJ(I,J)*TCIJ(I,J)**2
PIM=PIM/VCM
TCM=(TCM/VCM)**0.5
PCM=8.31451D-3*(0.28707D0-0.05559*PIM)*TCM/VCM
RETURN
END
SUBROUTINE PHASE
IMPLICIT REAL*8(A-H,O-Z)
COMMON/Z/Z/RM/RM/T/T/P/P/PCM/PCM/AI/AO,A1
IF(T.LT.250D0.OR.T.GT.340D0.OR.P.LE.0D0.OR.P.GT.12D0) THEN
Z=0D0
GO TO 134
ENDIF
PR=P/PCM
RO=9D3*P/(RM*T*(1.1*PR+0.7D0))
CALL FUN(RO)
CALL OMTAU(RO,T)
IF(Z.EQ.0D0) GO TO 134
Z=1.D0+AO
134 RETURN
END
С Подпрограмма, реализующая итерационный процесс определения
С плотности из уравнения состояния (метод Ньютона)
SUBROUTINE FUN(X)
IMPLICIT REAL*8(A-H,О-Z)
COMMON/P/P/RM/RM/T/T/AI/AO,A1
ITER=1
1 CONTINUE
NPRIZ=0
IF(ITER.NE.l) NPRIZ=1
CALL COMPL(X,T,NPRIZ)
Z=1.D0+AO
FX=1.D6*(P-(1.D-3*RM*T*Z*X))
F=1.D3*RM*T*(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 COMPL(X,T,NPRIZ)
RETURN
END
SUBROUTINE OMTAU(RO,T)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/PARCM/TCM,VCM/Z/Z
Z=-1D0
TR=T/TCM
ROR=RO*VCM
IF(TR.LT.1.05D0) Z=0D0
IF(ROR.LT.0.D0.OR.ROR.GT.3.D0) Z=0D0
RETURN
END
SUBROUTINE COMPL(RO,T,NPRIZ)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION B(10,8),BK(10)
COMMON/PARCM/TCM,VCM/B/B/AI/AO,A1
IF(NPRIZ.NE.0) GO TO 7
TR=T/TCM
DO 1 I=1,10
BK(I)=0
DO 1 J=1,8
1 BK(I)=BK(I)+B(I,J)/TR**(J-1)
7 ROR=RO*VCM
AO=0.D0
A1=0.D0
DO 33 I=1,10
D=BK(I)*ROR**I
AO=AO+D
33 A1=A1+(I+1)*D
RETURN
END
BLOCK DATA BDVNIC
IMPLICIT REAL*8(A-H,O-Z)
COMMON/PARCD/VCD(8),TCD(8),PIID(8)/ABIJ/AIJ(10,8),BIJ(10,8)
DATA TCD/190.67D0,305.57D0,369.96D0,425.4D0,407.96D0,
*125.65D0,304.11D0,373.18D0/
DATA VCD/163.03D0,205.53D0,218.54D0,226.69D0,225.64D0,
*315.36D0,466.74D0,349.37D0/
DATA PIID/0.0006467D0,0.1103D0,0.1764D0,0.2213D0,0.2162D0,
*0.04185D0,0.2203D0,0.042686D0/
DATA AIJ/.6087766D0,-.4596885D0,1.14934D0,-.607501D0,
*-.894094D0,1.144404D0,-.34579D0,-.1235682D0,.1098875D0,
*-.219306D-1,-1.832916D0,4.175759D0,-9.404549D0,10.62713D0,
*-3.080591D0,-2.122525D0,1.781466D0,-.4303578D0,-.4963321D-1,
*.347496D-1,1.317145D0,-10.73657D0,23.95808 D0,-31.47929D0,
* 18.42846D0,-4.092685D0,-. 1906595D0,.4015072D0,-.1016264D0,