-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathBLOCK4.FOR
324 lines (321 loc) · 11.9 KB
/
BLOCK4.FOR
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
C.............................................................PSEUDO
SUBROUTINE PSEUDO(NP,IFATAL,IOCODE)
C MACHINE DEPENDENT INCLUDE STATEMENT
INCLUDE 'PARAMS.FOR'
C GAS PROP BY CORRELATIONS
REAL KROT,KROGT,KRWT,KRGT,MUOT,MUWT,MUGT
COMMON /BUBBLE/ PBO,VSLOPE(LP8),BSLOPE(LP8),RSLOPE(LP8),PMAXT,
& IREPRS,MPGT(LP8),
& RHOSCO(LP8),RHOSCW(LP8),RHOSCG(LP8),MSAT(LP7),MPOT(LP8),
& MPWT(LP8),PBOT(LP1,LP2,LP3),PBOTN(LP1,LP2,LP3)
COMMON /SPVT/ SAT(LP7,LP9),KROT(LP7,LP9),KRWT(LP7,LP9),
& BGT(LP7,LP9),
& KRGT(LP7,LP9),ITHREE(LP7),RSOT(LP7,LP9),BWPT(LP7,LP9),
& PCOWT(LP7,LP9),PCGOT(LP7,LP9),KROGT(LP7,LP9),SWR(LP7),
& POT(LP7,LP9),MUOT(LP7,LP9),BOT(LP7,LP9),BOPT(LP7,LP9),
& RSOPT(LP7,LP9),PWT(LP7,LP9),MUWT(LP7,LP9),BWT(LP7,LP9),
& PGT(LP7,LP9),MUGT(LP7,LP9),
& BGPT(LP7,LP9),CRT(LP7,LP9),IPVT(LP1,LP2,LP3),IROCK(LP1,LP2,LP3),
& NROCK,NPVT,PSIT(LP7,LP9),PRT(LP7,LP9),WOROCK(LP7),GOROCK(LP7)
DIMENSION FRCI(12),PRSCI(12),RMWTI(12),TEMCI(12)
& ,ZEDD(LP8,LP9),CMPGD(LP8,LP9)
DATA NCOMP/12/
DATA TEMCI /672.37,
A547.57, 227.27, 343.04, 549.76, 665.68,
B734.65, 765.32, 828.77, 845.37, 913.37, 1023.89/
DATA PRSCI /1306.0,
A1071.0, 493.00, 667.80, 707.80, 616.30,
B529.10, 550.70, 490.40, 488.60, 436.90, 360.60/
200 FORMAT(/)
220 FORMAT(1H0,'RESERVOIR TEMPERATURE, DEGREES F =',F10.2,/
& 1H , 'GAS GRAVITY =',F10.4,/
& 1H , 'PSEUDO-CRITICAL TEMP., DEGREES R =',F10.2,/
& 1H , 'PSEUDO-CRITICAL PRESSURE, PSIA =',F10.2,/
& 1H , 'MOLE PERCENT - HYDROGEN SULFIDE =',F10.2,/
& 1H , 'MOLE PERCENT - CARBON DIOXIDE =',F10.2,/
& 1H , 'MOLE PERCENT - NITROGEN =',F10.2,//)
225 FORMAT(1HO,'PRESSURE',4X,'PSEUDO-PRESS',3X,'Z-FACTOR',4X,
&'GAS FVF',4X,'GAS COMP',2X,'GAS VISG',/,
&2X,'(PSIA)',5X,'(PSIA**2/CP)',14X,
&'(RB/SCF)',4X,'(1/PSIA)',5X,'(CP)',/)
241 FORMAT(1H ,F7.1,5X,E13.7,3X,F6.3,2(2X,E10.4),2X,F7.6)
READ(20,*) KODEA,MPGT(NP),TEM,SPG
READ(20,*) (FRCI(I),I=1,12)
IF(KODEA.LT.4) GO TO 560
READ(20,*) PRSCI(12),TEMCI(12),RMWTI(12)
560 CONTINUE
NDATA=MPGT(NP)
NFIRST=1
NLAST=NDATA
PREF=14.7
DELP=(PMAXT-PREF)/(MPGT(NP)-1.)
CNCH2S=FRCI(1)*100.
CNCCO2=FRCI(2)*100.
CNCN2=FRCI(3)*100.
IF(KODEA.GT.2) GO TO 2
TEMPC=170.491 + 307.344*SPG
PRSPC=709.604 - 58.718*SPG
GO TO 51
2 TEMPC=0.0
PRSPC=0.0
DO 5 J=1,NCOMP
TEMPC=TEMPC+TEMCI(J)*FRCI(J)
PRSPC=PRSPC+PRSCI(J)*FRCI(J)
5 CONTINUE
51 CONTINUE
WRITE(IOCODE,201)
201 FORMAT(///,8X,4('*'),' REAL GAS PROPERTIES ',4('*'),//)
WRITE(IOCODE,203)
WRITE(IOCODE,204) FRCI(1),FRCI(2),FRCI(3)
WRITE(IOCODE,205)
WRITE(IOCODE,206) (FRCI(I),I=4,NCOMP)
203 FORMAT(/,1X,'GAS COMPOSITION (MOLE FRACTION):',//,
& 3X,'H2S',4X,'CO2',5X,'N2')
204 FORMAT(2X,F5.4,2(2X,F5.4),/)
205 FORMAT(4X,'C1',5X,'C2',5X,'C3',4X,'IC4',4X,'NC4',
& 4X,'IC5',4X,'NC5',5X,'C6',4X,'C7+',/)
206 FORMAT(1X,F6.4,8(2X,F5.4),/)
WRITE(IOCODE,200)
WRITE(IOCODE,220) TEM,SPG,TEMPC,PRSPC,CNCH2S,CNCCO2,CNCN2
WRITE(IOCODE,225)
PRSDEL=DELP
PRS=PREF
PRSSI=0.0
PRSSI1=0.0
PRSSI2=0.0
DO 20 I=1,NDATA
BGT(NP,I)=0.
IF(PRS.LE.0.0) GO TO 15
TEMPRD=(TEM+460.)/TEMPC
PRSPRD=PRS/PRSPC
CALL ZANDC(TEM,TEMPC,PRS,PRSPC,CNCH2S,
& CNCCO2,ZED,CMPG,IERR)
IFATAL=IERR+IFATAL
IF(IERR.GT.0) WRITE(IOCODE,1010) IERR
1010 FORMAT(/5X,5('-'),'ZANDC ERRORS (IERR =',I5,') HAVE',
& ' OCCURRED IN THE DEFAULT GAS PROPERTIES CALCULATION.',
& /10X,'CHECK INPUT GAS DATA AND USER-SPECIFIED OPTIONS.',/)
CALL VISCY(TEMPRD,PRSPRD,SPG,TEM,CNCH2S,
& CNCCO2,CNCN2,VISGR,VISG,IERR)
IFATAL=IERR+IFATAL
IF(IERR.GT.0) WRITE(IOCODE,1020) IERR
1020 FORMAT(/5X,5('-'),'VISCY ERRORS (IERR =',I5,') HAVE',
& ' OCCURRED IN THE DEFAULT GAS PROPERTIES CALCULATION.',
& /10X,'CHECK INPUT GAS DATA AND USER-SPECIFIED OPTIONS.',/)
PRSSI2=2.*PRS/(VISG*ZED)
PRSSI=(PRSSI1+PRSSI2)/2.*PRSDEL + PRSSI
15 IF(PRS.LE.PREF ) PRSSI=0.0
PSIT(NP,I)=PRSSI
ZEDD(NP,I)=ZED
CMPGD(NP,I)=CMPG
MUGT(NP,I)=VISG
PGT(NP,I)=PRS
PRS=PRS+PRSDEL
PRSSI1=PRSSI2
IF(PGT(NP,I).EQ.0.0) GO TO 16
C!!!!!!!!!!!!!!!!!!!!! BEGIN BUG FIX !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C!!!!!!!!!!!!!!!! CLAY KIMBRELL CALLED 06/09/94 - SAID THEY HAD FOUND BUG HERE-
C!!!!!!!!!!!!!!!! .004675 should be .005035 - I will just accept this and change
!!!!!!!!!!!!!!!!! it now - will not recompile today - 06/16/94
C** BGT FACTOR 0.004675 = 14.7/((460+60)*5.6146)
C-BUG BGT(NP,I)=(TEM+460.)*ZEDD(NP,I)*0.004675/PGT(NP,I)
BGT(NP,I)=(TEM+460.)*ZEDD(NP,I)*0.005035/PGT(NP,I)
C!!!!!!!!!!!!!!!!!!!!!!!! END BUG FIX !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16 CONTINUE
IF(I.LT.NFIRST.OR.I.GT.NLAST) GO TO 20
WRITE(IOCODE,241) PGT(NP,I),PSIT(NP,I),ZEDD(NP,I),
& BGT(NP,I),CMPGD(NP,I),MUGT(NP,I)
20 CONTINUE
DO 30 I=1,NDATA
BGT(NP,I)=BGT(NP,I)*5.6146
30 CONTINUE
RETURN
END
C...............................................................REPRS1
SUBROUTINE REPRS1(I,J,K)
C CALCULATION OF REPRESSURIZATION PVT PROPERTIES
C MACHINE DEPENDENT INCLUDE STATEMENT
INCLUDE 'PARAMS.FOR'
REAL KROT,KROGT,KRWT,KRGT,MUWT,MUOT,MUGT,KX,KY,KZ
& ,MUO,MUW,MUG,KRO,KRW,MBEWI,MBEGI,MCFGI
& ,MBEO,MBEW,MBEG,MCFG,MBEOI,MIN,MCFG1,MCFGT
COMMON /BUBBLE/ PBO,VSLOPE(LP8),BSLOPE(LP8),RSLOPE(LP8),PMAXT,
& IREPRS,MPGT(LP8),
& RHOSCO(LP8),RHOSCW(LP8),RHOSCG(LP8),MSAT(LP7),MPOT(LP8),
& MPWT(LP8),PBOT(LP1,LP2,LP3),PBOTN(LP1,LP2,LP3)
COMMON /SARRAY/ PN(LP1,LP2,LP3),IOCODE,IDMAX,
& SON(LP1,LP2,LP3),SWN(LP1,LP2,LP3),SGN(LP1,LP2,LP3),
& A1(LP1,LP2,LP3),A2(LP1,LP2,LP3),A3(LP1,LP2,LP3),
& SUM(LP1,LP2,LP3),GAM(LP1,LP2,LP3),QS(LP1,LP2,LP3)
COMMON /SPRTPS/ P(LP1,LP2,LP3),SO(LP1,LP2,LP3),SW(LP1,LP2,LP3),
& SG(LP1,LP2,LP3)
COMMON /SPVT/ SAT(LP7,LP9),KROT(LP7,LP9),KRWT(LP7,LP9),
& BGT(LP7,LP9),
& KRGT(LP7,LP9),ITHREE(LP7),RSOT(LP7,LP9),BWPT(LP7,LP9),
& PCOWT(LP7,LP9),PCGOT(LP7,LP9),KROGT(LP7,LP9),SWR(LP7),
& POT(LP7,LP9),MUOT(LP7,LP9),BOT(LP7,LP9),BOPT(LP7,LP9),
& RSOPT(LP7,LP9),PWT(LP7,LP9),MUWT(LP7,LP9),BWT(LP7,LP9),
& PGT(LP7,LP9),MUGT(LP7,LP9),
& BGPT(LP7,LP9),CRT(LP7,LP9),IPVT(LP1,LP2,LP3),IROCK(LP1,LP2,LP3),
& NROCK,NPVT,PSIT(LP7,LP9),PRT(LP7,LP9),WOROCK(LP7),GOROCK(LP7)
COMMON /SSOLN/ BO(LP1,LP2,LP3),BW(LP1,LP2,LP3),BG(LP1,LP2,LP3),
& QO(LP1,LP2,LP3),QW(LP1,LP2,LP3),QG(LP1,LP2,LP3),
& GOWT(LP1,LP2,LP3),GWWT(LP1,LP2,LP3),GGWT(LP1,LP2,LP3),
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
CCC & OW(LP4,LP2,LP3),OE(LP4,LP2,LP3),WW(LP4,LP2,LP3),WE(LP4,LP2,LP3),
CCC & OS(LP1,LP5,LP3),ON(LP1,LP5,LP3),WS(LP1,LP5,LP3),WN(LP1,LP5,LP3),
CCC & OT(LP1,LP2,LP6),OB(LP1,LP2,LP6),WT(LP1,LP2,LP6),WB(LP1,LP2,LP6),
& O1(LP4,LP2,LP3),W1(LP4,LP2,LP3),
& O2(LP1,LP5,LP3),W2(LP1,LP5,LP3),
& O3(LP1,LP2,LP6),W3(LP1,LP2,LP6),
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
& QOWG(LP1,LP2,LP3),VP(LP1,LP2,LP3),CT(LP1,LP2,LP3)
COMMON /TSTDAT/ IFATAL,IWARN
C REPRESSURIZATION ASSUMING EQUILIBRATION IN ENTIRE GRID BLOCK
PP=P(I,J,K)
IF(P(I,J,K).GT.PBOT(I,J,K)) PP=PBOT(I,J,K)
IPVTR=IPVT(I,J,K)
CALL INTERP(IPVTR,POT,BOT,MPOT(IPVTR),PP,BBO)
CALL INTERP(IPVTR,POT,RSOT,MPOT(IPVTR),PP,RSO)
CALL INTERP(IPVTR,PGT,BGT,MPGT(IPVTR),PP,BBG)
SGOSOG=0.
IF(SO(I,J,K).GT.0.) SGOSOG=SG(I,J,K)*BBO/(SO(I,J,K)*BBG)
RSONEW=RSO+SGOSOG
CALL INTERP(IPVTR,RSOT,POT,MPOT(IPVTR),RSONEW,PBONEW)
PBOT(I,J,K)=PBONEW
RETURN
END
C.............................................................GAUS1D
SUBROUTINE GAUS1D(IOCODE,IFATAL,NX,NY,NZ)
C MACHINE DEPENDENT INCLUDE STATEMENT
INCLUDE 'PARAMS.FOR'
C SOLVE 1D PROBLEMS
COMMON /COEF/ AW(LP1,LP2,LP3),AE(LP1,LP2,LP3),AN(LP1,LP2,LP3),
& AS(LP1,LP2,LP3),AB(LP1,LP2,LP3),AT(LP1,LP2,LP3),E(LP1,LP2,LP3),
& B(LP1,LP2,LP3)
COMMON /SPRTPS/ P(LP1,LP2,LP3),SO(LP1,LP2,LP3),SW(LP1,LP2,LP3),
& SG(LP1,LP2,LP3)
COMMON /TRIDI/ UM(LP15),AZL(LP15),BZL(LP15),CZL(LP15),DZL(LP15),
& UZL(LP15),X(LP15),BETA(LP15),GAMMA(LP15),W(LP15)
IF(NX.GT.1.OR.NY.GT.1.OR.NZ.GT.1) GO TO 50
I=1
J=1
K=1
IF(E(I,J,K).EQ.0.) GO TO 20
P(I,J,K)=B(I,J,K)/E(I,J,K)
RETURN
20 WRITE(IOCODE,30)
30 FORMAT(/5X,5('-'),'COEF E(1,1,1) = 0 FOR',
& ' NX=NY=NZ=1 CASE; PLEASE CHECK INPUT DATA.',/)
IFATAL=IFATAL+1
RETURN
50 CONTINUE
IF(NX.GT.1.AND.NY.GT.1) GO TO 3000
IF(NX.GT.1.AND.NZ.GT.1) GO TO 3000
IF(NY.GT.1.AND.NZ.GT.1) GO TO 3000
IF(NX.EQ.1) GO TO 1000
J=1
K=1
DO 100 I=1,NX
AZL(I)=AW(I,J,K)
BZL(I)=E(I,J,K)
CZL(I)=AE(I,J,K)
DZL(I)=B(I,J,K)
100 CONTINUE
C THIS ALGORITHM SOLVES THE TRIDIAGONAL SYSTEM
C GENERATED BY THE SYSTEM OF N EQUATIONS
C AZL(I)*U(I-1) + BZL(I)*U(I) + CZL(I)*U(I+1) = DZL(I)
BETA(1)=BZL(1)
GAMMA(1)=DZL(1)/BZL(1)
NXM=NX-1
C COMPUTE FORWARD SOLUTION.
DO 5010 ITRI=1,NXM
W(ITRI)=CZL(ITRI)/BETA(ITRI)
ITRIP=ITRI+1
5010 BETA(ITRIP)=BZL(ITRIP)-AZL(ITRIP)*W(ITRI)
DO 5020 ITRI=2,NX
ITRIM=ITRI-1
5020 GAMMA(ITRI)=(DZL(ITRI)-AZL(ITRI)*GAMMA(ITRIM))/BETA(ITRI)
C COMPUTE BACK SOLUTION
UZL(NX)=GAMMA(NX)
DO 5030 JTRI=1,NXM
ITRI=NX-JTRI
ITRIP=ITRI+1
5030 UZL(ITRI)=GAMMA(ITRI)-W(ITRI)*UZL(ITRIP)
DO 300 I=1,NX
P(I,J,K)=UZL(I)
300 CONTINUE
RETURN
1000 IF(NZ.EQ.1) GO TO 2000
I=1
J=1
DO 1100 K=1,NZ
AZL(K)=AT(I,J,K)
BZL(K)=E(I,J,K)
CZL(K)=AB(I,J,K)
DZL(K)=B(I,J,K)
1100 CONTINUE
C THIS ALGORITHM SOLVES THE TRIDIAGONAL SYSTEM
C GENERATED BY THE SYSTEM OF N EQUATIONS
C AZL(K)*U(K-1) + BZL(K)*U(K) + CZL(K)*U(K+1) = DZL(K)
BETA(1)=BZL(1)
GAMMA(1)=DZL(1)/BZL(1)
NZM=NZ-1
C COMPUTE FORWARD SOLUTION.
DO 5110 ITRI=1,NZM
W(ITRI)=CZL(ITRI)/BETA(ITRI)
ITRIP=ITRI+1
5110 BETA(ITRIP)=BZL(ITRIP)-AZL(ITRIP)*W(ITRI)
DO 5120 ITRI=2,NZ
ITRIM=ITRI-1
5120 GAMMA(ITRI)=(DZL(ITRI)-AZL(ITRI)*GAMMA(ITRIM))/BETA(ITRI)
C COMPUTE BACK SOLUTION
UZL(NZ)=GAMMA(NZ)
DO 5130 JTRI=1,NZM
ITRI=NZ-JTRI
ITRIP=ITRI+1
5130 UZL(ITRI)=GAMMA(ITRI)-W(ITRI)*UZL(ITRIP)
DO 1300 K=1,NZ
P(I,J,K)=UZL(K)
1300 CONTINUE
RETURN
2000 CONTINUE
I=1
K=1
DO 2100 J=1,NY
AZL(J)=AS(I,J,K)
BZL(J)=E(I,J,K)
CZL(J)=AN(I,J,K)
DZL(J)=B(I,J,K)
2100 CONTINUE
C THIS ALGORITHM SOLVES THE TRIDIAGONAL SYSTEM
C GENERATED BY THE SYSTEM OF N EQUATIONS
C AZL(J)*U(J-1) + BZL(J)*U(J) + CZL(J)*U(J+1) = DZL(J)
BETA(1)=BZL(1)
GAMMA(1)=DZL(1)/BZL(1)
NYM=NY-1
C COMPUTE FORWARD SOLUTION.
DO 5210 ITRI=1,NYM
W(ITRI)=CZL(ITRI)/BETA(ITRI)
ITRIP=ITRI+1
5210 BETA(ITRIP)=BZL(ITRIP)-AZL(ITRIP)*W(ITRI)
DO 5220 ITRI=2,NY
ITRIM=ITRI-1
5220 GAMMA(ITRI)=(DZL(ITRI)-AZL(ITRI)*GAMMA(ITRIM))/BETA(ITRI)
C COMPUTE BACK SOLUTION
UZL(NY)=GAMMA(NY)
DO 5230 JTRI=1,NYM
ITRI=NY-JTRI
ITRIP=ITRI+1
5230 UZL(ITRI)=GAMMA(ITRI)-W(ITRI)*UZL(ITRIP)
DO 2300 J=1,NY
P(I,J,K)=UZL(J)
2300 CONTINUE
RETURN
3000 CONTINUE
IFATAL=IFATAL+1
WRITE(IOCODE,3100) NX,NY,NZ
3100 FORMAT(/5X,5('-'),'GAUS1D INCORRECTLY CALLED',
& ' WHEN NX,NY,NZ = ',3I5,/)
RETURN
END