home *** CD-ROM | disk | FTP | other *** search
/ Hall of Fame / HallofFameCDROM.cdr / oilfield / spe-46-2.lzh / BLOCK6.FOR < prev    next >
Text File  |  1988-07-14  |  40KB  |  1,090 lines
  1. $DO66
  2. C.................................................................TABLE
  3.       SUBROUTINE TABLE(IOCODE,II,JJ,KK)
  4. C      MACHINE DEPENDENT INCLUDE STATEMENT
  5. $INCLUDE:'PARAMS.FOR'
  6. C      READ ROCK AND PVT DATA
  7.       REAL KROT,KROGT,KRWT,KRGT,MUOT,MUWT,MUGT,KX,KY,KZ
  8.       COMMON /BUBBLE/ PBO,VSLOPE(LP8),BSLOPE(LP8),RSLOPE(LP8),PMAXT,
  9.      & IREPRS,MPGT(LP8),
  10.      & RHOSCO(LP8),RHOSCW(LP8),RHOSCG(LP8),MSAT(LP7),MPOT(LP8),
  11.      & MPWT(LP8),PBOT(LP1,LP2,LP3),PBOTN(LP1,LP2,LP3)
  12.       COMMON /SPARM/ KX(LP1,LP2,LP3),KY(LP1,LP2,LP3),KZ(LP1,LP2,LP3),
  13.      & EL(LP1,LP2,LP3),TX(LP4,LP2,LP3),TY(LP1,LP5,LP3),TZ(LP1,LP2,LP6),
  14.      & PDAT(LP1,LP2,LP3),PDATUM,GRAD
  15.       COMMON /SPVT/ SAT(LP7,LP9),KROT(LP7,LP9),KRWT(LP7,LP9),
  16.      & BGT(LP7,LP9),
  17.      & KRGT(LP7,LP9),ITHREE(LP7),RSOT(LP7,LP9),BWPT(LP7,LP9),
  18.      & PCOWT(LP7,LP9),PCGOT(LP7,LP9),KROGT(LP7,LP9),SWR(LP7),
  19.      & POT(LP7,LP9),MUOT(LP7,LP9),BOT(LP7,LP9),BOPT(LP7,LP9),
  20.      & RSOPT(LP7,LP9),PWT(LP7,LP9),MUWT(LP7,LP9),BWT(LP7,LP9),
  21.      & RSWT(LP7,LP9),RSWPT(LP7,LP9),PGT(LP7,LP9),MUGT(LP7,LP9),
  22.      & BGPT(LP7,LP9),CRT(LP7,LP9),IPVT(LP1,LP2,LP3),IROCK(LP1,LP2,LP3),
  23.      & NROCK,NPVT,PSIT(LP7,LP9),PRT(LP7,LP9),WOROCK(LP7),GOROCK(LP7)
  24.       COMMON /TSTDAT/ IFATAL,IWARN
  25.       WRITE(IOCODE,111)
  26.       READ(20,1)
  27.       READ(20,*) NROCK,NPVT
  28.       WRITE(IOCODE,1050) NROCK,NPVT
  29.  1050 FORMAT(//,T15,'*** NUMBER OF REGIONS ***',/,
  30.      & T20,'ROCK',T30,I5,/,
  31.      & T20,'PVT',T30,I5,//)
  32.       IF(NROCK.LE.1) GO TO 2000
  33.       READ(20,1)
  34.       READ(20,*) NUMROK
  35.       IF(NUMROK.EQ.0) GO TO 1500
  36.       WRITE(IOCODE,1100)
  37.       DO 1200 L=1,NUMROK
  38.       READ(20,*) I1,I2,J1,J2,K1,K2,IVAL
  39.       WRITE(IOCODE,1150) I1,I2,J1,J2,K1,K2,IVAL
  40.       DO 1200 K=K1,K2
  41.       DO 1200 J=J1,J2
  42.       DO 1200 I=I1,I2
  43.       IROCK(I,J,K)=IVAL
  44.  1200 CONTINUE
  45.  1500 WRITE(IOCODE,1600)
  46.       DO 1700 K=1,KK
  47.       WRITE(IOCODE,1650) K
  48.       DO 1700 J=1,JJ
  49.       IF(NUMROK.EQ.0) READ(20,*) (IROCK(I,J,K),I=1,II)
  50.       WRITE(IOCODE,1675) (IROCK(I,J,K),I=1,II)
  51.  1700 CONTINUE
  52.  2000 CONTINUE
  53.       IF(NPVT.LE.1) GO TO 3000
  54.       READ(20,1)
  55.       READ(20,*) NUMPVT
  56.       IF(NUMPVT.EQ.0) GO TO 2500
  57.       WRITE(IOCODE,2100)
  58.       DO 2200 L=1,NUMPVT
  59.       READ(20,*) I1,I2,J1,J2,K1,K2,IVAL
  60.       WRITE(IOCODE,1150) I1,I2,J1,J2,K1,K2,IVAL
  61.       DO 2200 K=K1,K2
  62.       DO 2200 J=J1,J2
  63.       DO 2200 I=I1,I2
  64.       IPVT(I,J,K)=IVAL
  65.  2200 CONTINUE
  66.  2500 WRITE(IOCODE,2600)
  67.       DO 2700 K=1,KK
  68.       WRITE(IOCODE,1650) K
  69.       DO 2700 J=1,JJ
  70.       IF(NUMPVT.EQ.0) READ(20,*) (IPVT(I,J,K),I=1,II)
  71.       WRITE(IOCODE,1675) (IPVT(I,J,K),I=1,II)
  72.  2700 CONTINUE
  73.  3000 CONTINUE
  74.  1100 FORMAT(//T15,'***** ROCK REGION MODIFICATIONS *****',
  75.      & //T15,'   I1  I2  J1  J2  K1  K2  IROCK')
  76.  1150 FORMAT(15X,6I4,3X,I3)
  77.  1600 FORMAT(//T15,'***** ROCK REGION DISTRIBUTION *****',/)
  78.  1650 FORMAT(/1X,'K =',I4/)
  79.  1675 FORMAT(1X,20I6)
  80.  2100 FORMAT(//T15,'***** PVT REGION MODIFICATIONS *****',
  81.      & //T15,'   I1  I2  J1  J2  K1  K2  IPVT')
  82.  2600 FORMAT(//T15,'***** PVT REGION DISTRIBUTION *****',/)
  83. C****  RELATIVE PERMEABILITY & CAPILLARY PRESSURE TABLE
  84.       DO 10 NR=1,NROCK
  85.       READ(20,1)
  86.       WRITE(IOCODE,7) NR
  87.       DO 5 I=1, 25
  88.       READ(20,*) SAT(NR,I),KROT(NR,I),KRWT(NR,I),KRGT(NR,I),
  89.      & KROGT(NR,I),PCOWT(NR,I),PCGOT(NR,I)
  90.       WRITE(IOCODE,21)SAT(NR,I),KROT(NR,I),KRWT(NR,I),KRGT(NR,I),
  91.      & KROGT(NR,I),PCOWT(NR,I),PCGOT(NR,I)
  92.       IF(SAT(NR,I).GE.1.0001)GO TO 9
  93.     5 CONTINUE
  94.     9 MSAT(NR)=I
  95.       READ(20,1)
  96.       READ(20,*) ITHREE(NR),SWR(NR)
  97.       IF(ITHREE(NR).EQ.0) WRITE(IOCODE,92) SWR(NR)
  98.       IF(ITHREE(NR).GT.0) WRITE(IOCODE,94) SWR(NR)
  99.    92 FORMAT(//,10X,'THE THREE-PHASE REL. PERM. CALC. IS NOT USED.',
  100.      & /,10X,'IRREDUCIBLE WATER SATURATION IS ',F6.4)
  101.    94 FORMAT(//,10X,'THE THREE-PHASE REL. PERM. CALC. IS USED.',
  102.      & /,10X,'IRREDUCIBLE WATER SATURATION IS ',F6.4)
  103.    10 CONTINUE
  104. C  SAT DEPENDENT DATA CHECK
  105.       DO 4200 NR=1,NROCK
  106.       DO 4200 I=1,MSAT(NR)
  107.       IF(SAT(NR,I).GE.-0.1.AND.SAT(NR,I).LE.1.1) GO TO 4020
  108.       IFATAL=IFATAL+1
  109.       WRITE(IOCODE,4010) NR,I
  110.  4010 FORMAT(/5X,5('-'),'SAT ERROR FOR ROCK REGION',I5,
  111.      & ' AND ENTRY # ',I5)
  112.  4020 IF(KROT(NR,I).GE.0.0.AND.KROT(NR,I).LE.1.0) GO TO 4040
  113.       IFATAL=IFATAL+1
  114.       WRITE(IOCODE,4030) NR,I
  115.  4030 FORMAT(/5X,5('-'),'KROW ERROR FOR ROCK REGION',I5,
  116.      & ' AND ENTRY # ',I5)
  117.  4040 IF(KRWT(NR,I).GE.0.0.AND.KRWT(NR,I).LE.1.0) GO TO 4060
  118.       IFATAL=IFATAL+1
  119.       WRITE(IOCODE,4050) NR,I
  120.  4050 FORMAT(/5X,5('-'),'KRW ERROR FOR ROCK REGION',I5,
  121.      & ' AND ENTRY # ',I5)
  122.  4060 IF(KRGT(NR,I).GE.0.0.AND.KRGT(NR,I).LE.1.0) GO TO 4080
  123.       IFATAL=IFATAL+1
  124.       WRITE(IOCODE,4070) NR,I
  125.  4070 FORMAT(/5X,5('-'),'KRG ERROR FOR ROCK REGION',I5,
  126.      & ' AND ENTRY # ',I5)
  127.  4080 IF(KROGT(NR,I).GE.0.0.AND.KROGT(NR,I).LE.1.0) GO TO 4100
  128.       IFATAL=IFATAL+1
  129.       WRITE(IOCODE,4090) NR,I
  130.  4090 FORMAT(/5X,5('-'),'KROG ERROR FOR ROCK REGION',I5,
  131.      & ' AND ENTRY # ',I5)
  132.  4100 IF(PCOWT(NR,I).GE.0.0) GO TO 4120
  133.       IFATAL=IFATAL+1
  134.       WRITE(IOCODE,4110) NR,I
  135.  4110 FORMAT(/5X,5('-'),'PCOW ERROR FOR ROCK REGION',I5,
  136.      & ' AND ENTRY # ',I5)
  137.  4120 IF(PCGOT(NR,I).GE.0.0) GO TO 4200
  138.       IFATAL=IFATAL+1
  139.       WRITE(IOCODE,4130) NR,I
  140.  4130 FORMAT(/5X,5('-'),'PCGO ERROR FOR ROCK REGION',I5,
  141.      & ' AND ENTRY # ',I5)
  142.  4200 CONTINUE
  143.       DO 4300 NR=1,NROCK
  144.       IF(SAT(NR,1).LT.0.0) GO TO 4220
  145.       IFATAL=IFATAL+1
  146.       WRITE(IOCODE,4210) NR
  147.  4210 FORMAT(/5X,5('-'),'FIRST SAT ENTRY FOR ROCK REGION',
  148.      & I5,' SHOULD EQUAL -0.10.',/10X,
  149.      & 'ADJUST REL PERM AND CAP PRESS CURVES ALSO.',/)
  150.  4220 IF(SAT(NR,MSAT(NR)).GE.1.05) GO TO 4240
  151.       IFATAL=IFATAL+1
  152.       WRITE(IOCODE,4230) NR
  153.  4230 FORMAT(/5X,5('-'),'LAST SAT ENTRY FOR ROCK REGION',
  154.      & I5,' SHOULD EQUAL 1.10.',/10X,
  155.      & 'ADJUST REL PERM AND CAP PRESS CURVES ALSO.',/)
  156. C      CHECK 3-PHASE SWR DATA
  157.  4240 IF(SWR(NR).GE.0.0.AND.SWR(NR).LE.1.0) GO TO 4300
  158.       IFATAL=IFATAL+1
  159.       WRITE(IOCODE,4250) NR
  160.  4250 FORMAT(/5X,5('-'),'3-PHASE SWR ERROR FOR ROCK REGION',I5)
  161.  4300 CONTINUE
  162. C****  BUBBLE POINT & MAXIMUM PRESSURES
  163.       DO 400 NP=1,NPVT
  164. C      INITIALIZE BUBBLE POINT PRESSURE ARRAY.
  165.       READ(20,1)
  166.       READ(20,*) PBO,PBODAT,PBGRAD
  167.       WRITE(IOCODE,3100) NP,PBO,PBODAT,PBGRAD
  168.  3100 FORMAT(1H1/T15,'PVT REGION # ',I5,//,
  169.      & //T10,'BUBBLE POINT PARAMETERS:',
  170.      & /,T5,'BUBBLE POINT PRESSURE (PSI)',T35,F10.3,
  171.      & /,T5,'BUBBLE POINT DEPTH (FT)',T35,F10.3,
  172.      & /,T5,'BUBBLE POINT GRADIENT (PSI/FT)',T35,F10.3)
  173.       DO 3500 K=1,KK
  174.       DO 3500 J=1,JJ
  175.       DO 3500 I=1,II
  176.       IF(IPVT(I,J,K).NE.NP) GO TO 3500
  177.       PBOT(I,J,K)=PBO+(PBODAT-EL(I,J,K))*PBGRAD
  178.       PBOTN(I,J,K)=PBOT(I,J,K)
  179. C      CHECK BUBBLE POINT PRESSURE DATA
  180.       IF(PBOT(I,J,K).GE.0.0) GO TO 3500
  181.       IFATAL=IFATAL+1
  182.       WRITE(IOCODE,3300) I,J,K
  183.  3300 FORMAT(/5X,5('-'),'BUBBLE POINT PRESSURE ERROR IN',
  184.      & ' GRID BLOCK IJK = ',3I5)
  185.  3500 CONTINUE
  186. C****  OIL DATA TABLE
  187.       READ(20,1)
  188.       WRITE(IOCODE,8)
  189.       READ(20,*) VSLOPE(NP),BSLOPE(NP),RSLOPE(NP),PMAXT,IREPRS
  190.       WRITE(IOCODE,31)VSLOPE(NP),BSLOPE(NP),RSLOPE(NP),PMAXT,IREPRS
  191. C      CHECK UNDERSATURATED SLOPES
  192.       IF(VSLOPE(NP).GE.0.0) GO TO 3520
  193.       IWARN=IWARN+1
  194.       WRITE(IOCODE,3510) NP
  195.  3510 FORMAT(/5X,5('-'),'VSLOPE FOR PVT REGION',I5,
  196.      & ' IS NEGATIVE')
  197.  3520 IF(BSLOPE(NP).LE.0.0) GO TO 3540
  198.       IWARN=IWARN+1
  199.       WRITE(IOCODE,3530) NP
  200.  3530 FORMAT(/5X,5('-'),'BSLOPE FOR PVT REGION',I5,
  201.      & ' IS POSITIVE')
  202.  3540 IF(RSLOPE(NP).EQ.0.0) GO TO 3560
  203.       IWARN=IWARN+1
  204.       WRITE(IOCODE,3550) NP
  205.  3550 FORMAT(/5X,5('-'),'RSLOPE FOR PVT REGION',I5,
  206.      & 'IS NOT ZERO')
  207.  3560 CONTINUE
  208.       READ(20,1)
  209.       WRITE(IOCODE,12)
  210.       DO 15 I=1, 25
  211.       READ(20,*) POT(NP,I),MUOT(NP,I),BOT(NP,I),RSOT(NP,I)
  212.       WRITE(IOCODE,23)POT(NP,I),MUOT(NP,I),BOT(NP,I),RSOT(NP,I)
  213.       RSOT(NP,I)=0.17809*RSOT(NP,I)
  214.       IF(POT(NP,I).GE.PMAXT)GO TO 20
  215. 15    CONTINUE
  216. 20    MPOT(NP)=I
  217. C      POT  DEPENDENT DATA CHECK
  218.       DO 5200 I=1,MPOT(NP)
  219.       IF(I.GT.1) GO TO 5105
  220.       IF(POT(NP,1).GE.0.0.AND.POT(NP,1).LE.PMAXT) GO TO 5120
  221.       IFATAL=IFATAL+1
  222.       WRITE(IOCODE,5110) NP,I
  223.       GO TO 5120
  224.  5105 IF(POT(NP,I).GE.POT(NP,I-1).AND.POT(NP,I).LE.PMAXT) GO TO 5120
  225.       IFATAL=IFATAL+1
  226.       WRITE(IOCODE,5110) NP,I
  227.  5110 FORMAT(/5X,5('-'),'POT ERROR FOR PVT REGION',I5,
  228.      & ' AND ENTRY # ',I5)
  229.  5120 IF(MUOT(NP,I).GE.0.0) GO TO 5151
  230.       IFATAL=IFATAL+1
  231.       WRITE(IOCODE,5130) NP,I
  232.  5130 FORMAT(/5X,5('-'),'MUO ERROR FOR PVT REGION',I5,
  233.      & ' AND ENTRY # ',I5)
  234.  5151 IF(BOT(NP,I).GE.0.0) GO TO 5160
  235.       IFATAL=IFATAL+1
  236.       WRITE(IOCODE,5150) NP,I
  237.  5150 FORMAT(/5X,5('-'),'BO ERROR FOR PVT REGION',I5,
  238.      & ' AND ENTRY # ',I5)
  239.  5160 IF(RSOT(NP,I).GE.0.0) GO TO 5200
  240.       IFATAL=IFATAL+1
  241.       WRITE(IOCODE,5170) NP,I
  242.  5170 FORMAT(/5X,5('-'),'RSO ERROR FOR PVT REGION',I5,
  243.      & ' AND ENTRY # ',I5)
  244.  5200 CONTINUE
  245. C****  WATER DATA TABLE
  246.       READ(20,1)
  247.       WRITE(IOCODE,13)
  248.       DO 25 I=1, 25
  249.       READ(20,*) PWT(NP,I),MUWT(NP,I),BWT(NP,I),RSWT(NP,I)
  250.       WRITE(IOCODE,23)PWT(NP,I),MUWT(NP,I),BWT(NP,I),RSWT(NP,I)
  251.       RSWT(NP,I)=0.17809*RSWT(NP,I)
  252.       IF(PWT(NP,I).GE.PMAXT)GO TO 30
  253. 25    CONTINUE
  254. 30    MPWT(NP)=I
  255. C      PWT DEPENDENT DATA CHECK
  256.       DO 5400 I=1,MPWT(NP)
  257.       IF(I.GT.1) GO TO 5305
  258.       IF(PWT(NP,1).GE.0.0.AND.PWT(NP,1).LE.PMAXT) GO TO 5320
  259.       IFATAL=IFATAL+1
  260.       WRITE(IOCODE,5310) NP,I
  261.       GO TO 5320
  262.  5305 IF(PWT(NP,I).GE.PWT(NP,I-1).AND.PWT(NP,I).LE.PMAXT) GO TO 5320
  263.       IFATAL=IFATAL+1
  264.       WRITE(IOCODE,5310) NP,I
  265.  5310 FORMAT(/5X,5('-'),'PWT ERROR FOR PVT REGION',I5,
  266.      & ' AND ENTRY # ',I5)
  267.  5320 IF(MUWT(NP,I).GE.0.0) GO TO 5353
  268.       IFATAL=IFATAL+1
  269.       WRITE(IOCODE,5330) NP,I
  270.  5330 FORMAT(/5X,5('-'),'MUW ERROR FOR PVT REGION',I5,
  271.      & ' AND ENTRY # ',I5)
  272.  5353 IF(BWT(NP,I).GE.0.0) GO TO 5360
  273.       IFATAL=IFATAL+1
  274.       WRITE(IOCODE,5350) NP,I
  275.  5350 FORMAT(/5X,5('-'),'BW ERROR FOR PVT REGION',I5,
  276.      & ' AND ENTRY # ',I5)
  277.  5360 IF(RSWT(NP,I).GE.0.0) GO TO 5400
  278.       IFATAL=IFATAL+1
  279.       WRITE(IOCODE,5370) NP,I
  280.  5370 FORMAT(/5X,5('-'),'RSW ERROR FOR PVT REGION',I5,
  281.      & ' AND ENTRY # ',I5)
  282.  5400 CONTINUE
  283. C****  GAS & ROCK DATA TABLE
  284.       READ(20,1)
  285.       READ(20,*) KGCOR
  286.       READ(20,1)
  287.       IF(KGCOR.EQ.1) GO TO 50
  288.       WRITE(IOCODE,14)
  289.       DO 35 I=1,999
  290.       READ(20,*) PGT(NP,I),MUGT(NP,I),BGT(NP,I),PSIT(NP,I),CRT(NP,I)
  291.       WRITE(IOCODE,24)PGT(NP,I),MUGT(NP,I),BGT(NP,I),PSIT(NP,I),
  292.      & CRT(NP,I)
  293.       PRT(NP,I)=PGT(NP,I)
  294.       IF(PGT(NP,I).GE.PMAXT)GO TO 40
  295. 35    CONTINUE
  296. 40    MPGT(NP)=I
  297.       GO TO 99
  298.    50 CONTINUE
  299.       NPHASE=NP
  300.       CALL PSEUDO(NPHASE,IFATAL,IOCODE)
  301.       WRITE(IOCODE,51)
  302.    51 FORMAT(//,T6,2X,'**** ROCK COMPRESSIBILITY ****',
  303.      & /T5,1X,'PRESSURE',4X,'ROCK COMP.',
  304.      & /T5,2X,'(PSI)',6X,'(1/PSI)',/)
  305.       READ(20,1)
  306.       DO 55 IG=1,MPGT(NP)
  307.       IF(IG.GT.1.AND.PRT(NP,1).EQ.PMAXT) GO TO 53
  308.       READ(20,*) PRT(NP,IG),CRT(NP,IG)
  309.       WRITE(IOCODE,52) PRT(NP,IG),CRT(NP,IG)
  310.    52 FORMAT(3X,F10.1,2X,E10.3)
  311.       GO TO 55
  312.    53 CONTINUE
  313.       PRT(NP,IG)=PGT(NP,IG)
  314.       CRT(NP,IG)=CRT(NP,1)
  315.       IF(IG.EQ.MPGT(NP)) PRT(NP,1)=PGT(NP,1)
  316.    55 CONTINUE
  317.    99 CONTINUE
  318. C      PGT DEPENDENT DATA CHECK
  319.       DO 5600 I=1,MPGT(NP)
  320.       IF(I.GT.1) GO TO 5505
  321.       IF(PGT(NP,1).GE.0.0.AND.PGT(NP,1).LE.PMAXT) GO TO 5520
  322.       IFATAL=IFATAL+1
  323.       WRITE(IOCODE,5510) NP,I
  324.       GO TO 5520
  325.  5505 IF(PGT(NP,I).GE.PGT(NP,I-1).AND.
  326.      & (PGT(NP,I)-PMAXT).LE.0.1) GO TO 5520
  327.       IFATAL=IFATAL+1
  328.       WRITE(IOCODE,5510) NP,I
  329.  5510 FORMAT(/5X,5('-'),'PGT ERROR FOR PVT REGION',I5,
  330.      & ' AND ENTRY # ',I5)
  331.  5520 IF(MUGT(NP,I).GE.0.0) GO TO 5555
  332.       IFATAL=IFATAL+1
  333.       WRITE(IOCODE,5530) NP,I
  334.  5530 FORMAT(/5X,5('-'),'MUG ERROR FOR PVT REGION',I5,
  335.      & ' AND ENTRY # ',I5)
  336.  5555 IF(BGT(NP,I).GE.0.0) GO TO 5560
  337.       IFATAL=IFATAL+1
  338.       WRITE(IOCODE,5550) NP,I
  339.  5550 FORMAT(/5X,5('-'),'BG ERROR FOR PVT REGION',I5,
  340.      & ' AND ENTRY # ',I5)
  341.  5560 IF(CRT(NP,I).GE.0.0) GO TO 5580
  342.       IFATAL=IFATAL+1
  343.       WRITE(IOCODE,5570) NP,I
  344.  5570 FORMAT(/5X,5('-'),'CR ERROR FOR PVT REGION',I5,
  345.      & ' AND ENTRY # ',I5)
  346. 5580  IF(PRT(NP,I).GE.0.0.AND.(PRT(NP,I)-PMAXT).LE.0.1) GO TO 5590
  347.        IFATAL=IFATAL+1
  348.        WRITE(IOCODE,5585) NP,I
  349.  5585  FORMAT(/5X,5('-'),'ROCK PRES ERROR FOR PVT REGION',I5,
  350.      & ' AND ENTRY # ',I5)
  351.  5590 IF(PSIT(NP,I).GE.0.0) GO TO 5596
  352.       IFATAL=IFATAL+1
  353.       WRITE(IOCODE,5595) NP,I
  354.  5595 FORMAT(/5X,5('-'),'PSEUDO-PRES ERROR FOR PVT REGION',I5,
  355.      & ' AND ENTRY # ',I5)
  356.  5596 IF(I.EQ.1) GO TO 5600
  357.       IF(PRT(NP,I).GE.PRT(NP,I-1).AND.(PRT(NP,I)-PMAXT).LE.0.1)
  358.      &  GO TO 5600
  359.       IFATAL=IFATAL+1
  360.       WRITE(IOCODE,5585) NP,I
  361.  5600 CONTINUE
  362. C****  OIL, WATER, & GAS DENSITIES AT STOCK TANK CONDITIONS
  363.       READ(20,1)
  364.       WRITE(IOCODE,6)
  365.       READ(20,*)RHOSCO(NP),RHOSCW(NP),RHOSCG(NP)
  366.       WRITE(IOCODE,4)RHOSCO(NP),RHOSCW(NP),RHOSCG(NP)
  367. C      CHECK DENSITY DATA
  368.       IF(RHOSCO(NP).GT.0.0) GO TO 5720
  369.       IFATAL=IFATAL+1
  370.       WRITE(IOCODE,5710) NP
  371.  5710 FORMAT(/5X,5('-'),'OIL DENSITY FOR PVT REGION',I5,
  372.      & ' IN ERROR')
  373.  5720 IF(RHOSCW(NP).GT.0.0) GO TO 5740
  374.       IFATAL=IFATAL+1
  375.       WRITE(IOCODE,5730) NP
  376.  5730 FORMAT(/5X,5('-'),'WATER DENSITY FOR PVT REGION',I5,
  377.      & ' IN ERROR')
  378.  5740 IF(RHOSCG(NP).GT.0.0) GO TO 5760
  379.       IFATAL=IFATAL+1
  380.       WRITE(IOCODE,5750) NP
  381.  5750 FORMAT(/5X,5('-'),'GAS DENSITY  FOR PVT REGION',I5,
  382.      & ' IN ERROR')
  383.  5760 CONTINUE
  384. C**** CALCULATE SLOPES DBO/DP, DRSO/DP, DBW/DP, DRSW/DP, DBG/DP
  385. C**** NOTE: SLOPE AT POINT "I" IS BASED ON POINTS "I" AND "I-1"
  386.       WRITE(IOCODE,222)
  387.       WRITE(IOCODE,223)
  388.       DO 100 I=2,MPOT(NP)
  389.       DIV=1.0/(POT(NP,I)-POT(NP,I-1))
  390.       BOPT(NP,I)=(BOT(NP,I)-BOT(NP,I-1))*DIV
  391.       RSOPT(NP,I)=(RSOT(NP,I)-RSOT(NP,I-1))*DIV
  392. 100   WRITE(IOCODE,224)POT(NP,I),BOT(NP,I),BOPT(NP,I),RSOT(NP,I)
  393.      & ,RSOPT(NP,I)
  394.       WRITE(IOCODE,333)
  395.       DO 200 I=2,MPWT(NP)
  396.       DIV=1.0/(PWT(NP,I)-PWT(NP,I-1))
  397.       BWPT(NP,I)=(BWT(NP,I)-BWT(NP,I-1))*DIV
  398.       RSWPT(NP,I)=(RSWT(NP,I)-RSWT(NP,I-1))*DIV
  399. 200   WRITE(IOCODE,224)PWT(NP,I),BWT(NP,I),BWPT(NP,I),RSWT(NP,I)
  400.      & ,RSWPT(NP,I)
  401.       WRITE(IOCODE,444)
  402.       DO 300 I=2,MPGT(NP)
  403.       BGPT(NP,I)=(BGT(NP,I)-BGT(NP,I-1))/(PGT(NP,I)-PGT(NP,I-1))
  404. 300   WRITE(IOCODE,445)PGT(NP,I),BGT(NP,I),BGPT(NP,I)
  405. C      CHECK FLUID COMPRESSIBILITIES
  406.       DO 5820 I=2,MPOT(NP)
  407.       PP=POT(NP,I)
  408.       IPVTR=NP
  409.       CALL INTERP(IPVTR,PGT,BGT,MPGT(IPVTR),PP,BBG)
  410.       CALL INTCOM(IPVTR,POT,BOPT,MPOT(IPVTR),PP,BODER)
  411.       CALL INTCOM(IPVTR,POT,RSOPT,MPOT(IPVTR),PP,RSODER)
  412.       COTST=-(BODER-BBG*RSODER)
  413.       IF(COTST.GE.0.0) GO TO 5820
  414.       IFATAL=IFATAL+1
  415.       WRITE(IOCODE,5810) NP,I
  416.  5810 FORMAT(/5X,5('-'),'OIL COMP ERROR FOR PVT REGION',I5,
  417.      & ' AND ENTRY # ',I5)
  418.  5820 CONTINUE
  419.       DO 5840 I=2,MPWT(NP)
  420.       PP=PWT(NP,I)
  421.       IPVTR=NP
  422.       CALL INTERP(IPVTR,PGT,BGT,MPGT(IPVTR),PP,BBG)
  423.       CALL INTCOM(IPVTR,PWT,BWPT,MPWT(IPVTR),PP,BWDER)
  424.       CALL INTCOM(IPVTR,PWT,RSWPT,MPWT(IPVTR),PP,RSWDER)
  425.       CWTST=-(BWDER-BBG*RSWDER)
  426.       IF(CWTST.GE.0.0) GO TO 5840
  427.       IFATAL=IFATAL+1
  428.       WRITE(IOCODE,5830) NP,I
  429.  5830 FORMAT(/5X,5('-'),'WATER COMP ERROR FOR PVT REGION',I5,
  430.      & ' AND ENTRY # ',I5)
  431.  5840 CONTINUE
  432.       DO 5900 I=2,MPGT(NP)
  433.       PP=PGT(NP,I)
  434.       IPVTR=NP
  435.       CALL INTCOM(IPVTR,PGT,BGPT,MPGT(IPVTR),PP,BGDER)
  436.       CGTST=-BGDER
  437.       IF(CGTST.GE.0.0) GO TO 5900
  438.       IFATAL=IFATAL+1
  439.       WRITE(IOCODE,5850) NP,I
  440.  5850 FORMAT(/5X,5('-'),'GAS COMP ERROR FOR PVT REGION',I5,
  441.      & ' AND ENTRY # ',I5)
  442.  5900 CONTINUE
  443.   400 CONTINUE
  444. 1     FORMAT(40A2)
  445.     4 FORMAT(9X,8F10.4)
  446.     6 FORMAT(//,T8,'*** DENSITIES AT STD. CONDITIONS ***',/,
  447.      & 14X,'OIL',6X,'WATER',6X,'GAS',/,
  448.      & 10X,3(1X,'(LBM/SCF)'),/)
  449.     7 FORMAT(1H1/T15,'ROCK REGION # ',I5,//,
  450.      & 2X,'SATURATION',4X,'KROT',6X,'KRW',7X,'KRG',
  451.      & 8X,'KROG',5X,'PCOW',6X,'PCGO',/,
  452.      & T57,'(PSI)',T67,'(PSI)',/)
  453.     8 FORMAT(//,3X,'VSLOPE',5X,'BSLOPE',7X,'RSLOPE',
  454.      & 5X,'PMAX',2X,'REPRS',/,
  455.      & 2X,'(CP/PSI)',1X,'(RB/STB/PSI)',4X,
  456.      & '(1/PSI)',3X,'(PSI)',/)
  457.    21 FORMAT(1X,7F10.4)
  458. 23    FORMAT(3X,F10.1,1X,F8.4,2X,F8.4,2X,F8.2)
  459. 24    FORMAT(3X,F10.1,1X,F8.4,2X,E10.4,2X,E10.4,2X,E10.3)
  460.    12 FORMAT(//,T5,1X,'**** SATURATED OIL PVT PROPERTIES ****',
  461.      & /,T5,1X,'PRESSURE',2X,'VISCOSITY',4X,'FVF',4X,'SOLN. GAS',
  462.      & /,T5,2X,'(PSI)',7X,'(CP)',4X,'(RB/STB)',2X,'(SCF/STB)',/)
  463.    13 FORMAT(//,T5,5X,'**** WATER PVT PROPERTIES ***',
  464.      & /,T5,1X,'PRESSURE',2X,'VISCOSITY',4X,'FVF',4X,'SOLN. GAS',
  465.      & /,T5,2X,'(PSI)',7X,'(CP)',4X,'(RB/STB)',2X,'(SCF/STB)',/)
  466.    14 FORMAT(//,T6,2X,'****  GAS  PVT AND ROCK COMP. ****',
  467.      & /,T5,1X,'PRESSURE',2X,'VISCOSITY',4X,'FVF',5X,'PSEUDO-PRS',
  468.      & 2X,'ROCK COMP.',
  469.      & /,T5,2X,'(PSI)',7X,'(CP)',3X,'(RCF/SCF)',' (PSIA**2/CP)',
  470.      & 2X,'(1/PSI)',/)
  471. 111   FORMAT(1H1/T15,'***** EMPIRICAL DATA TABLES *****')
  472.    31 FORMAT(1X,E10.3,1X,E10.3,1X,2F10.2,1X,I5)
  473. 222   FORMAT(//T15,'***** SLOPES FOR COMPRESSIBILITY CALCULATIONS',
  474.      &' ****')
  475. 223   FORMAT(//T12,'PRESSURE',T25,'BO',T34,'DBO/DP',T47,'RSO',
  476.      & T55,'DRSO/DP'/,
  477.      & T13,'(PSI)',T22,'(RB/STB)',T31,'(RB/STB/PSI)',T45,
  478.      & '(CF/CF)',T55,'(1/PSI)',/)
  479. 445   FORMAT(10X,F9.1,1X,2E15.4)
  480. 333   FORMAT(//T12,'PRESSURE',T25,'BW',T34,'DBW/DP',T47,'RSW',
  481.      & T55,'DRSW/DP'/,
  482.      & T13,'(PSI)',T22,'(RB/STB)',T31,'(RB/STB/PSI)',T45,
  483.      & '(CF/CF)',T55,'(1/PSI)',/)
  484. 224   FORMAT( T12,F7.1,F10.4,2X,E11.4,F9.1,2X,E11.4)
  485. 444   FORMAT(//12X,'PRESSURE',9X,'BG',12X,'DBG/DP',/,
  486.      & T14,'(PSI)',T27,'(RCF/SCF)',T40,'(RCF/SCF/PSI)',/)
  487.       RETURN
  488.       END
  489. C........................................................................TRANS
  490.       SUBROUTINE TRANS(II,JJ,KK)
  491. C      MACHINE DEPENDENT INCLUDE STATEMENT
  492. $INCLUDE:'PARAMS.FOR'
  493. C      TRANSMISSIBILITY CALC
  494.       REAL KX,KY,KZ
  495.       COMMON /TSTDAT/ IFATAL,IWARN
  496.       COMMON /SARRAY/ PN(LP1,LP2,LP3),IOCODE,IDMAX,
  497.      & SON(LP1,LP2,LP3),SWN(LP1,LP2,LP3),SGN(LP1,LP2,LP3),
  498.      & A1(LP1,LP2,LP3),A2(LP1,LP2,LP3),A3(LP1,LP2,LP3),
  499.      & SUM(LP1,LP2,LP3),GAM(LP1,LP2,LP3),QS(LP1,LP2,LP3)
  500.       COMMON /SPARM/ KX(LP1,LP2,LP3),KY(LP1,LP2,LP3),KZ(LP1,LP2,LP3),
  501.      & EL(LP1,LP2,LP3),TX(LP4,LP2,LP3),TY(LP1,LP5,LP3),TZ(LP1,LP2,LP6),
  502.      & PDAT(LP1,LP2,LP3),PDATUM,GRAD
  503.       COMMON /SSOLN/ BO(LP1,LP2,LP3),BW(LP1,LP2,LP3),BG(LP1,LP2,LP3),
  504.      & QO(LP1,LP2,LP3),QW(LP1,LP2,LP3),QG(LP1,LP2,LP3),
  505.      & GOWT(LP1,LP2,LP3),GWWT(LP1,LP2,LP3),GGWT(LP1,LP2,LP3),
  506.      & OW(LP4,LP2,LP3),OE(LP4,LP2,LP3),WW(LP4,LP2,LP3),WE(LP4,LP2,LP3),
  507.      & OS(LP1,LP5,LP3),ON(LP1,LP5,LP3),WS(LP1,LP5,LP3),WN(LP1,LP5,LP3),
  508.      & OT(LP1,LP2,LP6),OB(LP1,LP2,LP6),WT(LP1,LP2,LP6),WB(LP1,LP2,LP6),
  509.      & QOWG(LP1,LP2,LP3),VP(LP1,LP2,LP3),CT(LP1,LP2,LP3)
  510.       COMMON /VECTOR/ DX(LP1,LP2,LP3),DY(LP1,LP2,LP3),DZ(LP1,LP2,LP3),
  511.      & DZNET(LP1,LP2,LP3),IQN1(LP11),IQN2(LP11),IQN3(LP11),IHEDIN(80)
  512.       DO 30 K=1,KK
  513.       DO 30 J=1,JJ
  514.       DO 30 I=1,II
  515.       FACX=1.0
  516.       FACY=1.0
  517.       FACZ=1.0
  518.       IF(I.EQ.1.OR.I.EQ.II) GO TO 20
  519.       XDENOM=2.*DX(I,J,K)+DX(I+1,J,K)+DX(I-1,J,K)
  520.       IF(XDENOM.GT.0.0) FACX=4.*DX(I,J,K)/XDENOM
  521.    20 IF(J.EQ.1.OR.J.EQ.JJ) GO TO 24
  522.       YDENOM=2.*DY(I,J,K)+DY(I,J+1,K)+DY(I,J-1,K)
  523.       IF(YDENOM.GT.0.0) FACY=4.*DY(I,J,K)/YDENOM
  524.    24 IF(K.EQ.1.OR.K.EQ.KK) GO TO 28
  525.       ZDENOM=2.*DZNET(I,J,K)+DZNET(I,J,K+1)+DZNET(I,J,K-1)
  526.       IF(ZDENOM.GT.0.0) FACZ=4.*DZNET(I,J,K)/ZDENOM
  527.    28 CONTINUE
  528.       A1(I,J,K)=FACX*KX(I,J,K)*DY(I,J,K)*DZNET(I,J,K)
  529.       A2(I,J,K)=FACY*KY(I,J,K)*DX(I,J,K)*DZNET(I,J,K)
  530.       A3(I,J,K)=FACZ*KZ(I,J,K)*DX(I,J,K)*DY(I,J,K)
  531.    30 CONTINUE
  532.       IF(II.EQ.1)GO TO 501
  533.       DO 50 K=1,KK
  534.       DO 50 J=1,JJ
  535.       DO 50 I=2,II
  536.       XDENOM=DX(I-1,J,K)*A1(I,J,K)+DX(I,J,K)*A1(I-1,J,K)
  537.       IF(XDENOM.EQ.0.0) GO TO 50
  538.       TX(I,J,K)=.012656*A1(I-1,J,K)*A1(I,J,K)/XDENOM
  539. 50    CONTINUE
  540. 501   CONTINUE
  541.       IF(JJ.EQ.1)GO TO 551
  542.       DO 55 K=1,KK
  543.       DO 55 J=2,JJ
  544.       DO 55 I=1,II
  545.       YDENOM=DY(I,J-1,K)*A2(I,J,K)+DY(I,J,K)*A2(I,J-1,K)
  546.       IF(YDENOM.EQ.0.0) GO TO 55
  547.       TY(I,J,K)=.012656*A2(I,J-1,K)*A2(I,J,K)/YDENOM
  548. 55    CONTINUE
  549. 551   CONTINUE
  550.       IF(KK.EQ.1)GO TO 601
  551.       DO 60 K=2,KK
  552.       DO 60 J=1,JJ
  553.       DO 60 I=1,II
  554.       ZDENOM=DZNET(I,J,K-1)*A3(I,J,K)+DZNET(I,J,K)*A3(I,J,K-1)
  555.       IF(ZDENOM.EQ.0.0) GO TO 60
  556.       TZ(I,J,K)=.012656*A3(I,J,K-1)*A3(I,J,K)/ZDENOM
  557. 60    CONTINUE
  558. 601   CONTINUE
  559. C**  COMPUTE BLOCK PV
  560.       DO 62 K=1,KK
  561.       DO 62 J=1,JJ
  562.       DO 62 I=1,II
  563.       VP(I,J,K)=VP(I,J,K)*DX(I,J,K)*DY(I,J,K)*DZNET(I,J,K)
  564.    62 CONTINUE
  565. C**  DEFINE 0 TRANS. FOR ALL BOUNDARIES OF 0 PV BLOCK.
  566.       DO 68 K=1,KK
  567.       DO 68 J=1,JJ
  568.       DO 68 I=1,II
  569.       IF(VP(I,J,K).GT.0.) GO TO 68
  570.       TX(I,J,K)=0.
  571.       TX(I+1,J,K)=0.
  572.       TY(I,J,K)=0.
  573.       TY(I,J+1,K)=0.
  574.       TZ(I,J,K)=0.
  575.       TZ(I,J,K+1)=0.
  576.    68 CONTINUE
  577. C****** TRANSMISSIBILITY MODIFICATIONS
  578.       READ(20,69)
  579.    69 FORMAT(40A2)
  580.       READ(20,*) NUMTX,NUMTY,NUMTZ,ITCODE
  581.       IF(NUMTX.EQ.0) GO TO 100
  582.       WRITE(IOCODE,31)
  583.    31 FORMAT(//T15,'**********TRANSMISSIBILITY (TX) NODE MODIFICATIONS',
  584.      & '**********',//T15,
  585.      & '   I1  I2  J1  J2  K1  K2  NEW TX VALUE')
  586.       DO 275 L=1,NUMTX
  587.       READ(20,*) I1,I2,J1,J2,K1,K2,REGVAL
  588.       WRITE(IOCODE,32) I1,I2,J1,J2,K1,K2,REGVAL
  589.       DO 275 K=K1,K2
  590.       DO 275 J=J1,J2
  591.       DO 275 I=I1,I2
  592.       IF(VP(I,J,K).GT.0.) TX(I,J,K)=REGVAL
  593.   275 CONTINUE
  594.    32 FORMAT(15X,6I4,2X,E10.4)
  595.   100 IF(ITCODE.NE.1)GO TO 8531
  596.       WRITE(IOCODE,44)
  597.    44 FORMAT(//T15,'**********TRANSMISSIBILITY (TX) DISTRIBUTION'
  598.      & ,'**********'/)
  599.       DO 853 K=1,KK
  600.       WRITE(IOCODE,38)K
  601.    38 FORMAT(/1X,'K =',I2/)
  602.       DO 854 J=1,JJ
  603.   854 WRITE(IOCODE,72)(TX(I,J,K),I=1,II)
  604.    72 FORMAT(1X,15F8.1)
  605. 853   CONTINUE
  606. 8531  CONTINUE
  607.       IF(NUMTY.EQ.0) GO TO 110
  608.       WRITE(IOCODE,34)
  609.    34 FORMAT(//T15,'**********TRANSMISSIBILITY (TY) NODE MODIFICATIONS',
  610.      & '**********',//T15,
  611.      & '   I1  I2  J1  J2  K1  K2  NEW TY VALUE')
  612.       DO 276 L=1,NUMTY
  613.       READ(20,*) I1,I2,J1,J2,K1,K2,REGVAL
  614.       WRITE(IOCODE,32) I1,I2,J1,J2,K1,K2,REGVAL
  615.       DO 276 K=K1,K2
  616.       DO 276 J=J1,J2
  617.       DO 276 I=I1,I2
  618.       IF(VP(I,J,K).GT.0.) TY(I,J,K)=REGVAL
  619.   276 CONTINUE
  620.   110 IF(ITCODE.NE.1) GO TO 8551
  621.       WRITE(IOCODE,47)
  622.    47 FORMAT(//T15,'**********TRANSMISSIBILITY (TY) DISTRIBUTION'
  623.      & ,'**********'/)
  624.       DO 855 K=1,KK
  625.       WRITE(IOCODE,38)K
  626.       DO 856 J=1,JJ
  627.   856 WRITE(IOCODE,72)(TY(I,J,K),I=1,II)
  628. 855   CONTINUE
  629. 8551  CONTINUE
  630.       IF(NUMTZ.EQ.0) GO TO 120
  631.       WRITE(IOCODE,37)
  632.    37 FORMAT(//T15,'**********TRANSMISSIBILITY (TZ) NODE MODIFICATIONS',
  633.      & '**********',//T15,
  634.      & '   I1  I2  J1  J2  K1  K2  NEW TZ VALUE')
  635.       DO 277 L=1,NUMTZ
  636.       READ(20,*) I1,I2,J1,J2,K1,K2,REGVAL
  637.       WRITE(IOCODE,32) I1,I2,J1,J2,K1,K2,REGVAL
  638.       DO 277 K=K1,K2
  639.       DO 277 J=J1,J2
  640.       DO 277 I=I1,I2
  641.       IF(VP(I,J,K).GT.0.) TZ(I,J,K)=REGVAL
  642.   277 CONTINUE
  643.   120 IF(ITCODE.NE.1) GO TO 8571
  644.       WRITE(IOCODE,48)
  645.    48 FORMAT(//T15,'**********TRANSMISSIBILITY (TZ) DISTRIBUTION'
  646.      & ,'**********'/)
  647.       DO 857 K=1,KK
  648.       WRITE(IOCODE,38)K
  649.       DO 858 J=1,JJ
  650.   858 WRITE(IOCODE,72)(TZ(I,J,K),I=1,II)
  651. 857   CONTINUE
  652. 8571  CONTINUE
  653. C      TRANSMISSIBILITY ARRAY CHECK
  654.       DO 900 K=1,KK
  655.       DO 900 J=1,JJ
  656.       DO 900 I=1,II
  657.       IF(TX(I,J,K).GE.0.0) GO TO 894
  658.       IFATAL=IFATAL+1
  659.       WRITE(IOCODE,893) I,J,K
  660.   893 FORMAT(/,5X,5('-'),'GRID BLOCK TX ERROR AT IJK =',3I5)
  661.   894 IF(TY(I,J,K).GE.0.0) GO TO 896
  662.       IFATAL=IFATAL+1
  663.       WRITE(IOCODE,895) I,J,K
  664.   895 FORMAT(/,5X,5('-'),'GRID BLOCK TY ERROR AT IJK =',3I5)
  665.   896 IF(TZ(I,J,K).GE.0.0) GO TO 900
  666.       IFATAL=IFATAL+1
  667.       WRITE(IOCODE,897) I,J,K
  668.   897 FORMAT(/,5X,5('-'),'GRID BLOCK TZ ERROR AT IJK =',3I5)
  669.   900 CONTINUE
  670.       RETURN
  671.       END
  672. C........................................................................TRIKR
  673.       SUBROUTINE TRIKRO(IREG,SO,SW,RKRO)
  674. C      MACHINE DEPENDENT INCLUDE STATEMENT
  675. $INCLUDE:'PARAMS.FOR'
  676. C      3-PHASE REL PERM CALC
  677.       REAL KROT,KROGT,KROW,KROG,KROWR,KRG,KRW,KRGT,KRWT
  678.       REAL MUOT,MUGT,MUWT
  679.       COMMON /BUBBLE/ PBO,VSLOPE(LP8),BSLOPE(LP8),RSLOPE(LP8),PMAXT,
  680.      & IREPRS,MPGT(LP8),
  681.      & RHOSCO(LP8),RHOSCW(LP8),RHOSCG(LP8),MSAT(LP7),MPOT(LP8),
  682.      & MPWT(LP8),PBOT(LP1,LP2,LP3),PBOTN(LP1,LP2,LP3)
  683.       COMMON /SARRAY/ PN(LP1,LP2,LP3),IOCODE,IDMAX,
  684.      & SON(LP1,LP2,LP3),SWN(LP1,LP2,LP3),SGN(LP1,LP2,LP3),
  685.      & A1(LP1,LP2,LP3),A2(LP1,LP2,LP3),A3(LP1,LP2,LP3),
  686.      & SUM(LP1,LP2,LP3),GAM(LP1,LP2,LP3),QS(LP1,LP2,LP3)
  687.       COMMON /SPVT/ SAT(LP7,LP9),KROT(LP7,LP9),KRWT(LP7,LP9),
  688.      & BGT(LP7,LP9),
  689.      & KRGT(LP7,LP9),ITHREE(LP7),RSOT(LP7,LP9),BWPT(LP7,LP9),
  690.      & PCOWT(LP7,LP9),PCGOT(LP7,LP9),KROGT(LP7,LP9),SWR(LP7),
  691.      & POT(LP7,LP9),MUOT(LP7,LP9),BOT(LP7,LP9),BOPT(LP7,LP9),
  692.      & RSOPT(LP7,LP9),PWT(LP7,LP9),MUWT(LP7,LP9),BWT(LP7,LP9),
  693.      & RSWT(LP7,LP9),RSWPT(LP7,LP9),PGT(LP7,LP9),MUGT(LP7,LP9),
  694.      & BGPT(LP7,LP9),CRT(LP7,LP9),IPVT(LP1,LP2,LP3),IROCK(LP1,LP2,LP3),
  695.      & NROCK,NPVT,PSIT(LP7,LP9),PRT(LP7,LP9),WOROCK(LP7),GOROCK(LP7)
  696.       SOWR=1.0-SWR(IREG)
  697.       SL=SO+SW
  698.       SG=1.-SL
  699.       IF(SG.LE.0.001) GO TO 10
  700.       CALL INTERP(IREG,SAT,KROT,MSAT(IREG),1.-SW,KROW)
  701.       CALL INTERP(IREG,SAT,KROGT,MSAT(IREG),SL,KROG)
  702.       CALL INTERP(IREG,SAT,KROT,MSAT(IREG),SOWR,KROWR)
  703.       CALL INTERP(IREG,SAT,KRWT,MSAT(IREG),SW,KRW)
  704.       CALL INTERP(IREG,SAT,KRGT,MSAT(IREG),SG,KRG)
  705.       RKRO=((KROW+KRW)*(KROG+KRG))/KROWR
  706.      & -(KRW+KRG)
  707. C     WRITE(IOCODE,5) KROWR,KROW,KROG,KRW,KRG,RKRO
  708. C 3   FORMAT(5X,'KROWR,KROW,KROG,KRW,KRG,RKRO = ',5F10.4)
  709.       IF(RKRO.GT.1.0) RKRO=1.0
  710.       IF(RKRO.LT.0.0) RKRO=0.0
  711. C     WRITE(IOCODE,15) SOWR,SL,SG,SO,KROG,RKRO
  712.       RETURN
  713. C      NO FREE GAS IN THE BLOCK
  714.    10 CALL INTERP(IREG,SAT,KROT,MSAT(IREG),SO,RKRO)
  715. C     WRITE(IOCODE,15) SOWR,SL,SG,SO,KROG,RKRO
  716. C  15 FORMAT(1X,'SOWR,SL,SG,SO,KROG,RKRO = ',6F10.4)
  717.       RETURN
  718.       END
  719. C........................................................................UINIT
  720.       SUBROUTINE UINITL(KPI,II,JJ,KK,PI,ET0,CUMPO,MBEO,
  721.      &CUMPW,MBEW,CUMPG,MBEG,SOI,SWI,SGI,
  722.      &WOC,GOC,PGOC,CUMIW,CUMIG)
  723. C      MACHINE DEPENDENT INCLUDE STATEMENT
  724. $INCLUDE:'PARAMS.FOR'
  725. C      P AND SAT INITIALIZATION
  726.       DIMENSION WOC(LP7),GOC(LP7)
  727.       REAL MBEO,MBEW,MBEG,KX,KY,KZ,KROT,KRWT,KRGT,KROGT,MUOT,MUWT,MUGT
  728.       COMMON /BUBBLE/ PBO,VSLOPE(LP8),BSLOPE(LP8),RSLOPE(LP8),PMAXT,
  729.      & IREPRS,MPGT(LP8),
  730.      & RHOSCO(LP8),RHOSCW(LP8),RHOSCG(LP8),MSAT(LP7),MPOT(LP8),
  731.      & MPWT(LP8),PBOT(LP1,LP2,LP3),PBOTN(LP1,LP2,LP3)
  732.       COMMON /SARRAY/ PN(LP1,LP2,LP3),IOCODE,IDMAX,
  733.      & SON(LP1,LP2,LP3),SWN(LP1,LP2,LP3),SGN(LP1,LP2,LP3),
  734.      & A1(LP1,LP2,LP3),A2(LP1,LP2,LP3),A3(LP1,LP2,LP3),
  735.      & SUM(LP1,LP2,LP3),GAM(LP1,LP2,LP3),QS(LP1,LP2,LP3)
  736.       COMMON /SPARM/ KX(LP1,LP2,LP3),KY(LP1,LP2,LP3),KZ(LP1,LP2,LP3),
  737.      & EL(LP1,LP2,LP3),TX(LP4,LP2,LP3),TY(LP1,LP5,LP3),TZ(LP1,LP2,LP6),
  738.      & PDAT(LP1,LP2,LP3),PDATUM,GRAD
  739.       COMMON /SPRTPS/ P(LP1,LP2,LP3),SO(LP1,LP2,LP3),SW(LP1,LP2,LP3),
  740.      & SG(LP1,LP2,LP3)
  741.       COMMON /SPVT/ SAT(LP7,LP9),KROT(LP7,LP9),KRWT(LP7,LP9),
  742.      & BGT(LP7,LP9),
  743.      & KRGT(LP7,LP9),ITHREE(LP7),RSOT(LP7,LP9),BWPT(LP7,LP9),
  744.      & PCOWT(LP7,LP9),PCGOT(LP7,LP9),KROGT(LP7,LP9),SWR(LP7),
  745.      & POT(LP7,LP9),MUOT(LP7,LP9),BOT(LP7,LP9),BOPT(LP7,LP9),
  746.      & RSOPT(LP7,LP9),PWT(LP7,LP9),MUWT(LP7,LP9),BWT(LP7,LP9),
  747.      & RSWT(LP7,LP9),RSWPT(LP7,LP9),PGT(LP7,LP9),MUGT(LP7,LP9),
  748.      & BGPT(LP7,LP9),CRT(LP7,LP9),IPVT(LP1,LP2,LP3),IROCK(LP1,LP2,LP3),
  749.      & NROCK,NPVT,PSIT(LP7,LP9),PRT(LP7,LP9),WOROCK(LP7),GOROCK(LP7)
  750.       COMMON /SSOLN/ BO(LP1,LP2,LP3),BW(LP1,LP2,LP3),BG(LP1,LP2,LP3),
  751.      & QO(LP1,LP2,LP3),QW(LP1,LP2,LP3),QG(LP1,LP2,LP3),
  752.      & GOWT(LP1,LP2,LP3),GWWT(LP1,LP2,LP3),GGWT(LP1,LP2,LP3),
  753.      & OW(LP4,LP2,LP3),OE(LP4,LP2,LP3),WW(LP4,LP2,LP3),WE(LP4,LP2,LP3),
  754.      & OS(LP1,LP5,LP3),ON(LP1,LP5,LP3),WS(LP1,LP5,LP3),WN(LP1,LP5,LP3),
  755.      & OT(LP1,LP2,LP6),OB(LP1,LP2,LP6),WT(LP1,LP2,LP6),WB(LP1,LP2,LP6),
  756.      & QOWG(LP1,LP2,LP3),VP(LP1,LP2,LP3),CT(LP1,LP2,LP3)
  757.       COMMON /TSTDAT/ IFATAL,IWARN
  758.       ET0=0.0
  759.       CUMPO=0.0
  760.       CUMPW=0.0
  761.       CUMPG=0.0
  762.       CUMIW=0.0
  763.       CUMIG=0.0
  764.       MBEO=0.0
  765.       MBEW=0.0
  766.       MBEG=0.0
  767.       READ(20,69)
  768.       READ(20,*) KPI,KSI,PDATUM,GRAD
  769.       DO 250 NR=1,NROCK
  770.       IF(KPI.NE.0) GO TO 250
  771.       READ(20,*) PI,WOC(NR),PGOC,GOC(NR)
  772.       DO 200 K=1,KK
  773.       DO 200 J=1,JJ
  774.       DO 200 I=1,II
  775.       BPT=PBOT(I,J,K)
  776.       IPVTR=IPVT(I,J,K)
  777.       IF(IROCK(I,J,K).NE.NR) GO TO 200
  778.       IF(EL(I,J,K).LT.GOC(NR)) GO TO 175
  779.       IF(EL(I,J,K).GT.WOC(NR))GO TO 150
  780.       CALL INTPVT(IPVTR,BPT,BSLOPE(IPVTR),POT,BOT,MPOT(IPVTR),PI,BBO)
  781.       CALL INTPVT(IPVTR,BPT,RSLOPE(IPVTR),POT,RSOT,MPOT(IPVTR),PI,RSO)
  782.       RHOO=(RHOSCO(IPVTR)+RSO*RHOSCG(IPVTR))/BBO
  783.       PN(I,J,K)=PI+RHOO*(EL(I,J,K)-WOC(NR))/144.
  784.       GO TO 200
  785.   150 CALL INTERP(IPVTR,PWT,BWT,MPWT(IPVTR),PI,BBW)
  786.       CALL INTERP(IPVTR,PWT,RSWT,MPWT(IPVTR),PI,RSW)
  787.       RHOW=(RHOSCW(IPVTR)+RSW*RHOSCG(IPVTR))/BBW
  788.       PN(I,J,K)=PI+RHOW*(EL(I,J,K)-WOC(NR))/144.
  789.       GO TO 200
  790.   175 CALL INTERP(IPVTR,PGT,BGT,MPGT(IPVTR),PGOC,BBG)
  791.       RHOG=RHOSCG(IPVTR)/BBG
  792.       PN(I,J,K)=PGOC+RHOG*(EL(I,J,K)-GOC(NR))/144.
  793.   200 CONTINUE
  794.   250 CONTINUE
  795.       IF(KPI.EQ.0) GO TO 9010
  796.       DO 3010 K=1,KK
  797.       DO 3005 J=1,JJ
  798.  3005 READ(20,*)(PN(I,J,K),I=1,II)
  799. 3010  CONTINUE
  800. 9010  CONTINUE
  801. C**** INITIALIZE N+1 PRESSURE ARRAY.
  802.       DO 3012 I=1,II
  803.       DO 3012 J=1,JJ
  804.       DO 3012 K=1,KK
  805.  3012 P(I,J,K)=PN(I,J,K)
  806.       IF(KSI.NE.0)GO TO 5600
  807.       DO 40 NR=1,NROCK
  808.       READ(20,*) SOI,SWI,SGI
  809.       DO 30 K=1,KK
  810.       DO 30 J=1,JJ
  811.       DO 30 I=1,II
  812.       IF(IROCK(I,J,K).NE.NR) GO TO 30
  813.       SON(I,J,K)=SOI
  814.       IF(EL(I,J,K).GT.WOC(NR)) SON(I,J,K)=0.
  815.       IF(EL(I,J,K).LT.GOC(NR)) SON(I,J,K)=0.
  816.       SWN(I,J,K)=SWI
  817.       IF(EL(I,J,K).GT.WOC(NR)) SWN(I,J,K)=1.
  818.       IF(EL(I,J,K).LT.GOC(NR)) SWN(I,J,K)=SWR(NR)
  819.       SGI=1.0-SOI-SWI
  820.       SGN(I,J,K)=SGI
  821.       SO(I,J,K)=SON(I,J,K)
  822.       SW(I,J,K)=SWN(I,J,K)
  823.       SG(I,J,K)=SGN(I,J,K)
  824.       IF(SG(I,J,K).LT.0.0) SG(I,J,K)=0.0
  825.    30 CONTINUE
  826.    40 CONTINUE
  827.       GO TO 4000
  828.  5600 DO 3011 K=1,KK
  829.       DO 3006 J=1,JJ
  830. 3006  READ(20,*)(SO(I,J,K),I=1,II)
  831. 3011  CONTINUE
  832.       DO 3020 K=1,KK
  833.       DO 3007 J=1,JJ
  834. 3007  READ(20,*)(SW(I,J,K),I=1,II)
  835. 3020  CONTINUE
  836.       DO 3030 K=1,KK
  837.       DO 3030 J=1,JJ
  838.       DO 3030 I=1,II
  839.       SG(I,J,K)=1.0-SO(I,J,K)-SW(I,J,K)
  840.       IF(SG(I,J,K).LT.0.0) SG(I,J,K)=0.0
  841.       SON(I,J,K)=SO(I,J,K)
  842.       SWN(I,J,K)=SW(I,J,K)
  843.       SGN(I,J,K)=SG(I,J,K)
  844.  3030 CONTINUE
  845. 69    FORMAT(40A2)
  846.  4000 CONTINUE
  847.       DO 4200 K=1,KK
  848.       DO 4200 J=1,JJ
  849.       DO 4200 I=1,II
  850.       IF(VP(I,J,K).GT.0.0) GO TO 4100
  851.       P(I,J,K)=0.0
  852.       SO(I,J,K)=0.0
  853.       SW(I,J,K)=0.0
  854.       SG(I,J,K)=0.0
  855.       PN(I,J,K)=0.0
  856.       SON(I,J,K)=0.0
  857.       SWN(I,J,K)=0.0
  858.       SGN(I,J,K)=0.0
  859.  4100 CONTINUE
  860. C      CHECK INITIAL ARRAYS
  861.       IF(VP(I,J,K).LE.0.) GO TO 4200
  862.       IF(P(I,J,K).GE.0.0) GO TO 4120
  863.       IFATAL=IFATAL+1
  864.       WRITE(IOCODE,4110) I,J,K
  865.  4110 FORMAT(/5X,5('-'),'INIT P ERROR AT GRID BLOCK IJK = ',3I5)
  866.  4120 IF(SO(I,J,K).GE.0.0.AND.SO(I,J,K).LE.1.0) GO TO 4140
  867.       IFATAL=IFATAL+1
  868.       WRITE(IOCODE,4130) I,J,K
  869.  4130 FORMAT(/5X,5('-'),'INIT SO ERROR AT GRID BLOCK IJK = ',3I5)
  870.  4140 IF(SW(I,J,K).GE.0.0.AND.SW(I,J,K).LE.1.0) GO TO 4160
  871.       IFATAL=IFATAL+1
  872.       WRITE(IOCODE,4150) I,J,K
  873.  4150 FORMAT(/5X,5('-'),'INIT SW ERROR AT GRID BLOCK IJK = ',3I5)
  874.  4160 IF(SG(I,J,K).GE.0.0.AND.SG(I,J,K).LE.1.0) GO TO 4180
  875.       IFATAL=IFATAL+1
  876.       WRITE(IOCODE,4170) I,J,K
  877.  4170 FORMAT(/5X,5('-'),'INIT SG ERROR AT GRID BLOCK IJK = ',3I5)
  878.  4180 SUMSAT=SO(I,J,K)+SW(I,J,K)+SG(I,J,K)
  879.       ERR=(1.0-SUMSAT)
  880.       IF(ABS(ERR).LE.0.0001) GO TO 4200
  881.        WRITE(IOCODE,10) SO(I,J,K) ,SW(I,J,K),SG(I,J,K)
  882. 10     FORMAT(/,5X,' SO=',F10.6,' SW=',F10.6,' SG=',F10.6)
  883.       IFATAL=IFATAL+1
  884.       WRITE(IOCODE,4190) I,J,K, SUMSAT
  885.  4190 FORMAT(/5X,5('-'),'INIT SAT SUM ERROR AT GRID BLOCK IJK = ',3I5,'
  886.      & SUMSAT = ',F10.5)
  887.  4200 CONTINUE
  888.       RETURN
  889.       END
  890. C......................................................VISCY
  891.       SUBROUTINE VISCY(TEMPRD,PRSPRD,SPG,TEM,CNCH2S,CNCCO2,
  892.      &  CNCN2,VISGR,VISG,IERR)
  893. C      MACHINE DEPENDENT INCLUDE STATEMENT
  894. $INCLUDE:'PARAMS.FOR'
  895. C      CORRELATION GAS VISCOSITY
  896.       DIMENSION TEMTBL(13),PRSTBL(22),VISTBL(22,13),VISTB1(22,7),
  897.      &  VISTB2(22,6)
  898.       EQUIVALENCE (VISTB1(1,1),VISTBL(1,1)),(VISTB2(1,1),VISTBL(1,8))
  899.       DATA TEMTBL/1.05,
  900.      &1.10,1.15,1.20,1.30,1.40,1.50,1.60,1.75,2.00,2.25,2.50,3.0/
  901.       DATA PRSTBL/0.1,
  902.      &0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2,
  903.      & 1.4, 1.6, 1.8, 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0, 20.0/
  904.       DATA VISTB1 /1.000
  905.      A,1.012,1.025,1.050,1.075,1.100,1.145,1.195,1.285,1.415,1.760,
  906.      B2.285,2.865,3.290,3.650,4.760,5.500,6.460,7.150,7.680,8.650,9.370,
  907.      C1.000,1.011,1.023,1.043,1.065,1.086,1.120,1.150,1.195,1.255,1.435,
  908.      D1.700,2.070,2.465,2.800,3.850,4.655,5.720,6.500,7.060,8.100,8.880,
  909.      E1.000,1.010,1.021,1.036,1.055,1.073,1.095,1.120,1.145,1.175,1.280,
  910.      F1.420,1.590,1.850,2.160,3.225,3.975,5.030,5.820,6.385,7.410,8.180,
  911.      G1.000,1.009,1.019,1.030,1.045,1.060,1.070,1.085,1.110,1.135,1.195,
  912.      H1.285,1.425,1.570,1.750,2.600,3.350,4.380,5.125,5.740,6.750,7.500,
  913.      I1.000,1.008,1.017,1.027,1.040,1.054,1.063,1.075,1.100,1.120,1.155,
  914.      J1.215,1.285,1.360,1.460,2.020,2.560,3.500,4.185,4.755,5.790,6.500,
  915.      K1.000,1.007,1.015,1.024,1.035,1.048,1.056,1.067,1.089,1.110,1.135,
  916.      L1.185,1.235,1.280,1.335,1.690,2.110,2.790,3.380,3.860,4.790,5.410,
  917.      M1.000,1.006,1.013,1.021,1.030,1.042,1.049,1.059,1.078,1.100,1.120,
  918.      N1.150,1.185,1.220,1.260,1.500,1.785,2.325,2.820,3.230,4.060,4.610/
  919.       DATA VISTB2 /1.000
  920.      A,1.005,1.011,1.018,1.025,1.036,1.042,1.051,1.067,1.070,1.095,
  921.      B1.120,1.150,1.180,1.215,1.385,1.595,2.030,2.425,2.770,3.490,4.025,
  922.      C1.000,1.004,1.009,1.015,1.021,1.030,1.035,1.043,1.056,1.065,1.090,
  923.      D1.110,1.125,1.145,1.165,1.280,1.435,1.770,2.095,2.375,2.990,3.500,
  924.      E1.000,1.003,1.007,1.012,1.017,1.024,1.028,1.035,1.045,1.055,1.060,
  925.      F1.070,1.080,1.095,1.110,1.205,1.290,1.500,1.725,1.955,2.480,2.925,
  926.      G1.000,1.002,1.005,1.009,1.013,1.018,1.021,1.027,1.034,1.040,1.045,
  927.      H1.055,1.065,1.075,1.085,1.145,1.210,1.340,1.485,1.665,2.085,2.460,
  928.      I1.000,1.001,1.003,1.006,1.009,1.012,1.015,1.019,1.023,1.025,1.030,
  929.      J1.040,1.050,1.060,1.065,1.105,1.155,1.245,1.360,1.485,1.830,2.150,
  930.      K1.000,1.000,1.001,1.003,1.005,1.007,1.009,1.011,1.013,1.015,1.020,
  931.      L1.025,1.030,1.035,1.040,1.060,1.085,1.140,1.205,1.265,1.495,1.750/
  932.       LOPRS = 2
  933.       IHIPRS = 21
  934.       LOTEM = 3
  935.       IHITEM = 11
  936.       IERR = 0
  937.       IF (TEMPRD .LT. 1.05 .OR. TEMPRD .GT.  3.00) GO TO 340
  938.       IF (PRSPRD .LT. 0.01 .OR. PRSPRD .GT. 20.00) GO TO 340
  939.       IF (SPG    .LT. 0.55 .OR. SPG    .GT.  1.50) GO TO 340
  940.       IF (TEM    .LT. 40.0 .OR. TEM    .GT. 400.0) GO TO 340
  941.       DO 100 J = LOTEM, IHITEM
  942.       IF (TEMTBL(J) .GE. TEMPRD) GO TO 120
  943.   100 CONTINUE
  944.       J      = IHITEM + 1
  945.   120 IF (PRSTBL(LOPRS - 1) .LT. PRSPRD) GO TO 160
  946.       I      = 1
  947.       ICALL  = 1
  948.       GO TO 320
  949.   140 VISGR  = 1.00 + (PRSPRD * (VISI - 1.0))
  950.       GO TO 300
  951.   160 IF(PRSTBL(IHIPRS+1).GT.PRSPRD) GO TO 200
  952.       I=IHIPRS+1
  953.       ICALL=2
  954.       GO TO 320
  955.   180 VISGR=VISI
  956.       GO TO 300
  957.   200 DO 220 I=LOPRS,IHIPRS
  958.       IF(PRSTBL(I).GE.PRSPRD) GO TO 240
  959.   220 CONTINUE
  960.       I=IHIPRS+1
  961.   240 ICALL=3
  962.       GO TO 320
  963.   260 VISJ=VISI
  964.       I=I-1
  965.       ICALL=4
  966.       GO TO 320
  967.   280 I=I+1
  968.       VISGR=VISI + (((PRSPRD-PRSTBL(I-1))/(PRSTBL(I)-PRSTBL(I-1)))
  969.      &  *(VISJ-VISI))
  970.   300 VISGU=(0.126585E-01) - (0.611823E-02)*SPG + (0.164574E-02)*SPG
  971.      & *SPG + (0.164574E-04)*TEM - (0.719221E-06)*SPG*TEM
  972.      & - (0.609046E-06)*SPG*SPG*TEM
  973.       CORH2S=(0.000113*CNCH2S*SPG - 0.000038*CNCH2S + 0.000001)
  974.      & *(1./(1.+SPG)) + 0.000001
  975.       CORCO2=(0.000134*CNCCO2*SPG - 0.000004*CNCCO2 + 0.000004*SPG)
  976.      & *(1./(1.+SPG))-0.000003
  977.       CORN2=(0.000170*CNCN2*SPG + 0.000021*CNCN2 + 0.000010*SPG)
  978.      & *(1./(1.+SPG)) - 0.000006
  979.       VISGA=VISGU + CORH2S + CORCO2 + CORN2
  980.       VISG=VISGR*VISGA
  981.       GO TO 360
  982.   320 CALL XLGR4(TEMPRD,TEMTBL(J-2),TEMTBL(J-1),TEMTBL(J),
  983.      &  TEMTBL(J+1),VISI,VISTBL(I,J-2),VISTBL(I,J-1),
  984.      &  VISTBL(I,J),VISTBL(I,J+1))
  985.       GO TO (140,180,260,280),ICALL
  986.   340 IERR=1
  987.       VISGR=0.0
  988.       VISG=0.0
  989.   360 RETURN
  990.       END
  991. C......................................................XLGR4
  992.       SUBROUTINE XLGR4(X,X1,X2,X3,X4,Y,Y1,Y2,Y3,Y4)
  993. C      MACHINE DEPENDENT INCLUDE STATEMENT
  994. $INCLUDE:'PARAMS.FOR'
  995. C      LAGRANGIAN INTERP
  996.       A1=X1-X2
  997.       A2=X1-X3
  998.       A3=X1-X4
  999.       A4=X2-X3
  1000.       A5=X2-X4
  1001.       A6=X3-X4
  1002.       B1=X-X1
  1003.       B2=X-X2
  1004.       B3=X-X3
  1005.       B4=X-X4
  1006.       Y= B2/A1 *B3/A2 *B4/A3* Y1 -
  1007.      &   B1/A1 *B3/A4 *B4/A5* Y2 +
  1008.      &   B1/A2 *B2/A4 *B4/A6 *Y3 -
  1009.      &   B1/A3 *B2/A5 *B3/A6* Y4
  1010.       RETURN
  1011.       END
  1012. C......................................................ZANDC
  1013.       SUBROUTINE ZANDC(TEMP,TEMPC,PRS,PRSPC,CNCH2S,
  1014.      & CNCCO2,ZED,CMPG,IERR)
  1015. C      MACHINE DEPENDENT INCLUDE STATEMENT
  1016. $INCLUDE:'PARAMS.FOR'
  1017. C      Z-FACTOR AND GAS COMP CALC USING CORRELATIONS
  1018.       DIMENSION A(8)
  1019.       DATA A/0.31506237,
  1020.      &-1.04670990, -0.57832729, 0.53530771,
  1021.      &-0.61232032, -0.10488813, 0.68157001, 0.68446549/
  1022.       FRCA=(CNCH2S+CNCCO2)/100.
  1023.       FRCB=CNCH2S/100.
  1024.       EPS=120.*(FRCA**0.9 -FRCA**1.6) + 15.*(FRCB**0.5
  1025.      &  -FRCB**4.0)
  1026.       TEMPCA=TEMPC-EPS
  1027.       PRSPCA=PRSPC*TEMPCA/(TEMPC+FRCB*(1.-FRCB)*EPS)
  1028.       TEMPRD=(TEMP+460.)/TEMPCA
  1029.       PRSPRD=PRS/PRSPCA
  1030.       IERR=0
  1031.       IF(TEMPRD.LT.1.05.OR.TEMPRD.GT.3.0) GO TO 140
  1032.       IF(PRSPRD.LT.0.0.OR.PRSPRD.GT.15.0) GO TO 140
  1033.       IF(FRCA.LT.0.0.OR.FRCA.GT.0.85) GO TO 140
  1034.       ITER=0
  1035.       T1=A(1)*TEMPRD +A(2) +A(3)/(TEMPRD*TEMPRD)
  1036.       T2=A(4)*TEMPRD + A(5)
  1037.       T3=A(5)*A(6)
  1038.       T4=A(7)/(TEMPRD*TEMPRD)
  1039.       T5=A(8)
  1040.       DENRD=1.0
  1041.       DO 120 ITER=1,10
  1042.       DENRD2=DENRD*DENRD
  1043.       DENRD3=DENRD2*DENRD
  1044.       DENRD4=DENRD2*DENRD2
  1045.       DENRD5=DENRD2*DENRD3
  1046.       P=(TEMPRD +T1*DENRD + T2*DENRD2 + T3*DENRD5)
  1047.      &  *DENRD + T4*DENRD3*(1. + T5*DENRD2)*
  1048.      & EXP(-T5*DENRD2)
  1049.       DP=TEMPRD + 2.*T1*DENRD + 3.*T2*DENRD2
  1050.      & + 6.*T3*DENRD5 + T4*DENRD2*EXP
  1051.      &(-T5*DENRD2)*(3. + 3.*T5*DENRD2
  1052.      & - 2.*T5*T5*DENRD4)
  1053.       DENRD1=DENRD-(P-0.270*PRSPRD)/DP
  1054.       IF(DENRD1.GT.0.) GO TO 100
  1055.       DENRD1=0.5*DENRD
  1056.   100 IF(DENRD1.LT.2.2) GO TO 110
  1057.       DENRD1=DENRD + 0.9*(2.2-DENRD)
  1058.   110 IF(ABS(DENRD-DENRD1).LT.0.00001) GO TO 130
  1059.   120 DENRD=DENRD1
  1060.   130 ZED=0.270*PRSPRD/(DENRD1*TEMPRD)
  1061.       DENRD=DENRD1
  1062.       DENRD2=DENRD*DENRD
  1063.       DENRD4=DENRD2*DENRD2
  1064.       DZED=T1/TEMPRD + 2.*T2/TEMPRD*DENRD
  1065.      & + 5.*T3*DENRD4/TEMPRD + (1.+T5*DENRD2
  1066.      &  - T5*T5*DENRD4) * 2.*T4/TEMPRD*DENRD
  1067.      & *EXP(-T5*DENRD2)
  1068.       CMPPRD=1./PRSPRD -0.270*DZED/(ZED*ZED*TEMPRD
  1069.      & * (1.+DENRD/ZED*DZED))
  1070.       CMPG=CMPPRD/PRSPCA
  1071.       GO TO 150
  1072.   140 IERR=1
  1073.       ZED=0.
  1074.       CMPG=0.
  1075.   150 RETURN
  1076.       END
  1077. C.........................................................FPTD
  1078.       FUNCTION FPTD(A0,A1,A2,A3,X)
  1079. C      CALC. DIM PRESSURE AT TIME N+1 FOR CARTER-TRACY
  1080.       FPTD=A0+A1*X+A2*ALOG(X)+A3*ALOG(X)*ALOG(X)
  1081.       RETURN
  1082.       END
  1083. C..........................................................FDPTD
  1084.       FUNCTION FDPTD(A1,A2,A3,X)
  1085. C      CALC. DERIVATIVE OF DIM P WRT DIM TIME AT N+1
  1086. C      FOR CARTER-TRACY
  1087.       FDPTD=A1+(A2/X)+(2.*A3*ALOG(X)/X)
  1088.       RETURN
  1089.       END
  1090.