home *** CD-ROM | disk | FTP | other *** search
/ Hall of Fame / HallofFameCDROM.cdr / oilfield / spe-46-2.lzh / MAIN.FOR < prev    next >
Text File  |  1988-08-10  |  44KB  |  1,157 lines
  1. $DO66
  2.        PROGRAM BOASTII
  3. C      BOAST II: BLACK OIL APPLIED SIMULATION TOOL
  4. C      J.R. FANCHI
  5. C      APRIL 1987       RELEASE 1.2
  6. C      MACHINE DEPENDENT INCLUDE STATEMENT (SEE APPENDIX I OF THE
  7. C      BOAST II USER'S MANUAL FOR FURTHER COMMENTS)
  8.       PARAMETER (LP1=13,LP2=10,LP3=4)
  9.       PARAMETER (LP7=3,LP8=3,LP9=25,LP10=3,LP11=50,LP12=1000)
  10.       PARAMETER (LP14=1,LP15=13,LP17=5)
  11.       PARAMETER (LP4=LP1+1,LP5=LP2+1,LP6=LP3+1,LP13=LP1+LP2+LP3)
  12.       PARAMETER (LP19=LP4*LP2*LP3,LP20=LP1*LP5*LP3)
  13.       PARAMETER (LP21=LP1*LP2*LP6,LP22=LP1*LP2*LP3,LP23=LP11*LP3)
  14. C
  15.       CHARACTER*10 INDAT,OUTDAT,RESIN,RESOUT
  16.       CHARACTER*5 WELNAM
  17.       REAL KROT,KROGT,KRWT,KRGT,MUWT,MUOT,MUGT,KX,KY,KZ
  18.      & ,MUO,MUW,MUG,KRO,KRW,MBEWI,MBEGI,MCFGI
  19.      & ,MBEO,MBEW,MBEG,MCFG,MBEOI,MIN,MCFG1,MCFGT
  20. C
  21.        COMMON /BUBBLE/ PBO,VSLOPE(LP8),BSLOPE(LP8),RSLOPE(LP8),PMAXT,
  22.      & IREPRS,MPGT(LP8),
  23.      & RHOSCO(LP8),RHOSCW(LP8),RHOSCG(LP8),MSAT(LP7),MPOT(LP8),
  24.      & MPWT(LP8),PBOT(LP1,LP2,LP3),PBOTN(LP1,LP2,LP3)
  25.       COMMON /COEF/ AW(LP1,LP2,LP3),AE(LP1,LP2,LP3),AN(LP1,LP2,LP3),
  26.      & AS(LP1,LP2,LP3),AB(LP1,LP2,LP3),AT(LP1,LP2,LP3),E(LP1,LP2,LP3),
  27.      & B(LP1,LP2,LP3)
  28.       COMMON /RUNSUM/ ITSNO(LP12),STIME(LP12),SOPROD(LP12),
  29.      & SGPROD(LP12),
  30.      & SWPROD(LP12),SGOR(LP12),SWOR(LP12),SGINJ(LP12),SWINJ(LP12)
  31.       COMMON /RUN2/SPVWTP(LP12),SOCUMP(LP12),SWCUMP(LP12),SGCUMP(LP12),
  32.      & SGCUMI(LP12),SWCUMI(LP12),SAQUIR(LP12),SAQUIC(LP12)
  33.       COMMON /SAQUI/ IAQOPT,CPIAQ1(LP1,LP2,LP3),CPIAQ2(LP1,LP2,LP3),
  34.      & CPI1(LP1,LP2,LP3),CPI2(LP1,LP2,LP3),EWAQ(LP1,LP2,LP3),
  35.      & CUMAQW(LP1,LP2,LP3),
  36.      & QWAQ(LP1,LP2,LP3),CUMEW(LP1,LP2,LP3),QWAQR(LP7),CUMAQR(LP7)
  37.      & ,IAQREG(LP1,LP2,LP3),PAQ(LP1,LP2,LP3),PIAQ(LP1,LP2,LP3)
  38.       COMMON /SARRAY/ PN(LP1,LP2,LP3),IOCODE,IDMAX,
  39.      & SON(LP1,LP2,LP3),SWN(LP1,LP2,LP3),SGN(LP1,LP2,LP3),
  40.      & A1(LP1,LP2,LP3),A2(LP1,LP2,LP3),A3(LP1,LP2,LP3),
  41.      & SUM(LP1,LP2,LP3),GAM(LP1,LP2,LP3),QS(LP1,LP2,LP3)
  42.       COMMON /SLIMIT/ GORT(LP11),WORT(LP11),ILIMOP(LP11),
  43.      & GORL(LP11),WORL(LP11),QOC(LP11,LP3),QWC(LP11,LP3),QGC(LP11,LP3)
  44.       COMMON /SPARM/ KX(LP1,LP2,LP3),KY(LP1,LP2,LP3),KZ(LP1,LP2,LP3),
  45.      & EL(LP1,LP2,LP3),TX(LP4,LP2,LP3),TY(LP1,LP5,LP3),TZ(LP1,LP2,LP6),
  46.      & PDAT(LP1,LP2,LP3),PDATUM,GRAD
  47.       COMMON /SPOST/ KOPR,KGPR,KWPR,KGOR,KWOR,KGIR,KWIR,KRESP,
  48.      & KAIR,KAIC,KCOP,KCGP,KCWP,KCGI,KCWI,ITSMAX
  49.       COMMON /SPRTPS/ P(LP1,LP2,LP3),SO(LP1,LP2,LP3),SW(LP1,LP2,LP3),
  50.      & SG(LP1,LP2,LP3)
  51.       COMMON /SPVT/ SAT(LP7,LP9),KROT(LP7,LP9),KRWT(LP7,LP9),
  52.      & BGT(LP7,LP9),
  53.      & KRGT(LP7,LP9),ITHREE(LP7),RSOT(LP7,LP9),BWPT(LP7,LP9),
  54.      & PCOWT(LP7,LP9),PCGOT(LP7,LP9),KROGT(LP7,LP9),SWR(LP7),
  55.      & POT(LP7,LP9),MUOT(LP7,LP9),BOT(LP7,LP9),BOPT(LP7,LP9),
  56.      & RSOPT(LP7,LP9),PWT(LP7,LP9),MUWT(LP7,LP9),BWT(LP7,LP9),
  57.      & RSWT(LP7,LP9),RSWPT(LP7,LP9),PGT(LP7,LP9),MUGT(LP7,LP9),
  58.      & BGPT(LP7,LP9),CRT(LP7,LP9),IPVT(LP1,LP2,LP3),IROCK(LP1,LP2,LP3),
  59.      & NROCK,NPVT,PSIT(LP7,LP9),PRT(LP7,LP9),WOROCK(LP7),GOROCK(LP7)
  60.       COMMON /SRATE/ PID(LP11,LP3),PWF(LP11,LP3),PWFC(LP11,LP3),
  61.      & KIP(LP11),LAYER(LP11),QVO(LP11),CUMG(LP11,LP3),
  62.      & GMO(LP11,LP3),GMW(LP11,LP3),GMG(LP11,LP3),
  63.      & QVW(LP11),QVG(LP11),QVT(LP11),CUMO(LP11,LP3),CUMW(LP11,LP3),
  64.      & IDWELL(LP11),ALIT(LP11),BLIT(LP11)
  65.       COMMON /SSOLN/ BO(LP1,LP2,LP3),BW(LP1,LP2,LP3),BG(LP1,LP2,LP3),
  66.      & QO(LP1,LP2,LP3),QW(LP1,LP2,LP3),QG(LP1,LP2,LP3),
  67.      & GOWT(LP1,LP2,LP3),GWWT(LP1,LP2,LP3),GGWT(LP1,LP2,LP3),
  68.      & OW(LP4,LP2,LP3),OE(LP4,LP2,LP3),WW(LP4,LP2,LP3),WE(LP4,LP2,LP3),
  69.      & OS(LP1,LP5,LP3),ON(LP1,LP5,LP3),WS(LP1,LP5,LP3),WN(LP1,LP5,LP3),
  70.      & OT(LP1,LP2,LP6),OB(LP1,LP2,LP6),WT(LP1,LP2,LP6),WB(LP1,LP2,LP6),
  71.      & QOWG(LP1,LP2,LP3),VP(LP1,LP2,LP3),CT(LP1,LP2,LP3)
  72.       COMMON /TRIDI/ UM(LP15),AZL(LP15),BZL(LP15),CZL(LP15),DZL(LP15),
  73.      & UZL(LP15),X(LP15),BETA(LP15),GAMMA(LP15),W(LP15)
  74.       COMMON /TSTDAT/ IFATAL,IWARN
  75.       COMMON /VECTOR/ DX(LP1,LP2,LP3),DY(LP1,LP2,LP3),DZ(LP1,LP2,LP3),
  76.      & DZNET(LP1,LP2,LP3),IQN1(LP11),IQN2(LP11),IQN3(LP11),IHEDIN(80)
  77.       COMMON /CHAR/ WELNAM(LP11)
  78. C
  79.       DIMENSION FTIO(50),OOIP(LP3),OWIP(LP3),ODGIP(LP3),OFGIP(LP3)
  80.       DIMENSION IRETYM(LP17),REDATE(LP17)
  81.       DIMENSION RKCOEF(7,7),THRUA(LP1,LP2,LP3)
  82.       DIMENSION WOC(LP7),GOC(LP7)
  83. C
  84.       DATA KRUNGE/1/
  85.       DATA RKCOEF
  86.      & /1.0,6*0.0,
  87.      & 0.13157894,1.0,5*0.0,
  88.      & 0.03937504,0.15625736,1.0,4*0.0,
  89.      & 0.01666967,0.053199139,0.16491829,1.0,3*0.0,
  90.      & 0.0085487791,0.024394839,0.059605275,0.16893119,1.0,2*0.0,
  91.      & .0049515953,.013195402,.028593749,.06308716,.17111215,1.0,0.0,
  92.      & .0031198738,.0079382925,.015998138,.031126418,.06518731,
  93.      & .17242757,1.0/
  94.       DATA ETI,FT,FTMAX/0.0,0.0,0.0/
  95.       DATA COP,CWP,CGP,CWI,CGI/0.0,0.0,0.0,0.0,0.0/
  96.       DATA IOMETH,IFTCOD,IRSCNT/0,0,1/
  97. C
  98.       DATA NX,NY,NZ,NTER,NTEP,NW,
  99.      & NTMAX,NRST,N1DIR,N23DIR,NITER,NTRI
  100.      & /LP1,LP2,LP3,LP9,LP9,LP11,LP12,LP17,LP15,
  101.      & LP14,LP22,LP15/
  102.       IDMAX=0.0
  103.       CMBEO=0.0
  104.       CMBEW=0.0
  105.       CMBEG=0.0
  106. C      OPEN FILES FOR I/O ON DEC 20/60
  107. C      MACHINE DEPENDENT OPEN STATEMENTS
  108.       WRITE(*,3)
  109.     3 FORMAT(1X,'ENTER INPUT DATA FILE NAME ... ')
  110.       READ(*,7) INDAT
  111.       OPEN(UNIT=20,MODE='READ',
  112.      & FILE=INDAT,STATUS='OLD')
  113.       WRITE(*,6)
  114.     6 FORMAT(/,1X,'ENTER OUTPUT DATA FILE NAME ... ')
  115.       READ(*,7) OUTDAT
  116.     7 FORMAT(2A10)
  117.       OPEN(UNIT=21,MODE='WRITE',STATUS='NEW',
  118.      & FILE=OUTDAT)
  119.       IOCODE=21
  120. C      WRITE BANNER
  121.       WRITE(IOCODE,3005)
  122.       WRITE(IOCODE,3007)
  123.       WRITE(IOCODE,3006)
  124.  3007 FORMAT(T30,'*',T50,'          BOAST II:',T98,'*',/T30,'*',
  125.      & T48,'BLACK OIL APPLIED SIMULATION TOOL',T98,'*',
  126.      & /T30,'*',T57,'(RELEASE 1.1)',T98,'*')
  127.  3005 FORMAT(//T30,69('*'),/T30,'*',T98,'*'/T30,'*',T98,'*')
  128.  3006 FORMAT(T30,'*',T98,'*'/T30,'*',T98,'*'/T30,69('*'),///)
  129.       READ(20,69) (IHEDIN(IH),IH=1,40)
  130.       WRITE(IOCODE,68) (IHEDIN(IH),IH=1,40)
  131. 69    FORMAT(40A2)
  132.    68 FORMAT(///,T24,82('*'),/,T24,('*'),T105,('*'),/,
  133.      & T24,('*'),40A2,T105,('*'),/,T24,('*'),T105,
  134.      & ('*'),/,T24,82('*'),//)
  135. C      WRITE DIMENSION INFO
  136.       WRITE(IOCODE,3100) NX,NY,NZ,NROCK,NTER,NPVT,NTEP,NW,
  137.      & NTMAX,NRST,N1DIR,N23DIR,NITER,NTRI,N23DIR
  138.  3100 FORMAT(T37,'REDIMENSIONING INFORMATION:',
  139.      & /T32,'MAX X-DIRECTION GRID BLOCKS',T92,I5,
  140.      & /T32,'MAX Y-DIRECTION GRID BLOCKS',T92,I5,
  141.      & /T32,'MAX Z-DIRECTION GRID BLOCKS',T92,I5,
  142.      & /T32,'MAX ROCK REGIONS',T92,I5,
  143.      & /T32,'MAX ROCK REGION TABLE ENTRIES',T92,I5,
  144.      & /T32,'MAX PVT REGIONS',T92,I5,
  145.      & /T32,'MAX PVT REGION TABLE ENTRIES',T92,I5,
  146.      & /T32,'MAX WELLS',T92,I5,
  147.      & /T32,'MAX TIME STEPS',T92,I5,
  148.      & /T32,'MAX RESTART RECORDS',T92,I5,
  149.      & /T32,'TOTAL BLOCKS USING 1D DIRECT SOLN METHODS',T92,I5,
  150.      & /T32,'TOTAL BLOCKS USING 2D OR 3D DIRECT SOLN METHODS',T92,I5,
  151.      & /T32,'TOTAL BLOCKS USING ITERATIVE SOLN METHOD',T92,I5,
  152.      & /T32,'MAX NO OF BLOCKS IN 1D FOR LSOR',T92,I5,
  153.      & /T32,'MAX NO OF BLOCKS IN 2D FOR L2SOR',T92,I5)
  154.       WRITE(IOCODE,70) (IHEDIN(IH),IH=1,40)
  155.    70 FORMAT(1H1/T24,82('*'),/,T24,('*'),T105,('*'),/,
  156.      & T24,('*'),40A2,T105,('*'),/,T24,('*'),T105,
  157.      & ('*'),/,T24,82('*'),//)
  158. C** RESTART OPTION READ
  159.       READ(20,69)
  160.       READ(20,*) IREOPT,IPLOTP
  161.       IF(IREOPT.GT.-1) READ(20,*) IRNUM,IRSTRT,NN,TMAX
  162.       IF(IREOPT.EQ.-1) WRITE(IOCODE,12)
  163.       IF(IREOPT.EQ.0) WRITE(IOCODE,13)
  164.       IF(IREOPT.EQ.1) WRITE(IOCODE,14) IRSTRT
  165.       IF(IREOPT.GT.-1)
  166.      & WRITE(IOCODE,11) IREOPT,IRNUM,IRSTRT,NN,TMAX
  167.    11 FORMAT(/30X,'RESTART CODES:',
  168.      & /25X,'RESTART OPTION (IREOPT)  ',T70,I5,
  169.      & /25X,'NUMBER OF RESTART ENTRIES (IRNUM)  ',T70,I5,
  170.      & /25X,'RESTART TIME STEP (IRSTRT)  ',T70,I5,
  171.      & /25X,'MAXIMUM NUMBER OF TIME STEPS (NN)  ',T70,I5,
  172.      & /25X,'MAXIMUM NUMBER OF DAYS (TMAX)',T70,F10.1/)
  173.    12 FORMAT(/T52,5('*'),' RESTART OPTION ',5('*'),
  174.      & //T45,'RESTART OPTION HAS NOT BEEN ACTIVATED.',//)
  175.    13 FORMAT(/T52,5('*'),' RESTART OPTION ',5('*'),
  176.      & //T45,'RESTART OPTION HAS BEEN ACTIVATED.',//)
  177.    14 FORMAT(/T52,5('*'),' RESTART OPTION ',5('*'),
  178.      & //T41,'THIS IS A RESTART RUN BEGINNING AT TIME STEP',
  179.      & I5,//)
  180.       IF(IREOPT.LT.0) GO TO 18
  181.       READ(20,7) RESIN,RESOUT
  182.       WRITE(IOCODE,3003) RESIN,RESOUT
  183.  3003 FORMAT(/25X,'INPUT RESTART FILE    ',T70,A10,
  184.      & /,25X,'INPUT RESTART FILE    ',T70,A10,/)
  185. C      MACHINE DEPENDENT OPEN STATEMENTS
  186.       OPEN(UNIT=24,FORM='BINARY',ACCESS='SEQUENTIAL',
  187.      & FILE=RESIN,STATUS='UNKNOWN')
  188.       OPEN(UNIT=25,FORM='BINARY',ACCESS='SEQUENTIAL',
  189.      & FILE=RESOUT,STATUS='UNKNOWN')
  190.       READ(20,*) (IRETYM(I),I=1,IRNUM)
  191.       READ(20,*) (REDATE(I),I=1,IRNUM)
  192.       WRITE(IOCODE,15) (IRETYM(I),I=1,IRNUM)
  193.    15 FORMAT(25X,'RESTART TIME STEPS:',5(6X,I5))
  194.       WRITE(IOCODE,16) (REDATE(I),I=1,IRNUM)
  195.    16 FORMAT(25X,'RESTART TIMES (DAYS):',5(1X,F10.2))
  196.       WRITE(IOCODE,33)
  197.   33  FORMAT(//)
  198.    18 CONTINUE
  199. C**** POST-PLOT PACKAGE CODES
  200.       IF(IPLOTP.EQ.1) READ(20,*) NPLINE,
  201.      & KOPR,KGPR,KWPR,KGOR,KWOR,
  202.      & KGIR,KWIR,KRESP,KAIR,KAIC,
  203.      & KCOP,KCGP,KCWP,KCGI,KCWI
  204.       IF(IPLOTP.GE.0) WRITE(IOCODE,3015)
  205.  3015 FORMAT(/T50,5('*'),' POST-PLOT OPTIONS ',5('*'),
  206.      & //T40,'POST-PLOT TABLE WILL BE WRITTEN.')
  207.       IF(IPLOTP.EQ.1) WRITE(IOCODE,3020) NPLINE,
  208.      & KOPR,KGPR,KWPR,KGOR,KWOR,
  209.      & KGIR,KWIR,KRESP,KAIR,KAIC,
  210.      & KCOP,KCGP,KCWP,KCGI,KCWI
  211.  3020 FORMAT(/T40,'PLOT LINES/TIME STEP = ',I5,
  212.      & /T40,'OIL PROD RATE CODE   = ',I5,
  213.      & /T40,'GAS PROD RATE CODE   = ',I5,
  214.      & /T40,'WATER PROD RATE CODE = ',I5,
  215.      & /T40,'PROD G/O RATIO CODE  = ',I5,
  216.      & /T40,'PROD W/O RATIO CODE  = ',I5,
  217.      & /T40,'GAS INJ RATE CODE    = ',I5,
  218.      & /T40,'WATER INJ RATE CODE  = ',I5,
  219.      & /T40,'PV WT AVG RES P CODE = ',I5,
  220.      & /T40,'AQ INFLUX RATE CODE  = ',I5,
  221.      & /T40,'AQ INFLUX CUM  CODE  = ',I5,
  222.      & /T40,'CUM OIL PROD CODE    = ',I5,
  223.      & /T40,'CUM GAS PROD CODE    = ',I5,
  224.      & /T40,'CUM WATER PROD CODE  = ',I5,
  225.      & /T40,'CUM GAS INJ CODE     = ',I5,
  226.      & /T40,'CUM WATER INJ CODE   = ',I5,//)
  227.       IF(IREOPT.GT.0) GO TO 20
  228. C**** ESTABLISH RESERVOIR AND BLOCK DIMENSIONS
  229.       CALL GRIDSZ(IOCODE,II,JJ,KK)
  230. C      GRID DIMENSION LIMIT CHECK
  231.       IF(II.LE.NX) GO TO 1520
  232.       IFATAL=IFATAL+1
  233.       WRITE(IOCODE,1510) NX
  234.  1510 FORMAT(/5X,5('-'),'NO. OF X-DIRECTION BLOCKS CANNOT',
  235.      & ' EXCEED ',I5)
  236.  1520 IF(JJ.LE.NY) GO TO 1540
  237.       IFATAL=IFATAL+1
  238.       WRITE(IOCODE,1530) NY
  239.  1530 FORMAT(/5X,5('-'),'NO. OF Y-DIRECTION BLOCKS CANNOT',
  240.      & ' EXCEED ',I5)
  241.  1540 IF(KK.LE.NZ) GO TO 1560
  242.       IFATAL=IFATAL+1
  243.       WRITE(IOCODE,1550) NZ
  244.  1550 FORMAT(/5X,5('-'),'NO. OF Z-DIRECTION BLOCKS CANNOT',
  245.      & ' EXCEED ',I5)
  246.  1560 CONTINUE
  247. C**** ESTABLISH POROSITY AND PERMEABILITY AT EACH ZONE
  248.       CALL PORPRM(IOCODE,II,JJ,KK)
  249. C**** CALCULATE INTERBLOCK TRANSMISSIBILITIES
  250.       CALL TRANS(II,JJ,KK)
  251. C**** EMPIRICAL DATA
  252.       CALL TABLE(IOCODE,II,JJ,KK)
  253. C**** ESTABLISH INITIAL CONDITIONS
  254.       CALL UINITL(KPI,II,JJ,KK,PI,ET0,CUMPO,MBEO,
  255.      &CUMPW,MBEW,CUMPG,MBEG,SOI,SWI,SGI,
  256.      &WOC,GOC,PGOC,CUMIW,CUMIG)
  257. C**** SOLUTION METHOD, DEBUG PRINT, AND TIME-STEP CONTROL PARAMETERS.
  258.       CALL CODES(IOCODE,KSM1,KSN1,KCO1,NN,
  259.      &FACT1,FACT2,TMAX,KSOL,MITER,OMEGA,TOL,TOL1,
  260.      & KSN,KSM,KCO,KCOFF,DSMAX,DPMAX,
  261.      & NUMDIS,IRK,THRUIN,
  262.      &WORMAX,GORMAX,PAMIN,PAMAX)
  263. C** AQUIFER DATA ENTRY
  264.       CALL AQUI(TMAX)
  265. C**** IF A FATAL INPUT ERROR OCCURS, DO NOT BEGIN RUN
  266.        IF(IFATAL.GE.1) GO TO 1006
  267.       IF(IWARN.GT.0) WRITE(IOCODE,19) IWARN
  268.    19 FORMAT(/5X,5('-'),'A TOTAL OF',I5,
  269.      & ' WARNINGS HAVE BEEN NOTED',5('-'),/)
  270.    20 READ(20,69)
  271.       D5615=1.0/5.615
  272.       D288=1.0/288.
  273.       D144=1.0/144.
  274. C      DETERMINE TIME STEP INDEX RANGE
  275.       NINIT=1
  276.       IF(IREOPT.NE.1) GO TO 24
  277.       NINIT=IRSTRT
  278.    24 NMAX=NN+1
  279. C      RECURRENT CALC LOOP
  280.       DO 1000 N=NINIT,NMAX
  281.       NLOOP=N
  282. C** RESTART ALGORITHM
  283.       IF(IREOPT.EQ.-1) GO TO 30
  284.       IF(IREOPT.EQ.0) GO TO 25
  285.       IF(N.LT.IRSTRT) GO TO 1000
  286.       IF(N.GT.IRSTRT) GO TO 25
  287.       IRS=0
  288.    22 IRS=IRS+1
  289.       READ (24) CPIAQ1,CPIAQ2,CPI1,CPI2,EWAQ,CUMAQW,
  290.      & QWAQ,CUMEW,IAQREG,PAQ,PIAQ,
  291.      & P,SO,SW,SG,PN,SON,SWN,SGN,
  292.      & CUMO,CUMW,CUMG,TX,TY,TZ,
  293.      & QOC,QWC,QGC,PID,PWF,PWFC,GMO,GMW,GMG,
  294.      & BO,BW,BG,VP,CT,EL,PBOT,PBOTN
  295. C**    VECTORS
  296.       READ (24) QWAQR,CUMAQR,
  297.      & ITSNO,STIME,SOPROD,SGPROD,SWPROD,SGOR,SWOR,
  298.      & SGINJ,SWINJ,SPVWTP,SOCUMP,SWCUMP,SGCUMP,
  299.      & SGCUMI,SWCUMI,SAQUIR,SAQUIC,
  300.      & GORT,WORT,GORL,WORL,
  301.      & ILIMOP,KIP,WELNAM,LAYER,IDWELL,
  302.      & QVO,QVG,QVW,QVT,
  303.      & IQN1,IQN2,IQN3,
  304.      & MSAT,MPOT,MPWT,MPGT,RHOSCO,RHOSCW,RHOSCG,
  305.      & VSLOPE,BSLOPE,RSLOPE,
  306.      & SAT,KROT,KRWT,KRGT,KROGT,PCOWT,PCGOT,
  307.      & POT,MUOT,BOT,BOPT,RSOT,RSOPT,
  308.      & PWT,MUWT,BWT,BWPT,RSWT,RSWPT,
  309.      & PGT,PSIT,MUGT,BGT,BGPT,PRT,CRT
  310. C**    SCALARS
  311.       READ (24) II,JJ,KK,KSM1,KSN1,KCO1,KSOL,MITER,
  312.      & NUMDIS,IRK,THRUIN,
  313.      & FACT1,FACT2,OMEGA,TOL,TOL1,WORMAX,GORMAX,PAMIN,PAMAX,
  314.      & NVQN,CWI,CGI,ETI,PBO,PMAXT,
  315.      & STBO,STBOI,STBW,STBWI,MCFG,MCFGI,MBEO,MBEW,MBEG,
  316.      & CMBEO,CMBEW,CMBEG,TOOIP,TOWIP,TOGIP,TODGIP,TOFGIP,
  317.      & DELT0,RESVOL,OP,WP,GP,WI,GI,PAVG0,PAVG,OPR,WPR,
  318.      & GPR,WIR,GIR,COP,CWP,CGP,MCFG1,MCFGT,CWOR,WOR,
  319.      & CGOR,GOR,DSMAX,DPMAX,KCOFF,KSN,KSM,KCO,
  320.      & SWR,PDATUM,GRAD,DSMC,DPMC,IAQOPT,
  321.      & NTSTEP,IREPRS,ITHREE,
  322.      & NROCK,NPVT,IROCK,IPVT
  323.       WRITE(IOCODE,31) N,NTSTEP,ETI
  324.    31 FORMAT(/T10,110('-'),/T10,
  325.      & 'RESTART DATA READ AT BEGINNING OF TIME STEP ',I5,
  326.      & '; NTSTEP =',I5,' AND ELAPSED TIME (ETI) =',F10.2,
  327.      & ' DAYS.',/T10,110('-'),/)
  328.       IF(N.EQ.NTSTEP) GO TO 25
  329.       IF(IRS.LT.5) GO TO 22
  330.       WRITE(IOCODE,23) IRS,N,ETI,NTSTEP
  331.    23 FORMAT(/T10,110('-'),/T10,
  332.      & 'ATTEMPT TO READ RESTART DATA FAILED AFTER',I3,
  333.      & ' TRIES AT TIME STEP',I5,' (TIME =',F10.2,' DAYS).',
  334.      & /T10,'LAST RESTART TIME STEP READ =',I5,
  335.      & /T10,110('-'),/)
  336.       GO TO 1006
  337.    25 CONTINUE
  338.       IF(IRETYM(IRSCNT).GT.0.AND.N.NE.IRETYM(IRSCNT)) GO TO 30
  339.       IF(REDATE(IRSCNT).GT.0.AND.ETI.NE.REDATE(IRSCNT)) GO TO 30
  340.       IF(IRETYM(IRSCNT).EQ.0.AND.REDATE(IRSCNT).EQ.0.0) GO TO 30
  341.       IF(IRSCNT.GT.5) GO TO 28
  342.       NTSTEP=N
  343. C**    ARRAYS
  344.       WRITE (25) CPIAQ1,CPIAQ2,CPI1,CPI2,EWAQ,CUMAQW,
  345.      & QWAQ,CUMEW,IAQREG,PAQ,PIAQ,
  346.      & P,SO,SW,SG,PN,SON,SWN,SGN,
  347.      & CUMO,CUMW,CUMG,TX,TY,TZ,
  348.      & QOC,QWC,QGC,PID,PWF,PWFC,GMO,GMW,GMG,
  349.      & BO,BW,BG,VP,CT,EL,PBOT,PBOTN
  350. C**    VECTORS
  351.       WRITE (25) QWAQR,CUMAQR,
  352.      & ITSNO,STIME,SOPROD,SGPROD,SWPROD,SGOR,SWOR,
  353.      & SGINJ,SWINJ,SPVWTP,SOCUMP,SWCUMP,SGCUMP,
  354.      & SGCUMI,SWCUMI,SAQUIR,SAQUIC,
  355.      & GORT,WORT,GORL,WORL,
  356.      & ILIMOP,KIP,WELNAM,LAYER,IDWELL,
  357.      & QVO,QVG,QVW,QVT,
  358.      & IQN1,IQN2,IQN3,
  359.      & MSAT,MPOT,MPWT,MPGT,RHOSCO,RHOSCW,RHOSCG,
  360.      & VSLOPE,BSLOPE,RSLOPE,
  361.      & SAT,KROT,KRWT,KRGT,KROGT,PCOWT,PCGOT,
  362.      & POT,MUOT,BOT,BOPT,RSOT,RSOPT,
  363.      & PWT,MUWT,BWT,BWPT,RSWT,RSWPT,
  364.      & PGT,PSIT,MUGT,BGT,BGPT,PRT,CRT
  365. C**    SCALARS
  366.       WRITE (25) II,JJ,KK,KSM1,KSN1,KCO1,KSOL,MITER,
  367.      & NUMDIS,IRK,THRUIN,
  368.      & FACT1,FACT2,OMEGA,TOL,TOL1,WORMAX,GORMAX,PAMIN,PAMAX,
  369.      & NVQN,CWI,CGI,ETI,PBO,PMAXT,
  370.      & STBO,STBOI,STBW,STBWI,MCFG,MCFGI,MBEO,MBEW,MBEG,
  371.      & CMBEO,CMBEW,CMBEG,TOOIP,TOWIP,TOGIP,TODGIP,TOFGIP,
  372.      & DELT0,RESVOL,OP,WP,GP,WI,GI,PAVG0,PAVG,OPR,WPR,
  373.      & GPR,WIR,GIR,COP,CWP,CGP,MCFG1,MCFGT,CWOR,WOR,
  374.      & CGOR,GOR,DSMAX,DPMAX,KCOFF,KSN,KSM,KCO,
  375.      & SWR,PDATUM,GRAD,DSMC,DPMC,IAQOPT,
  376.      & NTSTEP,IREPRS,ITHREE,
  377.      & NROCK,NPVT,IROCK,IPVT
  378.       WRITE(IOCODE,26) NTSTEP,ETI
  379.    26 FORMAT(/T10,100('-'),/T10,'RESTART RECORD WRITTEN',
  380.      & ' AT BEGINNING OF TIME STEP',I5,' (TIME =',F10.2,' DAYS)',
  381.      & /T10,100('-'),//)
  382.       IRSCNT=IRSCNT+1
  383.       GO TO 30
  384.    28 WRITE(IOCODE,29) IRSCNT,NTSTEP,ETI
  385.    29 FORMAT(/T10,100('-'),/T10,
  386.      & 'THE NUMBER OF RESTART WRITE (',I2,') EXCEEDS THE',
  387.      & ' ALLOWED NUMBER AT TIME STEP',I5,
  388.      & ' (TIME =',F10.2,' DAYS)',/T10,100('-'),/)
  389.       GO TO 1006
  390.    30 CONTINUE
  391. C      END OF RESTART
  392. C**** RECURRENT DATA.
  393.       IF(N.EQ.1) GO TO 40
  394.       IF(FT.LT.FTMAX) GO TO 46
  395.       IF(IOMETH.EQ.0) GO TO 40
  396.       IF(IFTCOD.EQ.IOMETH) GO TO 40
  397.       IFTCOD=IFTCOD+1
  398.       FTMAX=FTIO(IFTCOD)
  399.       IF(FTMAX.GT.TMAX) FTMAX=TMAX
  400.       GO TO 46
  401.    40 READ(20,69,END=1006)
  402.       READ(20,*) ICHANG,IWLCNG,IOMETH
  403.       IF(IOMETH.NE.0) GO TO 42
  404.       READ(20,*) IWLREP,ISUMRY
  405.       READ(20,*) IPMAP,ISOMAP,ISWMAP,ISGMAP,IPBMAP,IAQMAP
  406.       READ(20,*) DAY,DTMIN,DTMAX
  407.       DELT=DAY
  408.       FTMAX=ETI+ICHANG*DELT
  409.       IF(FTMAX.GT.TMAX) FTMAX=TMAX
  410.       IFREQW=0
  411.       IFREQS=0
  412.       GO TO 43
  413.    42 READ(20,*) (FTIO(I),I=1,IOMETH)
  414.       READ(20,*) IPMAP,ISOMAP,ISWMAP,ISGMAP,IPBMAP,IAQMAP
  415.       READ(20,*) DAY,DTMIN,DTMAX
  416.       DELT=DAY
  417.       IFTCOD=1
  418.       FTMAX=FTIO(IFTCOD)
  419.       IF(FTMAX.GT.TMAX) FTMAX=TMAX
  420.    43 CONTINUE
  421.       IF(IWLCNG.EQ.0) GO TO 45
  422.       CALL NODES(NVQN)
  423. C      TIME STEP SUMMARY HEADING ACTIVATED
  424.       ITSS=0
  425.    45 CONTINUE
  426.    46 IF(N.EQ.1)DELT0=DELT
  427.       IF(N.EQ.1)GO TO 1050
  428. C      AUTO TIME STEP ADJUSTMENT
  429.       IF(DSMC.LT.DSMAX.AND.DPMC.LT.DPMAX) DELT=DELT*FACT1
  430.       IF(DSMC.GT.DSMAX.OR.DPMC.GT.DPMAX) DELT=DELT*FACT2
  431.       IF(DELT.GT.DTMAX)DELT=DTMAX
  432.       IF(IRK.EQ.0) GO TO 1045
  433. C      R-K SOLN LOOP MAX.
  434.       DO 47 I=1,II
  435.       DO 47 J=1,JJ
  436.       DO 47 K=1,KK
  437.    47 THRUA(I,J,K)=0.0
  438.       DO 48 J=1,NVQN
  439.       IQ1=IQN1(J)
  440.       IQ2=IQN2(J)
  441.       IQ3=IQN3(J)
  442.       IJLOC=IDWELL(J)
  443. C MAY NEED TO CHANGE THIS BACK TO 481?
  444.       IF(IJLOC.EQ.0) GO TO 48
  445.       LAY=IQ3+(LAYER(J)-1)
  446.       DO 48 K=IQ3,LAY
  447.       BPT=PBOT(IQ1,IQ2,K)
  448.       PP=P(IQ1,IQ2,K)
  449.       IPVTR=IPVT(IQ1,IQ2,K)
  450.       CALL INTPVT(IPVTR,BPT,BSLOPE(IPVTR),POT,BOT,MPOT(IPVTR),PP,BBO)
  451.       CALL INTERP(IPVTR,PWT,BWT,MPWT(IPVTR),PP,BBW)
  452.       CALL INTERP(IPVTR,PGT,BGT,MPGT(IPVTR),PP,BBG)
  453.       THRUA(IQ1,IQ2,K)=THRUA(IQ1,IQ2,K) + DELT*
  454.      & (QOC(IJLOC,K)*BBO+QWC(IJLOC,K)*BBW+QGC(IJLOC,K)*BBG)
  455.    48 CONTINUE
  456.       THRU=0.
  457.       DO 49 J=1,NVQN
  458.       IQ1=IQN1(J)
  459.       IQ2=IQN2(J)
  460.       IQ3=IQN3(J)
  461.       LAY=IQ3+(LAYER(J)-1)
  462.       DO 49 K=IQ3,LAY
  463.       IF(VP(IQ1,IQ2,K).LE.0.0) GO TO 49
  464.       THRUPV=ABS(THRUA(IQ1,IQ2,K))/VP(IQ1,IQ2,K)
  465.       IF(THRU.LT.THRUPV) THRU=THRUPV
  466.    49 CONTINUE
  467.       KRUNGE=1.+SQRT(THRU/THRUIN)
  468.       IF(KRUNGE.LE.7) GO TO 1045
  469.       DELT=DELT*7./KRUNGE
  470.       KRUNGE=7
  471.  1045 CONTINUE
  472.       IF(DELT.LT.DTMIN)DELT=DTMIN
  473.       IF(ETI+DELT.GT.FTMAX)DELT=FTMAX-ETI
  474. 1050  FT=ETI+DELT
  475.       IF(ETI+DELT*0.5.GE.TMAX)GO TO 1006
  476. C****ITFLAG COUNTS THE NUMBER OF TIME STEP REPETITIONS.
  477.       ITFLAG=0
  478. C****REENTRY POINT FOR REPEATED TIME STEP.
  479.  1060 CONTINUE
  480.       DIV1=1.0/DELT
  481.       IF(N.GT.1.OR.ITFLAG.GT.0)GO TO 105
  482. C**  INITIAL FLUID CALCULATIONS
  483.       WRITE(IOCODE,80)
  484.    80 FORMAT(1H1/T15,'******* INITIAL FLUID VOLUMES *******',/)
  485.       RESVOL=0.0
  486.       SCFO=0.0
  487.       SCFG=0.0
  488.       SCFG1=0.0
  489.       DO 102 K=1,KK
  490.       OOIP(K)=0.
  491.       OWIP(K)=0.
  492.       ODGIP(K)=0.
  493.       OFGIP(K)=0.
  494.       DO 100 J=1,JJ
  495.       DO 100 I=1,II
  496.       PP=P(I,J,K)
  497.       BPT=PBOT(I,J,K)
  498.       RESVOL=RESVOL+VP(I,J,K)
  499.       IPVTR=IPVT(I,J,K)
  500. C**** NOTE: WE ARE ASSUMING INITIAL PHI IS AT INITIAL RESERVOIR PRESSURE
  501.       CALL INTPVT(IPVTR,BPT,RSLOPE(IPVTR),POT,RSOT,MPOT(IPVTR),PP,RSO)
  502.       CALL INTERP(IPVTR,PWT,RSWT,MPWT(IPVTR),PP,RSW)
  503.       CALL INTPVT(IPVTR,BPT,BSLOPE(IPVTR),POT,BOT,MPOT(IPVTR),
  504.      & PP,BO(I,J,K))
  505.       CALL INTERP(IPVTR,PWT,BWT,MPWT(IPVTR),PP,BW(I,J,K))
  506.       CALL INTERP(IPVTR,PGT,BGT,MPGT(IPVTR),PP,BG(I,J,K))
  507.       FF1=SO(I,J,K)/BO(I,J,K)
  508.       FF2=SW(I,J,K)/BW(I,J,K)
  509.       SCFO=SCFO+VP(I,J,K)*FF1
  510.       SCFW=SCFW+VP(I,J,K)*FF2
  511.       SCFG=SCFG+VP(I,J,K)*SG(I,J,K)/BG(I,J,K)
  512.       SCFG1=SCFG1+VP(I,J,K)*(RSO*FF1+RSW*FF2)
  513. C      CALC COMPRESSIBILITY
  514.       CALL INTCOM(IPVTR,POT,BOPT,MPOT(IPVTR),PP,BODER)
  515.       CALL INTCOM(IPVTR,POT,RSOPT,MPOT(IPVTR),PP,RSODER)
  516.       CALL INTCOM(IPVTR,PWT,BWPT,MPWT(IPVTR),PP,BWDER)
  517.       CALL INTCOM(IPVTR,PWT,RSWPT,MPWT(IPVTR),PP,RSWDER)
  518.       CALL INTCOM(IPVTR,PGT,BGPT,MPGT(IPVTR),PP,BGDER)
  519.       IF(PP.GT.PBOT(I,J,K))BODER=BSLOPE(IPVTR)
  520.       IF(PP.GT.PBOT(I,J,K))RSODER=RSLOPE(IPVTR)
  521.       CO=-(BODER-BG(I,J,K)*RSODER)/BO(I,J,K)
  522.       CW=-(BWDER-BG(I,J,K)*RSWDER)/BW(I,J,K)
  523.       CG=-BGDER/BG(I,J,K)
  524.       CALL INTERP(IPVTR,PRT,CRT,MPGT(IPVTR),PP,CR)
  525.       CT(I,J,K)=CR + CO*SO(I,J,K) +CW*SW(I,J,K) + CG*SG(I,J,K)
  526. C      CALC INITIAL FLUIDS IN PLACE
  527.       OOIP(K)=OOIP(K)+D5615*0.000001*SON(I,J,K)*VP(I,J,K)/BO(I,J,K)
  528.       OWIP(K)=OWIP(K)+D5615*0.000001*SWN(I,J,K)*VP(I,J,K)/BW(I,J,K)
  529.       ODGIP(K)=ODGIP(K)+0.001*0.000001*(RSO*SON(I,J,K)*
  530.      & VP(I,J,K)/BO(I,J,K)+RSW*SWN(I,J,K)*VP(I,J,K)/BW(I,J,K))
  531.       OFGIP(K)=OFGIP(K)+0.001*0.000001*SGN(I,J,K)*VP(I,J,K)/BG(I,J,K)
  532.       IF(KCO1.EQ.0)GO TO 100
  533.       WRITE(IOCODE,21)I,J,K,VP(I,J,K),CT(I,J,K),BO(I,J,K),SO(I,J,K),
  534.      &BW(I,J,K),SW(I,J,K),BG(I,J,K),SG(I,J,K)
  535.   100 CONTINUE
  536.       WRITE(IOCODE,110)K,OOIP(K),OWIP(K),ODGIP(K),OFGIP(K)
  537.   110 FORMAT(/,T15,'LAYER',I3,' INITIAL FLUID VOLUMES:',
  538.      & /,T20,'OIL IN PLACE (MILLION STB)',T65,F10.4,
  539.      & /,T20,'WATER IN PLACE (MILLION STB)',T65,F10.4,
  540.      & /,T20,'SOLUTION GAS IN PLACE (BILLION SCF)',T65,F10.4,
  541.      & /,T20,'FREE GAS IN PLACE (BILLION SCF)',T65,F10.4/)
  542.   102 CONTINUE
  543.       TOOIP=0.
  544.       TOWIP=0.
  545.       TODGIP=0.
  546.       TOFGIP=0.
  547.       DO 103 K=1,KK
  548.       TOOIP=TOOIP+OOIP(K)
  549.       TOWIP=TOWIP+OWIP(K)
  550.       TODGIP=TODGIP+ODGIP(K)
  551.   103 TOFGIP=TOFGIP+OFGIP(K)
  552.       TOGIP = TODGIP + TOFGIP
  553.       WRITE(IOCODE,115) TOOIP,TOWIP,TODGIP,TOFGIP
  554.   115 FORMAT(/,T15,'TOTAL INITIAL FLUID VOLUMES IN RESERVOIR:',
  555.      & /,T20,'OIL IN PLACE (MILLION STB)',T65,F10.4,
  556.      & /,T20,'WATER IN PLACE (MILLION STB)',T65,F10.4,
  557.      & /,T20,'SOLUTION GAS IN PLACE (BILLION SCF)',T65,F10.4,
  558.      & /,T20,'FREE GAS IN PLACE (BILLION SCF)',T65,F10.4/)
  559.       STBO=SCFO*D5615
  560.       STBW=SCFW*D5615
  561.       MCFG=SCFG*0.001
  562.       MCFG1=SCFG1*0.001
  563.       STBOI=STBO
  564.       STBWI=STBW
  565.       MCFGT=MCFG+MCFG1
  566.       IF(MCFGI.LE.1.D-7.AND.MCFGT.LE.1.D-7)MBEG=0.0
  567.          NLOOP=N
  568.   105 IF(N.EQ.1.AND.ITFLAG.LE.0)
  569.      & CALL PRTPS(NLOOP,II,JJ,KK,PAVG0,PAVG,CMBEO,CMBEW,CMBEG,
  570.      & COP,CWP,CWI,CGP,CGI,MBEO,MBEW,MBEG,DELT0,
  571.      & OPR,WPR,GPR,WIR,GIR,ETI,
  572.      &CWOR,CGOR,WOR,GOR,IPMAP,ISOMAP,ISWMAP,ISGMAP,IPBMAP,IAQMAP)
  573.       IF(N.EQ.NMAX)GO TO 1006
  574. C      BEGIN OUTER ITERATION LOOP FOR STABILISED IMPES
  575.       ITNCNT=0
  576.  4000 ITNCNT=ITNCNT+1
  577. C**** ESTABLISH RATES & CALCULATE BHFP(IF PID IS NONZERO)
  578.       IF(NVQN.EQ.0)GO TO 1160
  579.       CALL QRATE(II,JJ,KK,NVQN,GORMAX,WORMAX,ETI)
  580. 1160   CONTINUE
  581. C**** CALCULATE SEVEN DIAGONAL MATRIX FOR FLOW EQ SOLUTION
  582.       IF(NUMDIS.EQ.0) CALL SOLONE(II,JJ,KK,DIV1,D288,
  583.      &KSM,KSM1,N,NN,KCOFF)
  584.       IF(NUMDIS.EQ.1) CALL SOLTWO(II,JJ,KK,DIV1,D288,
  585.      &KSM,KSM1,N,NN,KCOFF)
  586. C**** MODIFY MATRIX ELEMENTS FOR WELLS UNDER IMPLICIT CONTROL.
  587.       IF(NVQN.EQ.0) GO TO 1170
  588.       CALL PRATEI(NVQN)
  589.  1170 CONTINUE
  590. C    AQUIFER MODEL: P EQ COEF. MODIFICATIONS
  591.       IF(IAQOPT.NE.0) CALL AQIN(II,JJ,KK,DELT,ETI)
  592. C**** CALCULATE NEW PRESSURE DISTRIBUTION
  593. C      SKIP P CALC WHEN OUTER ITER EXCEEDS 1
  594.       IF(ITNCNT.GT.1) GO TO 1175
  595.       IF(KSOL.EQ.1) CALL GAUS1D(IOCODE,IFATAL,II,JJ,KK)
  596.       IF(IFATAL.GE.1) GO TO 1006
  597.       IF(KSOL.EQ.2) CALL LSORX(II,JJ,KK,OMEGA,TOL,TOL1,MITER,
  598.      & DELT,DELT0,KSN,NLOOP,NLITER)
  599.       IF(KSOL.EQ.3) CALL LSORY(II,JJ,KK,OMEGA,TOL,TOL1,MITER,
  600.      & DELT,DELT0,KSN,NLOOP,NLITER)
  601.       IF(KSOL.EQ.4) CALL LSORZ(II,JJ,KK,OMEGA,TOL,TOL1,MITER,
  602.      & DELT,DELT0,KSN,NLOOP,NLITER)
  603. C**** CALCULATE IMPLICIT RATES.
  604.  1175 CONTINUE
  605.       IF(NVQN.EQ.0) GO TO 1180
  606.       CALL PRATEO(NVQN)
  607.  1180 CONTINUE
  608. C** AQUIFER MODEL RATES
  609.       IF(IAQOPT.NE.0) CALL AQOUT(II,JJ,KK,DELT)
  610. C**** CALCULATE NEW FLUID SATURATIONS
  611.       SCFO=0.0
  612.       SCFW=0.0
  613.       SCFG=0.0
  614.       SCFG1=0.0
  615.       RESVOL=0.0
  616.       DO 400 K=1,KK
  617.       DO 400 J=1,JJ
  618.       DO 400 I=1,II
  619. C      CALC NEW RES PORE VOLUME
  620.       PPN=PN(I,J,K)
  621.       PP=P(I,J,K)
  622.       BPT=PBOT(I,J,K)
  623.       IPVTR=IPVT(I,J,K)
  624.       CALL INTPVT(IPVTR,BPT,RSLOPE(IPVTR),POT,RSOT,MPOT(IPVTR),PP,RSO)
  625.       CALL INTERP(IPVTR,PWT,RSWT,MPWT(IPVTR),PP,RSW)
  626.       CALL INTERP(IPVTR,PRT,CRT,MPGT(IPVTR),PPN,CR)
  627.       CALL INTPVT(IPVTR,BPT,BSLOPE(IPVTR),POT,BOT,MPOT(IPVTR),PP,BBO)
  628.       CALL INTERP(IPVTR,PWT,BWT,MPWT(IPVTR),PP,BBW)
  629.       CALL INTERP(IPVTR,PGT,BGT,MPGT(IPVTR),PP,BBG)
  630.       VPP=VP(I,J,K) * (1.0+CR*(P(I,J,K)-PPN))
  631.       RESVOL=RESVOL+VPP
  632.       DP1=0.0
  633.       DP2=0.0
  634.       DP3=0.0
  635.       DP4=0.0
  636.       DP5=0.0
  637.       DP6=0.0
  638.       IF((I-1).GT.0)  DP1=P(I-1,J,K)-PP
  639.       IF((I+1).LE.II) DP2=P(I+1,J,K)-PP
  640.       IF((J-1).GT.0)  DP3=P(I,J-1,K)-PP
  641.       IF((J+1).LE.JJ) DP4=P(I,J+1,K)-PP
  642.       IF((K-1).GT.0)  DP5=P(I,J,K-1)-PP
  643.       IF((K+1).LE.KK) DP6=P(I,J,K+1)-PP
  644.       DAODP=OW(I,J,K)*DP1+OE(I,J,K)*DP2+OS(I,J,K)*DP3
  645.      &     +ON(I,J,K)*DP4+OT(I,J,K)*DP5+OB(I,J,K)*DP6
  646.       DAWDP=WW(I,J,K)*DP1+WE(I,J,K)*DP2+WS(I,J,K)*DP3
  647.      &     +WN(I,J,K)*DP4+WT(I,J,K)*DP5+WB(I,J,K)*DP6
  648. C*** SKIP SAT. CALC. FOR 0. PV BLOCKS.  ASSUME NO SAT CHANGE.
  649.       IF(VPP.EQ.0.0) GO TO 405
  650.       SW(I,J,K)=((DAWDP+GWWT(I,J,K)-(QW(I,J,K)+EWAQ(I,J,K)))*DELT
  651.      & +VP(I,J,K)*SWN(I,J,K)/BW(I,J,K)) * BBW/VPP
  652.       IF(SW(I,J,K).GT.1.0) GO TO 401
  653.       SO(I,J,K)=((DAODP+GOWT(I,J,K)-QO(I,J,K))*DELT
  654.      & +VP(I,J,K)*SON(I,J,K)/BO(I,J,K)) * BBO/VPP
  655.       SG(I,J,K)=1.0-SO(I,J,K)-SW(I,J,K)
  656.       IF(SG(I,J,K).GT.0.0) GO TO 405
  657.   401 IF(SW(I,J,K).GT.1.0) SW(I,J,K)=1.0
  658.       SG(I,J,K)=0.0
  659.       SO(I,J,K)=1.0-SW(I,J,K)
  660.   405 CONTINUE
  661.       IF(KCOFF.EQ.0) GO TO 397
  662.       RHO1=VPP*SO(I,J,K)/BBO
  663.       RHO2=VP(I,J,K)*SON(I,J,K)/BO(I,J,K)
  664.       DIFFO=RHO1-RHO2
  665.       RHW1=VPP*SW(I,J,K)/BBW
  666.       RHW2=VP(I,J,K)*SWN(I,J,K)/BW(I,J,K)
  667.       DIFFW=RHW1-RHW2
  668.       RHG1=VPP*SG(I,J,K)/BBG
  669.       RHG2=VP(I,J,K)*SGN(I,J,K)/BG(I,J,K)
  670.       DIFFG=RHG1-RHG2
  671.       WRITE(IOCODE,33)
  672.       WRITE(IOCODE,21)I,J,K,P(I,J,K),SO(I,J,K),SON(I,J,K),SW(I,J,K),
  673.      &SWN(I,J,K),SG(I,J,K),SGN(I,J,K),VPP
  674.       WRITE(IOCODE,21)I,J,K,GOWT(I,J,K),QO(I,J,K),GWWT(I,J,K),
  675.      &QW(I,J,K)
  676.       WRITE(IOCODE,21)I,J,K,DAODP,DAWDP,DELT
  677.       WRITE(IOCODE,21)I,J,K,RHO1,RHO2,DIFFO
  678.       WRITE(IOCODE,21)I,J,K,RHW1,RHW2,DIFFW
  679.       WRITE(IOCODE,21)I,J,K,RHG1,RHG2,DIFFG,GGWT(I,J,K),QG(I,J,K)
  680. 21    FORMAT(1X,3I3,8E15.6)
  681.   397 CONTINUE
  682.   400 CONTINUE
  683. C      NEW SET OF INTERPOLATED VALUES FOR NEXT R-K OUTER ITER.
  684.       IF(IRK.EQ.0) GO TO 4300
  685.       DO 4200 I=1,II
  686.       DO 4200 J=1,JJ
  687.       DO 4200 K=1,KK
  688.       SO(I,J,K)=SON(I,J,K)+RKCOEF(ITNCNT,KRUNGE)*(SO(I,J,K)-SON(I,J,K))
  689.       SW(I,J,K)=SWN(I,J,K)+RKCOEF(ITNCNT,KRUNGE)*(SW(I,J,K)-SWN(I,J,K))
  690.       SG(I,J,K)=SGN(I,J,K)+RKCOEF(ITNCNT,KRUNGE)*(SG(I,J,K)-SGN(I,J,K))
  691. C      INTERPOLATE BUBBLE POINT PRESSURE FOR NEXT R-K ITERATION
  692.       IF(IREPRS.EQ.1) GO TO 4200
  693.       IBP=I
  694.       JBP=J
  695.       KBP=K
  696.       IF(IREPRS.EQ.0) CALL REPRS1(IBP,JBP,KBP)
  697.  4200 CONTINUE
  698.       IF(ITNCNT.EQ.KRUNGE) GO TO 4300
  699.       GO TO 4000
  700.  4300 CONTINUE
  701. C      UPDATE FVF AND COMPRESSIBILITY
  702.       DO 450 I=1,II
  703.       DO 450 J=1,JJ
  704.       DO 450 K=1,KK
  705.       PPN=PN(I,J,K)
  706.       PP=P(I,J,K)
  707.       BPT=PBOT(I,J,K)
  708.       IPVTR=IPVT(I,J,K)
  709.       CALL INTPVT(IPVTR,BPT,RSLOPE(IPVTR),POT,RSOT,MPOT(IPVTR),PP,RSO)
  710.       CALL INTERP(IPVTR,PWT,RSWT,MPWT(IPVTR),PP,RSW)
  711.       CALL INTERP(IPVTR,PRT,CRT,MPGT(IPVTR),PPN,CR)
  712.       CALL INTPVT(IPVTR,BPT,BSLOPE(IPVTR),POT,BOT,MPOT(IPVTR),PP,BBO)
  713.       CALL INTERP(IPVTR,PWT,BWT,MPWT(IPVTR),PP,BBW)
  714.       CALL INTERP(IPVTR,PGT,BGT,MPGT(IPVTR),PP,BBG)
  715.       VPP=VP(I,J,K) * (1.0+CR*(P(I,J,K)-PPN))
  716.       VP(I,J,K)=VPP
  717.       BO(I,J,K)=BBO
  718.       BW(I,J,K)=BBW
  719.       BG(I,J,K)=BBG
  720.       FF1=SO(I,J,K)/BO(I,J,K)
  721.       FF2=SW(I,J,K)/BW(I,J,K)
  722.       SCFO=SCFO+VP(I,J,K)*FF1
  723.       SCFW=SCFW+VP(I,J,K)*FF2
  724.       SCFG=SCFG+VP(I,J,K)*SG(I,J,K)/BG(I,J,K)
  725.       SCFG1=SCFG1+VP(I,J,K)*(RSO*FF1+RSW*FF2)
  726.       CALL INTCOM(IPVTR,POT,BOPT,MPOT(IPVTR),PP,BODER)
  727.       CALL INTCOM(IPVTR,POT,RSOPT,MPOT(IPVTR),PP,RSODER)
  728.       CALL INTCOM(IPVTR,PWT,BWPT,MPWT(IPVTR),PP,BWDER)
  729.       CALL INTCOM(IPVTR,PWT,RSWPT,MPWT(IPVTR),PP,RSWDER)
  730.       CALL INTCOM(IPVTR,PGT,BGPT,MPGT(IPVTR),PP,BGDER)
  731.       IF(PP.GT.PBOT(I,J,K))BODER=BSLOPE(IPVTR)
  732.       IF(PP.GT.PBOT(I,J,K))RSODER=RSLOPE(IPVTR)
  733.       CO=-(BODER-BG(I,J,K)*RSODER)/BO(I,J,K)
  734.       CW=-(BWDER-BG(I,J,K)*RSWDER)/BW(I,J,K)
  735.       CG=-BGDER/BG(I,J,K)
  736.       CALL INTERP(IPVTR,PRT,CRT,MPGT(IPVTR),PP,CR)
  737.       CT(I,J,K)=CR + CO*SO(I,J,K) +CW*SW(I,J,K) + CG*SG(I,J,K)
  738.       IF(N.EQ.KCO)WRITE(IOCODE,21)I,J,K,CR,CO,RSO,CW,RSW,CG
  739.   450 CONTINUE
  740. C      AUTO. TIME STEP CONTROL CALC. OF PRESSURE AND SAT. MAXIMA.
  741.       PPM=0.
  742.       SOM=0.
  743.       SWM=0.
  744.       SGM=0.
  745.       DO 240 K=1,KK
  746.       DO 240 J=1,JJ
  747.       DO 240 I=1,II
  748.       DPO=P(I,J,K)-PN(I,J,K)
  749.       DSO=SO(I,J,K)-SON(I,J,K)
  750.       DSW=SW(I,J,K)-SWN(I,J,K)
  751.       DSG=SG(I,J,K)-SGN(I,J,K)
  752.       IF(ABS(DPO).LE.ABS(PPM))GO TO 210
  753.       PPM=DPO
  754.       IPM=I
  755.       JPM=J
  756.       KPM=K
  757.   210 IF(ABS(DSO).LE.ABS(SOM))GO TO 220
  758.       SOM=DSO
  759.       IOM=I
  760.       JOM=J
  761.       KOM=K
  762.   220 IF(ABS(DSW).LE.ABS(SWM))GO TO 230
  763.       SWM=DSW
  764.       IWM=I
  765.       JWM=J
  766.       KWM=K
  767.   230 IF(ABS(DSG).LE.ABS(SGM))GO TO 240
  768.       SGM=DSG
  769.       IGM=I
  770.       JGM=J
  771.       KGM=K
  772.   240 CONTINUE
  773.       DPMC=ABS(PPM)
  774.       DSMC=ABS(SOM)
  775. C
  776.       ISM=IOM
  777.       JSM=JOM
  778.       KSM=KOM
  779.       IF(DSMC.GT.ABS(SGM)) GO TO 245
  780.       DSMC=ABS(SGM)
  781.       ISM=IGM
  782.       JSM=JGM
  783.       KSM=KGM
  784.   245 IF(DSMC.GT.ABS(SWM)) GO TO 248
  785.       DSMC=ABS(SWM)
  786.       ISM=IWM
  787.       JSM=JWM
  788.       KSM=KWM
  789.   248 CONTINUE
  790. C****REPEAT TIME STEP?
  791.       IF(DSMC.LT.DSMAX.AND.DPMC.LT.DPMAX)GO TO 402
  792.       IF(DELT.LE.DTMIN.OR.FACT2.GE.1.0)GO TO 402
  793.       ITFLAG=ITFLAG+1
  794.       DELT=DELT*FACT2
  795.       IF(DELT.LT.DTMIN) DELT=DTMIN
  796.       FT=ETI+DELT
  797.       IF(FT.GT.FTMAX)DELT=FTMAX-ETI
  798. C****RESET VARIABLES.
  799.       DO 260 I=1,II
  800.       DO 260 J=1,JJ
  801.       DO 260 K=1,KK
  802.       P(I,J,K)=PN(I,J,K)
  803.       SO(I,J,K)=SON(I,J,K)
  804.       SW(I,J,K)=SWN(I,J,K)
  805.       SG(I,J,K)=SGN(I,J,K)
  806.       PBOT(I,J,K)=PBOTN(I,J,K)
  807.   260 CONTINUE
  808.       GO TO 1060
  809.   402 CONTINUE
  810. C****UNDERSATURATED GRID BLOCK SATURATION CALCULATION.
  811. C      SKIP CALC IF NO OIL IN BLOCK
  812.       DO 410 I=1,II
  813.       DO 410 J=1,JJ
  814.       DO 410 K=1,KK
  815.       IF(SO(I,J,K).LE.0.0) GO TO 410
  816.       IF(P(I,J,K).GT.PN(I,J,K)) GO TO 410
  817.       IF(P(I,J,K).LT.PBOT(I,J,K)) GO TO 410
  818.       IP=I+1
  819.       IM=I-1
  820.       JP=J+1
  821.       JM=J-1
  822.       KP=K+1
  823.       KM=K-1
  824.       IF(IP.GT.II) GO TO 412
  825.       IF(SGN(IP,J,K).GT.0.0001) GO TO 410
  826.   412 IF(IM.LT.1) GO TO 414
  827.       IF(SGN(IM,J,K).GT.0.0001) GO TO 410
  828.   414 IF(JP.GT.JJ) GO TO 416
  829.       IF(SGN(I,JP,K).GT.0.0001) GO TO 410
  830.   416 IF(JM.LT.1) GO TO 418
  831.       IF(SGN(I,JM,K).GT.0.0001) GO TO 410
  832.   418 IF(KP.GT.KK) GO TO 420
  833.       IF(SGN(I,J,KP).GT.0.0001) GO TO 410
  834.   420 IF(KM.LT.1) GO TO 422
  835.       IF(SGN(I,J,KM).GT.0.0001) GO TO 410
  836.   422 SG(I,J,K)=0.0
  837.       SO(I,J,K)=1.0-SW(I,J,K)
  838.   410 CONTINUE
  839. C**** REPRESSURIZATION ALGORITHM.
  840.       IF(IREPRS.EQ.1) GO TO 5011
  841.       DO 50 I=1,II
  842.       DO 50 J=1,JJ
  843.       DO 50 K=1,KK
  844.       IBP=I
  845.       JBP=J
  846.       KBP=K
  847.       IF(IREPRS.EQ.0) CALL REPRS1(IBP,JBP,KBP)
  848.    50 CONTINUE
  849. 5011  CONTINUE
  850. C**** UPDATE OLD FLUID VOLUMES FOR MATERIAL BALANCE.
  851.       STBOI=STBO
  852.       STBWI=STBW
  853.       MCFGI=MCFGT
  854. C**** UPDATE  NEW FLUID VOLUMES.
  855.       STBO=SCFO*D5615
  856.       STBW=SCFW*D5615
  857.       MCFG=SCFG*0.001
  858.       MCFG1=SCFG1*0.001
  859.       MCFGT=MCFG+MCFG1
  860. C***** DEBUG PRINT OF PRESENT AND FUTURE P,SO,SW,SG VALUES.
  861.       IF(KCOFF.EQ.0)GO TO 2901
  862.       DO 290 K=1,KK
  863.       DO 290 J=1,JJ
  864.       DO 290 I=1,II
  865.       WRITE(IOCODE,21)I,J,K,PN(I,J,K),SON(I,J,K),SWN(I,J,K),SGN(I,J,K)
  866.       WRITE(IOCODE,21)I,J,K,P(I,J,K),SO(I,J,K),SW(I,J,K),SG(I,J,K)
  867. 290   CONTINUE
  868. 2901  CONTINUE
  869. C     AQUIFER MODEL: CUMULATIVES
  870.       IF(IAQOPT.NE.0) CALL AQCUM(NLOOP)
  871. C**** WELL REPORT (ALL RATE & PRESSURE DATA APPLICABLE THIS STEP)
  872.       IJ=0
  873.       TOR=0.
  874.       TGR=0.
  875.       TWR=0.
  876.       TOC=0.
  877.       TGC=0.
  878.       TWC=0.
  879.       IFREQW=IFREQW+1
  880.       DO 2040 J=1,NVQN
  881.       IQ1=IQN1(J)
  882.       IQ2=IQN2(J)
  883.       IQ3=IQN3(J)
  884.       IJLOC=IDWELL(J)
  885.       IF(IJLOC.EQ.0) GO TO 2040
  886.       IJ=IJ+1
  887.       LAY=IQ3+(LAYER(J)-1)
  888.       DO 2040 K=IQ3,LAY
  889.       GOR=0.
  890.       WOR=0.
  891.       QOO=QOC(IJLOC,K)*D5615
  892.       QWW=QWC(IJLOC,K)*D5615
  893.       QGG=QGC(IJLOC,K)*0.001
  894.       CUMO(J,K)=CUMO(J,K)+QOO*DELT*0.001
  895.       CUMW(J,K)=CUMW(J,K)+QWW*DELT*0.001
  896.       CUMG(J,K)=CUMG(J,K)+QGG*DELT*0.001
  897.       IF(IOMETH.GT.0) GO TO 992
  898.       IF(IWLREP.GT.IFREQW.AND.FT.LT.FTMAX) GO TO 2040
  899.       GO TO 994
  900.   992 IF(FT.LT.FTMAX) GO TO 2040
  901.   994 CONTINUE
  902.       IF(IJ.NE.1.OR.K.NE.IQ3) GO TO 996
  903.       WRITE(IOCODE,70) (IHEDIN(IH),IH=1,40)
  904.       WRITE(IOCODE,591) N,FT
  905.       WRITE(IOCODE,592)
  906.   591 FORMAT(/,T45,41('*'),/,T45,'*',T85,'*',/,T45,
  907.      & '*  WELL REPORT: BOAST II (RELEASE 1.1)  *',/,
  908.      & T45,'*',T85,'*',/,T45,41('*'),//,
  909.      & T5,10('*'),' WELL REPORT FOR TIME STEP NUMBER ',I3,
  910.      &1X,'  ELAPSED TIME =',F9.3,' DAYS FROM BEGINNING OF SIMULATION ',
  911.      & 10('*'),//)
  912.   592 FORMAT(/,T56,'------ RATE ------',22X,'--- CUMULATIVE ---',
  913.      & /,6X,'   WELL     LOCATION',4X,'CALC    SPEC    SPEC',4X,
  914.      & 'OIL     GAS    WATER    GOR     WOR',5X,
  915.      & 'OIL     GAS    WATER',/,
  916.      & 8X,'#     ID',3X,'I  J  K    BHFP    BHFP     PI',
  917.      & 3X,'STB/D    MCF/D   STB/D',20X,'MSTB    MMCF    MSTB',/)
  918.   996 CONTINUE
  919.       IF(QOO.EQ.0.)GO TO 998
  920.       GOR=QGG*1000./QOO
  921.       WOR=QWW/QOO
  922.   998 CONTINUE
  923.       FLDCUM=ABS(CUMO(J,K))+ABS(CUMG(J,K))+ABS(CUMW(J,K))
  924. C      IF FLDCUM LE 0, SKIP OUTPUT
  925.       IF(FLDCUM.LE.0.) GO TO 999
  926.       WRITE(IOCODE,593) IDWELL(J),WELNAM(J),IQN1(J),IQN2(J),K,PWFC(J,K),
  927.      & PWF(J,K),PID(J,K),QOO,QGG,QWW,GOR,WOR,CUMO(J,K),CUMG(J,K),
  928.      & CUMW(J,K)
  929.   593 FORMAT(4X,I5,2X,A5,1X,3I3,2F8.0,F7.3,F9.1,F9.0,F9.1,F7.1,F7.3,F8.1,
  930.      &F8.0,F8.1)
  931.   999 CONTINUE
  932. C      TIME STEP SUMMARY HEADING ACTIVATED
  933.       ITSS=0
  934.       TOR=TOR+QOO
  935.       TGR=TGR+QGG
  936.       TWR=TWR+QWW
  937.       TOC=TOC+CUMO(J,K)
  938.       TGC=TGC+CUMG(J,K)
  939.       TWC=TWC+CUMW(J,K)
  940.  2040 CONTINUE
  941.       IF(IOMETH.GT.0) GO TO 2542
  942.       IF(IWLREP.GT.IFREQW.AND.FT.LT.FTMAX) GO TO 2552
  943.       IFREQW=0
  944.       GO TO 2544
  945.  2542 IF(FT.LT.FTMAX) GO TO 2552
  946.  2544 CONTINUE
  947.       WRITE(IOCODE,594)TOR,TGR,TWR,TOC,TGC,TWC
  948.   594 FORMAT(6X,108('-'),/,
  949.      &12X,'TOTALS',31X,3F9.1,14X,3F8.1,/)
  950.       WRITE(IOCODE,2551)
  951.  2551 FORMAT(///T5,49('*'),'  END OF WELL REPORT  ',49('*'),6(/))
  952.  2552 CONTINUE
  953. C*****CALCULATE MATERIAL BALANCE ERRORS & AVERAGE RESERVOIR PRESSURE
  954.       DELT0=DELT
  955.       ETI=ETI+DELT
  956.       CALL MATBAL(II,JJ,KK,STBO,STBOI,STBW,STBWI,TOWIP,TOOIP,TOGIP,
  957.      &MCFGI,MBEO,MBEW,MBEG,CMBEO,CMBEW,CMBEG,DELT0,RESVOL,
  958.      &OP,WP,GP,WI,GI,PAVG0,PAVG,OPR,WPR,GPR,WIR,GIR,D5615,
  959.      &COP,CWP,CGP,CWI,CGI,MCFGT,CWOR,WOR,CGOR,GOR)
  960.       IF(WOR.GT.WORMAX)GO TO 1002
  961.       IF(1000.0*GOR.GT.GORMAX)GO TO 1003
  962.       IF(PAVG.LT.PAMIN)GO TO 1004
  963.       IF(PAVG.GT.PAMAX)GO TO 1005
  964. C***** SUMMARY REPORT.
  965.       IFREQS=IFREQS+1
  966.       IF(IOMETH.GT.0) GO TO 2554
  967.       IF(ISUMRY.GT.IFREQS.AND.FT.LT.FTMAX) GO TO 2557
  968.       IFREQS=0
  969.       GO TO 2556
  970.  2554 IF(FT.LT.FTMAX) GO TO 2557
  971.  2556 NLP=N+1
  972. C      TIME STEP SUMMARY HEADING ACTIVATED
  973.       ITSS=0
  974.       WRITE(IOCODE,70) (IHEDIN(IH),IH=1,40)
  975.       CALL PRTPS(NLP,II,JJ,KK,PAVG0,PAVG,CMBEO,CMBEW,CMBEG,
  976.      &COP,CWP,CWI,CGP,CGI,MBEO,MBEW,MBEG,DELT0,
  977.      &OPR,WPR,GPR,WIR,GIR,ETI,
  978.      &CWOR,CGOR,WOR,GOR,IPMAP,ISOMAP,ISWMAP,ISGMAP,IPBMAP,IAQMAP)
  979.  2557 IF(N.NE.KCO.OR.KCO1.EQ.0)GO TO 500
  980.       DO 300 K=1,KK
  981.       DO 300 J=1,JJ
  982.       DO 300 I=1,II
  983.       WRITE(IOCODE,21)I,J,K,VP(I,J,K),CT(I,J,K),BO(I,J,K),SO(I,J,K),
  984.      &BW(I,J,K),SW(I,J,K),BG(I,J,K),SG(I,J,K)
  985. 300   CONTINUE
  986.   500 CONTINUE
  987.       IF(N.EQ.KSN)KSN=KSN+KSN1
  988.       IF(N.EQ.KSM)KSM=KSM+KSM1
  989.       IF(N.EQ.KCO)KCO=KCO+KCO1
  990. C**** UPDATE ARRAYS.
  991.       DO 1150 K=1,KK
  992.       DO 1150 J=1,JJ
  993.       DO 1150 I=1,II
  994.       QO(I,J,K)=0.0
  995.       QW(I,J,K)=0.0
  996.       QG(I,J,K)=0.0
  997.        PN(I,J,K)=P(I,J,K)
  998.        SON(I,J,K)=SO(I,J,K)
  999.        SWN(I,J,K)=SW(I,J,K)
  1000.        SGN(I,J,K)=SG(I,J,K)
  1001.        PBOTN(I,J,K)=PBOT(I,J,K)
  1002. 1150  CONTINUE
  1003. C**  REQUIRED ONE LINE TIME STEP SUMMARY
  1004.       IF(ITSS.GT.0) GO TO 1275
  1005.       ITSS=1
  1006.       WRITE(IOCODE,1200)
  1007.  1200 FORMAT(1H1/T54,29('*'),/,
  1008.      & T54,'*   TIME   STEP   SUMMARY   *',/,
  1009.      & T54,29('*'),//,
  1010.      & T2,'TIME STEP',15X,'PRODUCTION',21X,'INJECTION',
  1011.      & 5X,'PV WT   MATERIAL  BALANCES',2X,
  1012.      & 'MAX SATN CHANGE MAX PRES CHANGE  ITN ',/,
  1013.      & T2,9('-'),1X,40('-'),1X,16('-'),2X,' AVG',3X,
  1014.      & 20('-'),1X,15('-'),1X,15('-'),1X,5('-'))
  1015.       WRITE(IOCODE,1250)
  1016.  1250 FORMAT(T3,13X,'OIL     GAS     WATER    GOR   WATER',
  1017.      & 3X,' GAS     WATER   RES     OIL    GAS',
  1018.      & 3X,'WATER  I  J  K  DSMAX  I  J  K  DPMAX O   I',/,
  1019.      & T3,38X,'SCF/   /OIL',20X,'PRES',5X,'%',2(6X,'%'),/,
  1020.      & T1,2X,'NO. DAYS    STB/D  MSCF/D    STB/D',
  1021.      & 3X,' STB   RATIO  MSCF/D    STB/D   PSIA',/,
  1022.      & T2,4('-'),1X,4('-'),1X,8('-'),2(2X,7('-')),
  1023.      & 1X,6('-'),2X,5('-'),1X,8('-'),1X,7('-'),1X,6('-'),2X,3(6('-'),1X),
  1024.      & 3(2('-'),1X),6('-'),1X,3(2('-'),1X),6('-'),
  1025.      & 1X,5('-'),/)
  1026.  1275 GORS=0.
  1027.       WORS=0.
  1028.       IF(OPR.GT.0.) GORS=GPR*1000./OPR
  1029.       IF(OPR.GT.0.) WORS=WPR/OPR
  1030. C      TOTAL RUN SUMMARY SAVED INFORMATION
  1031.       ITSMAX=N
  1032.       ITSNO(N)=N
  1033.       STIME(N)=ETI
  1034.       SOPROD(N)=OPR
  1035.       SGPROD(N)=GPR
  1036.       SWPROD(N)=WPR
  1037.       SGOR(N)=GORS
  1038.       SWOR(N)=WORS
  1039.       SGINJ(N)=GIR
  1040.       SWINJ(N)=WIR
  1041.       SPVWTP(N)=PAVG
  1042.       NM=N-1
  1043.       IF(NM.EQ.0) NM=1
  1044.       SOCUMP(N)=SOPROD(N)*DELT/1000. + SOCUMP(NM)
  1045.       SWCUMP(N)=SWPROD(N)*DELT/1000. + SWCUMP(NM)
  1046.       SGCUMP(N)=SGPROD(N)*DELT/1000. + SGCUMP(NM)
  1047.       SWCUMI(N)=SWINJ(N)*DELT/1000. + SWCUMI(NM)
  1048.       SGCUMI(N)=SGINJ(N)*DELT/1000. + SGCUMI(NM)
  1049.       WRITE(IOCODE,1280) N,ETI,OPR,GPR,WPR,GORS,WORS,GIR,WIR,PAVG,
  1050.      & MBEO,MBEG,MBEW,ISM,JSM,KSM,DSMC,IPM,JPM,KPM,DPMC,
  1051.      & ITNCNT,NLITER
  1052.  1280 FORMAT(I4,F6.0,1X,F8.1,1X,F8.1,1X,F8.1,F7.1,           
  1053.      & 1X,F5.1,1X,F9.1,1X,F7.0,1X,F6.0,3(1X,F6.3),3(1X,I2),1X,F6.3,1X,
  1054.      & 3(1X,I2),1X,F6.2,1X,I1,1X,I3)
  1055. C      END OF TIME?
  1056.       IF(FT.GE.TMAX) GO TO 1006
  1057. 1000  CONTINUE
  1058. 1002  WRITE(IOCODE,2902)
  1059.       GO TO 1006
  1060. 1003  WRITE(IOCODE,2903)
  1061.       GO TO 1006
  1062. 1004  WRITE(IOCODE,2904)
  1063.       GO TO 1006
  1064. 1005  WRITE(IOCODE,2905)
  1065. C******  EXIT POINT
  1066. 1006  CONTINUE
  1067.       IF(IFATAL.LT.1) GO TO 2020
  1068.       WRITE(IOCODE,1010) IFATAL
  1069.  1010 FORMAT(//5X,5('-'),'FATAL INPUT ERRORS HAVE BEEN',
  1070.      & ' DETECTED (# =',I9,'). RUN TERMINATED.',5('-'),/)
  1071.       GO TO 2250
  1072.  2020 CONTINUE
  1073. C      POST-PLOT PACKAGE CALL
  1074.       IF(IPLOTP.GE.0) CALL POSTP(NPLINE,IOCODE)
  1075.  2250 CONTINUE
  1076.       WRITE(IOCODE,2300)
  1077.  2300 FORMAT(///,T4,40('*'),
  1078.      & ' BOAST II (RELEASE 1.1) SIMULATION TERMINATED ',
  1079.      & 40('*'),//)
  1080. 2902  FORMAT(/T15,'MAXIMUM WOR HAS BEEN EXCEEDED --- SIMULATION',
  1081.      &' IS BEING TERMINATED'//)
  1082. 2903  FORMAT(/T15,'MAXIMUM GOR HAS BEEN EXCEEDED --- SIMULATION',
  1083.      &' IS BEING TERMINATED'//)
  1084. 2904  FORMAT(/T15,'MINIMUM AVERAGE RESERVOIR PRESSURE WAS NOT',
  1085.      &' ACHIEVED --- SIMULATION IS BEING TERMINATED'//)
  1086.  2905 FORMAT(/T15,'MAXIMUM AVERAGE RESERVOIR PRESSURE
  1087.      1 HAS BEEN EXCEEDED --- SIMULATION IS BEING TERMINATED'//)
  1088. C      CLOSE I/O FILES FOR DEC 20/60.
  1089. C      MACHINE DEPENDENT CLOSE STATEMENTS
  1090.       CLOSE(UNIT=20)
  1091.       CLOSE(UNIT=21)
  1092.       CLOSE(UNIT=24)
  1093.       CLOSE(UNIT=25)
  1094.       STOP
  1095.       END
  1096. C..................................................BLOCK DATA
  1097.       BLOCK DATA 
  1098.       PARAMETER (LP1=13,LP2=10,LP3=4)
  1099.       PARAMETER (LP7=3,LP8=3,LP9=25,LP10=3,LP11=50,LP12=1000)
  1100.       PARAMETER (LP14=1,LP15=13,LP17=5)
  1101.       PARAMETER (LP4=LP1+1,LP5=LP2+1,LP6=LP3+1,LP13=LP1+LP2+LP3)
  1102.       PARAMETER (LP19=LP4*LP2*LP3,LP20=LP1*LP5*LP3)
  1103.       PARAMETER (LP21=LP1*LP2*LP6,LP22=LP1*LP2*LP3,LP23=LP11*LP3)
  1104.       REAL KROT,KROGT,KRWT,KRGT,MUWT,MUOT,MUGT,KX,KY,KZ
  1105.      & ,MUO,MUW,MUG,KRO,KRW,MBEWI,MBEGI,MCFGI
  1106.      & ,MBEO,MBEW,MBEG,MCFG,MBEOI,MIN,MCFG1,MCFGT
  1107.       COMMON /COEF/ AW(LP1,LP2,LP3),AE(LP1,LP2,LP3),AN(LP1,LP2,LP3),
  1108.      & AS(LP1,LP2,LP3),AB(LP1,LP2,LP3),AT(LP1,LP2,LP3),E(LP1,LP2,LP3),
  1109.      & B(LP1,LP2,LP3)
  1110.       COMMON /SAQUI/ IAQOPT,CPIAQ1(LP1,LP2,LP3),CPIAQ2(LP1,LP2,LP3),
  1111.      & CPI1(LP1,LP2,LP3),CPI2(LP1,LP2,LP3),EWAQ(LP1,LP2,LP3),
  1112.      & CUMAQW(LP1,LP2,LP3),
  1113.      & QWAQ(LP1,LP2,LP3),CUMEW(LP1,LP2,LP3),QWAQR(LP7),CUMAQR(LP7)
  1114.      & ,IAQREG(LP1,LP2,LP3),PAQ(LP1,LP2,LP3),PIAQ(LP1,LP2,LP3)   
  1115.         COMMON /RUNSUM/ ITSNO(LP12),STIME(LP12),SOPROD(LP12),
  1116.      & SGPROD(LP12),
  1117.      & SWPROD(LP12),SGOR(LP12),SWOR(LP12),SGINJ(LP12),SWINJ(LP12)
  1118.       COMMON /RUN2/SPVWTP(LP12),SOCUMP(LP12),SWCUMP(LP12),SGCUMP(LP12),
  1119.      & SGCUMI(LP12),SWCUMI(LP12),SAQUIR(LP12),SAQUIC(LP12)
  1120.       COMMON /SPARM/ KX(LP1,LP2,LP3),KY(LP1,LP2,LP3),KZ(LP1,LP2,LP3),
  1121.      & EL(LP1,LP2,LP3),TX(LP4,LP2,LP3),TY(LP1,LP5,LP3),TZ(LP1,LP2,LP6),
  1122.      & PDAT(LP1,LP2,LP3),PDATUM,GRAD
  1123.       COMMON /SPVT/ SAT(LP7,LP9),KROT(LP7,LP9),KRWT(LP7,LP9),
  1124.      & BGT(LP7,LP9),
  1125.      & KRGT(LP7,LP9),ITHREE(LP7),RSOT(LP7,LP9),BWPT(LP7,LP9),
  1126.      & PCOWT(LP7,LP9),PCGOT(LP7,LP9),KROGT(LP7,LP9),SWR(LP7),
  1127.      & POT(LP7,LP9),MUOT(LP7,LP9),BOT(LP7,LP9),BOPT(LP7,LP9),
  1128.      & RSOPT(LP7,LP9),PWT(LP7,LP9),MUWT(LP7,LP9),BWT(LP7,LP9),
  1129.      & RSWT(LP7,LP9),RSWPT(LP7,LP9),PGT(LP7,LP9),MUGT(LP7,LP9),
  1130.      & BGPT(LP7,LP9),CRT(LP7,LP9),IPVT(LP1,LP2,LP3),IROCK(LP1,LP2,LP3),
  1131.      & NROCK,NPVT,PSIT(LP7,LP9),PRT(LP7,LP9),WOROCK(LP7),GOROCK(LP7)
  1132.       COMMON /SRATE/ PID(LP11,LP3),PWF(LP11,LP3),PWFC(LP11,LP3),
  1133.      & KIP(LP11),LAYER(LP11),QVO(LP11),CUMG(LP11,LP3),
  1134.      & GMO(LP11,LP3),GMW(LP11,LP3),GMG(LP11,LP3),
  1135.      & QVW(LP11),QVG(LP11),QVT(LP11),CUMO(LP11,LP3),CUMW(LP11,LP3),
  1136.      & IDWELL(LP11),ALIT(LP11),BLIT(LP11)
  1137.       COMMON /SSOLN/ BO(LP1,LP2,LP3),BW(LP1,LP2,LP3),BG(LP1,LP2,LP3),
  1138.      & QO(LP1,LP2,LP3),QW(LP1,LP2,LP3),QG(LP1,LP2,LP3),
  1139.      & GOWT(LP1,LP2,LP3),GWWT(LP1,LP2,LP3),GGWT(LP1,LP2,LP3),
  1140.      & OW(LP4,LP2,LP3),OE(LP4,LP2,LP3),WW(LP4,LP2,LP3),WE(LP4,LP2,LP3),
  1141.      & OS(LP1,LP5,LP3),ON(LP1,LP5,LP3),WS(LP1,LP5,LP3),WN(LP1,LP5,LP3),
  1142.      & OT(LP1,LP2,LP6),OB(LP1,LP2,LP6),WT(LP1,LP2,LP6),WB(LP1,LP2,LP6),
  1143.      & QOWG(LP1,LP2,LP3),VP(LP1,LP2,LP3),CT(LP1,LP2,LP3)
  1144.       COMMON /TSTDAT/ IFATAL,IWARN
  1145.       DATA IFATAL,IWARN/0,0/
  1146.        DATA NROCK,NPVT /LP7,LP8/
  1147.       DATA TX,TY,TZ/LP19*0.0,LP20*0.0,LP21*0.0/
  1148.       DATA    CUMO,CUMW,CUMG/LP23*0.0,LP23*0.0,LP23*0.0/
  1149.       DATA            AW,AS,AT/LP22*0.0,LP22*0.0,LP22*0.0/
  1150.       DATA            AE,AN,AB/LP22*0.0,LP22*0.0,LP22*0.0/
  1151.       DATA  OW,OE,WW,WE/LP19*0.0,LP19*0.0,LP19*0.0,LP19*0.0/
  1152.       DATA  OS,ON,WS,WN/LP20*0.0,LP20*0.0,LP20*0.0,LP20*0.0/
  1153.       DATA  OT,OB,WT,WB/LP21*0.0,LP21*0.0,LP21*0.0,LP21*0.0/
  1154.       DATA     IROCK,IPVT/LP22*1,LP22*1/
  1155.       DATA     CPIAQ1,CPIAQ2/LP22*0.,LP22*0./
  1156.       END
  1157.