IF(NL1.EQ.6)
*WRITE(1,F)PI(J),ZPP(1),ZPP(2),ZPP(3),ZPP(4),ZPP(5),ZPP(6)
NW=NW+1
133CONTINUE
IF(NW.EQ.20.AND.NYST.EQ.0) THEN
IF(J.EQ.NP.AND.I+NL-1.EQ.NT) GO TO 29
WRITE(*,7)
PAUSE ’ Для продолжения вывода нажмите <ВВОД> ’
WRITE(*,100)
NW=0
WRITE(1,7)
IF(NL.GT.1)WRITE(1,17)LIN2(1),(LIN1(K),K=1,NL-1)
IF(NL.EQ.l) WRITE(1,17)LIN2(1)
WRITE(1,19)AT(NL)
IF(NL.GT.1)WRITE(1,21)LIN4,(LIN2(K),K=1,NL-1)
IF(NL.EQ.l) WRITE(1,21)LIN4
WRITE(1,23)(TI(K),K=I,I+NL-1)
WRITE(1,17)(LIN3(K),K=1,NL)
NW=NW+6
ENDIF
IF(NW.EQ.54.AND.NYST.NE.0) THEN
IF(J.EQ.NP.AND.I+NL-1.EQ.NT) GO TO 29
WRITE(1,7)
WRITE(*,7)
IF(NYST.EQ.l) PAUSE
*’ Для продолжения вывода вставьте бумагу и нажмите <ВВОД> ’
NW=0
IF(NL.GT.1) WRITE(1,17)LIN2(1),(LIN1(K),K=1,NL-1)
IF(NL.EQ.l) WRITE(1,17)LIN2(1)
WRITE(1,19)AT(NL)
IF(NL.GT.1) WRITE(1,21)LIN4,(LIN2(K),K=1,NL-1)
IF(NL.EQ.l) WRITE(1,21)LIN4
WRITE(1,23)(TI(K),K=I,I+NL-1)
WRITE(1,17)(LIN3(K),K=1,NL)
NW=NW+6
ENDIF
25CONTINUE
15CONTINUE
29CLOSE(l)
WRITE(*,7)
PAUSE ’ Вывод завершен, для продолжения работы нажмите <ВВОД> ’
WRITE(*,66)
66FORMAT(/’ Назначить другое устройство вывода ?’,
*’, 0 - нет, 1 - да ’)
READ(*,*)NBOLB
IF(NBOLB.EQ.l) GO TO 22
RETURN
END
C**********************************************************
С**
С* Подпрограмма расчета коэффициента сжимаемости природного*
С*газа по модифицированному методу NX19.*
C**********************************************************
SUBROUTINE NX19(YA,YY)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/NCONT/NCONT/YA/Y(2)/RON/RON
Y(1)=YA
Y(2)=YY
CALL PTCONT
IF(NCONT.EQ.l) GO TO 134
CALL EA
CALL PHASEA
134RETURN
END
SUBROUTINE PTCONT
IMPLICIT REAL*8(A-H,O-Z)
COMMON/NCONT/NCONT/Z/Z/P/P/T/T/YA/Y(2)/RON/RON
NCONT=0
IF(RON.LT.0.66D0.OR.RON.GT.1D0) NCONT=1
IF(Y(1).GT.0.2D0.OR.Y(2).GT.0.15D0) NCONT=l
IF(P.LE.0.D0.OR.T.LE.0.D0) NCONT=1
IF(T.LT.250.D0.OR.T.GT.340.D0) NCONT=1
IF(P.GT.12.D0) NCONT=1
IF(NCONT.EQ.1) Z=0D0
RETURN
END
SUBROUTINE EA
IMPLICIT REAL*8(A-H,O-Z)
COMMON/T/T/YA/Y(2)/RON/RON/P/P/PT/PA,TA/BI/B1,B2/T0/T0
PCM=2.9585*(1.608D0-0.05994*RON+Y(2)-.392*Y(1))
TCM=88.25*(0.9915D0+1.759*RON-Y(2)-1.681*Y(1))
PA=0.6714*P/PCM+0.0147
TA=0.71892*T/TCM+0.0007
DTA=TA-1.09D0
F=0D0
IF(PA.GE.0D0.AND.PA.LT.2D0.AND.DTA.GE.0D0.AND.DTA.LT.0.3D0)
F=75D-5*PA**2.3/DEXP(20.*DTA)+
*11D-4*DTA**0.5*(PA*(2.17D0-PA+1.4*DTA**0.5))**2
IF(PA.GE.0D0.AND.PA.LT.1.3D0.AND.DTA.GE.-0.25D0.AND.DTA.LT.0D0)
*F=75D-5*PA**2.3*(2D0-DEXP(20.*DTA))+
*1.317*PA*(1.69D0-PA**2)*DTA**4
IF(PA.GE.1.3D0.AND.PA.LT.2D0.AND.DTA.GE.-0.21D0.AND.DTA.LT.0D0)
*F=75D-5*PA**2.3*(2D0-DEXP(20.*DTA))+
*0.455*(1.3D0-PA)*(1.69*2.D0**1.25-PA**2)*(DTA*(0.03249D0+
*18.028*DTA**2)+DTA**2*(2.0167D0+DTA**2*(42.844D0+200.*DTA**2)))
T1=:TA**5/(TA**2*(6.60756*TA-4.42646D0)+3.22706D0)
T0=(TA**2*(1.77218D0-0.8879*TA)+0.305131D0)*T1/TA**4
B1=2.*T1/3.-TO**2
B0=T0*(T1-T0**2)+0.1*T1*PA*(F-1D0)
B2=(B0+(B0**2+B1**3)**0.5)**(1D0/3D0)
RETURN
END
SUBROUTINE PHASEA
IMPLICIT REAL*8(A-H,O-Z)
COMMON/Z/Z/PT/PA,TA/BI/B1,B2/T0/T0
Z=(1D0+0.00132/TA**3.25)**2*0.1*PA/(B1/B2-B2+T0)
RETURN
END
C*************************************************************
С**
С*Подпрограмма расчета коэффициента сжимаемости природного*
С*газа по модифицированному уравнению состояния GERG-91.*
С**
C*************************************************************
$NOTRUNCATE
SUBROUTINE GERG2(ICALC,YA,YY)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/T/T1/P/PRESS/RON/RON/Z/Z
COMMON/XBLOK/X1,X2,X3,X11,X12,X13,X22,X23,X33
COMMON/MBLOK/GM2,GM3,FA,FB,TO,R
DATABMO/.0838137D0/,BM1/-.00851644D0/,WD0/134.2153D0/,
*WD1/1067.943D0/
Z=-1D0
IF(ICALC.EQ.2) GO TO 3
X2=YA
X3=YY
IF(RON.LT.0.66D0.OR.RON.GT.1D0)Z=0D0
IF(X2.LT.0D0.OR.X2.GT.0.2D0)Z=0D0
IF(X3.LT.0D0.OR.X3.GT.0.15D0) Z=0D0
IF(Z.EQ.0D0) GO TO 133
X1=1D0-X2-X3
Х11=Х1*Х1
X12=X1*X2
X13=X1*X3
X22=X2*X2
X23=X2*X3
X33=X3*X3
Z=1D0-(.0741*RON-.006D0-.063*YA-.0575*YY)**2
BMNG=24.05525*Z*RON
Y1=1D0-YA-YY
BMY=(BMNG-28.0135*YA-44.01*YY)/Y1
СРасчет теплоты сгорания эквивалентного углеводорода (Н)
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)
133RETURN
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))
123CONTINUE
75CONTINUE
CALL PARMI1
3CALL PHASE1
133RETURN
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
100YI(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
7SUM=SUM+YI(I)/ZNI(I)
DO 9 I=1,25
9YI(I)=YI(I)/ZNI(I)/SUM
5YI(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)
11CONTINUE
CALL MOLDOl(YI)
DO 13 I=1,NC
13Y(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
1YS=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
1GIJ(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.
1CONTINUE
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)
3FM=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))
5KM=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
7CONTINUE
DO 13 N=8,53
13C1(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
1B=B+AN(N)/T**UN(N)*B1(N)
DO 3 N=8,53
3C(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
134RETURN
END
СПодпрограмма, реализующая итерационный процесс определения
Сплотности из уравнения состояния (метод Ньютона)
SUBROUTINE FUN1(X)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/P/P/RM/RM/T/T/AI1/AO,A1
ITER=1
1CONTINUE
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
4CALL 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))
123CONTINUE
75CONTINUE
CALL PARMIX(DIJ,LIJ,TC,VC,PII,PIM)
DO 27 I=1,10
DO 27 J=l,8
27B(I,J)=AIJ(I,J)+BIJ(I,J)*PIM
1CALL PHASE
133RETURN
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,