home *** CD-ROM | disk | FTP | other *** search
/ Hall of Fame / HallofFameCDROM.cdr / oilfield / spe-46-2.lzh / BLOCK2.FOR < prev    next >
Text File  |  1988-07-20  |  30KB  |  889 lines
  1. $DO66
  2. C.................................................................CODES
  3.       SUBROUTINE CODES(IOCODE,KSM1,KSN1,KCO1,NN,
  4.      &FACT1,FACT2,TMAX,KSOL,MITER,OMEGA,TOL,TOL1,
  5.      & KSN,KSM,KCO,KCOFF,DSMAX,DPMAX,
  6.      & NUMDIS,IRK,THRUIN,
  7.      &WORMAX,GORMAX,PAMIN,PAMAX)
  8. C      MACHINE DEPENDENT INCLUDE STATEMENT
  9. $INCLUDE:'PARAMS.FOR'
  10.       REAL KROT,KRWT,KRGT,KROGT,MUOT,MUWT,MUGT
  11.       COMMON /SPVT/ SAT(LP7,LP9),KROT(LP7,LP9),KRWT(LP7,LP9),
  12.      & BGT(LP7,LP9),
  13.      & KRGT(LP7,LP9),ITHREE(LP7),RSOT(LP7,LP9),BWPT(LP7,LP9),
  14.      & PCOWT(LP7,LP9),PCGOT(LP7,LP9),KROGT(LP7,LP9),SWR(LP7),
  15.      & POT(LP7,LP9),MUOT(LP7,LP9),BOT(LP7,LP9),BOPT(LP7,LP9),
  16.      & RSOPT(LP7,LP9),PWT(LP7,LP9),MUWT(LP7,LP9),BWT(LP7,LP9),
  17.      & RSWT(LP7,LP9),RSWPT(LP7,LP9),PGT(LP7,LP9),MUGT(LP7,LP9),
  18.      & BGPT(LP7,LP9),CRT(LP7,LP9),IPVT(LP1,LP2,LP3),IROCK(LP1,LP2,LP3),
  19.      & NROCK,NPVT,PSIT(LP7,LP9),PRT(LP7,LP9),WOROCK(LP7),GOROCK(LP7)
  20. C      READ DEBUG, MATRIX SOLUTION AND RUN CONTROL INFO
  21.       COMMON /TSTDAT/ IFATAL,IWARN
  22.       READ(20,69)
  23.       READ(20,*) KSN1,KSM1,KCO1,KCOFF
  24. C EVERY KSM1 TH STEP SOLUTION MATRIX WILL BE WRITTEN
  25. C EVERY KSN1 TH STEP LSORX DATA WILL BE WRITTEN
  26. C EVERY KCO1 TH STEP COMPRESSIBILITIES & VOLUME FACTORS WILL BE WRITTEN
  27.       READ(20,69)
  28.       READ(20,*)NN,FACT1,FACT2,TMAX,WORMAX,GORMAX,PAMIN,PAMAX
  29.       WRITE(IOCODE,59)NN,FACT1,FACT2,TMAX,WORMAX,GORMAX,PAMIN,PAMAX
  30.       IF(WORMAX.NE.0.0) GO TO 20
  31.       READ(20,*) (WOROCK(I),I=1,NROCK)
  32.       WRITE(IOCODE,5)
  33.     5 FORMAT(//T15,'ROCK REGION SPECIFIED WOR MAXIMA:')
  34.       DO 10 I=1,NROCK
  35.       WRITE(IOCODE,8) I,WOROCK(I)
  36.     8 FORMAT(T20,'ROCK REGION ',I3,' WOR MAX (STB/STB) = ',F8.0)
  37.    10 CONTINUE
  38.    20 IF(GORMAX.NE.0.0) GO TO 40
  39.       READ(20,*) (GOROCK(I),I=1,NROCK)
  40.       WRITE(IOCODE,25)
  41.    25 FORMAT(//T15,'ROCK REGION SPECIFIED GOR MAXIMA:')
  42.       DO 30 I=1,NROCK
  43.       WRITE(IOCODE,28) I,GOROCK(I)
  44.    28 FORMAT(T20,'ROCK REGION ',I3,' GOR MAX (SCF/STB) = ',F8.0)
  45.    30 CONTINUE
  46.    40 CONTINUE
  47. C NN--MAXIMUM NUMBER OF TIME-STEPS
  48. C FACT1--FACTOR FOR INCREASING TIME-STEP SIZE
  49. C FACT2--FACTOR FOR DECREASING TIME-STEP SIZE
  50. C TMAX--MAXIMUM SIMULATION TIME
  51.       READ(20,69)
  52.       READ(20,*)KSOL,MITER,OMEGA,TOL,TOL1,DSMAX,DPMAX
  53.       IF(KSOL.EQ.1) WRITE(IOCODE,73)
  54.       IF(KSOL.EQ.2) WRITE(IOCODE,75) MITER,OMEGA,TOL,TOL1
  55.       IF(KSOL.EQ.3) WRITE(IOCODE,76) MITER,OMEGA,TOL,TOL1
  56.       IF(KSOL.EQ.4) WRITE(IOCODE,77) MITER,OMEGA,TOL,TOL1
  57.       IF(KSOL.GT.4) THEN
  58.         WRITE(IOCODE,*) 'IMPROPER SOLUTION METHOD SPECIFIED,KSOL=',KSOL
  59.         IFATAL=IFATAL + 1
  60.       ENDIF
  61.       WRITE(IOCODE,99) DSMAX,DPMAX
  62.    73 FORMAT(//T15,'SOLUTION METHOD IS BAND.')
  63.    75 FORMAT(//T15,'SOLUTION METHOD IS LSORX:',
  64.      &/T20,'MAXIMUM NUMBER OF ITERATIONS      (MITER) = ',5X,I5,
  65.      &/T20,'INITIAL ACCELERATION PARAMETER    (OMEGA) = ',F10.4,
  66.      &/T20,'MAXIMUM PRESSURE RESIDUAL           (TOL) = ',F10.4,
  67.      &/T20,'PARAMETER FOR CHANGING OMEGA       (TOL1) = ',F10.4)
  68.    76 FORMAT(//T15,'SOLUTION METHOD IS LSORY:',
  69.      &/T20,'MAXIMUM NUMBER OF ITERATIONS      (MITER) = ',5X,I5,
  70.      &/T20,'INITIAL ACCELERATION PARAMETER    (OMEGA) = ',F10.4,
  71.      &/T20,'MAXIMUM PRESSURE RESIDUAL           (TOL) = ',F10.4,
  72.      &/T20,'PARAMETER FOR CHANGING OMEGA       (TOL1) = ',F10.4)
  73.    77 FORMAT(//T15,'SOLUTION METHOD IS LSORZ:',
  74.      &/T20,'MAXIMUM NUMBER OF ITERATIONS      (MITER) = ',5X,I5,
  75.      &/T20,'INITIAL ACCELERATION PARAMETER    (OMEGA) = ',F10.4,
  76.      &/T20,'MAXIMUM PRESSURE RESIDUAL           (TOL) = ',F10.4,
  77.      &/T20,'PARAMETER FOR CHANGING OMEGA       (TOL1) = ',F10.4)
  78.    99 FORMAT(/T15,'AUTOMATIC TIME STEP CRITERIA:',
  79.      &/T20,'MAXIMUM ALLOWED SATURATION CHANGE (DSMAX) = ',F10.4,
  80.      &/T20,'MAXIMUM ALLOWED PRESSURE CHANGE   (DPMAX) = ',F10.4/)
  81.    59 FORMAT(1H1/T15,'RUN CONTROL PARAMETERS:',/,
  82.      & T20,'MAXIMUM NUMBER OF TIME-STEPS =',I5/
  83.      &T20,'FACTOR FOR INCREASING DELT =',F10.3,3X,
  84.      &'WHEN DSMAX AND DPMAX NOT EXCEEDED.',/,
  85.      &T20,'FACTOR FOR DECREASING DELT =',F10.3,3X,
  86.      &'WHEN DSMAX OR DPMAX IS EXCEEDED.',/,
  87.      &T20,'MAXIMUM SIMULATION TIME =',F11.3/,
  88.      &T20,'MAXIMUM RESERVOIR WOR/TIME-STEP =',F8.0,' STB/STB'/
  89.      &T20,'MAXIMUM RESERVOIR GOR/TIME-STEP =',F8.0,' SCF/STB'/
  90.      &T20,'MINIMUM AVERAGE RESERVOIR PRESSURE/TIME-STEP =',F8.0/
  91.      &T20,'MAXIMUM AVERAGE RESERVOIR PRESSURE/TIME-STEP =',F8.0)
  92.       KSN=KSN1
  93.       KSM=KSM1
  94.       KCO=KCO1
  95. 69    FORMAT(40A2)
  96. C      NUM DIS AND IMPES/R-K SOLN CONTROLS
  97.       READ(20,69)
  98.       READ(20,*) NUMDIS,IRK,THRUIN
  99.       IF(IRK.EQ.0.AND.NUMDIS.EQ.0) WRITE(IOCODE,200)
  100.   200 FORMAT(//T15,'IMPES FORMULATION SELECTED; ',
  101.      & ' SINGLE-POINT UPSTREAM WEIGHTING.')
  102.       IF(IRK.EQ.0.AND.NUMDIS.EQ.1) WRITE(IOCODE,220)
  103.   220 FORMAT(//T15,'IMPES FORMULATION SELECTED; ',
  104.      & ' TWO-POINT UPSTREAM WEIGHTING.')
  105.       IF(IRK.GT.0.AND.NUMDIS.EQ.0) WRITE(IOCODE,240)
  106.      & THRUIN
  107.   240 FORMAT(//T15,'STABILISED IMPES FORMULATION:',
  108.      & /T20,'SINGLE-POINT UPSTREAM WEIGHTING',
  109.      & /T20,'USER SPEC THROUGHPUT, FRACTION  (THRUIN) =',5X,F10.4,/)
  110.       IF(IRK.GT.0.AND.NUMDIS.EQ.1) WRITE(IOCODE,260)
  111.      & THRUIN
  112.   260 FORMAT(//T15,'STABILISED IMPES FORMULATION:',
  113.      & /T20,'TWO-POINT UPSTREAM WEIGHTING',
  114.      & /T20,'USER SPEC THROUGHPUT, FRACTION  (THRUIN) =',5X,F10.4,/)
  115. C      CHECK CODE DATA
  116. C      DEBUG
  117.       IF(KSN1.EQ.0) GO TO 1020
  118.       IWARN=IWARN+1
  119.       WRITE(IOCODE,1010)
  120.  1010 FORMAT(/5X,5('-'),'LSOR DEBUG OUTPUT ON')
  121.  1020 IF(KSM1.EQ.0) GO TO 1040
  122.       IWARN=IWARN+1
  123.       WRITE(IOCODE,1030)
  124.  1030 FORMAT(/5X,5('-'),'SOLN METHOD DEBUG OUTPUT ON')
  125.  1040 IF(KCO1.EQ.0) GO TO 1060
  126.       IWARN=IWARN+1
  127.       WRITE(IOCODE,1050)
  128.  1050 FORMAT(/5X,5('-'),'COMP AND FVF DEBUG OUTPUT ON')
  129.  1060 IF(KCOFF.EQ.0) GO TO 1400
  130.       IWARN=IWARN+1
  131.       WRITE(IOCODE,1070)
  132.  1070 FORMAT(/5X,5('-'),'DEN AND SAT DEBUG OUTPUT ON')
  133. C      RUN CONTROL
  134.  1400 CONTINUE
  135.       IF(NN.GE.1) GO TO 1420
  136.       IFATAL=IFATAL+1
  137.       WRITE(IOCODE,1410)
  138.  1410 FORMAT(/5X,5('-'),'MAX # OF TIME STEPS ERROR')
  139.  1420 IF(FACT1.GE.1.0) GO TO 1440
  140.       IFATAL=IFATAL+1
  141.       WRITE(IOCODE,1430)
  142.  1430 FORMAT(/5X,5('-'),'FACT1 ERROR')
  143.  1440 IF(FACT2.GT.0.0.AND.FACT2.LE.1.0) GO TO 1460
  144.       IFATAL=IFATAL+1
  145.       WRITE(IOCODE,1450)
  146.  1450 FORMAT(/5X,5('-'),'FACT2 ERROR')
  147.  1460 IF(TMAX.GT.0.0) GO TO 1480
  148.       IFATAL=IFATAL+1
  149.       WRITE(IOCODE,1470)
  150.  1470 FORMAT(/5X,5('-'),'TMAX ERROR')
  151.  1480 IF(WORMAX.GE.0.0) GO TO 1500
  152.       IFATAL=IFATAL+1
  153.       WRITE(IOCODE,1490)
  154.  1490 FORMAT(/5X,5('-'),'WORMAX ERROR')
  155.  1500 IF(GORMAX.GE.0.0) GO TO 1520
  156.       IFATAL=IFATAL+1
  157.       WRITE(IOCODE,1510)
  158.  1510 FORMAT(/5X,5('-'),'GORMAX ERROR')
  159.  1520 IF(PAMIN.GT.0.0) GO TO 1540
  160.       IFATAL=IFATAL+1
  161.       WRITE(IOCODE,1530)
  162.  1530 FORMAT(/5X,5('-'),'PAMIN ERROR')
  163.  1540 IF(PAMAX.GT.PAMIN) GO TO 2000
  164.       IFATAL=IFATAL+1
  165.       WRITE(IOCODE,1550)
  166.  1550 FORMAT(/5X,5('-'),'PAMAX ERROR')
  167. C      SOLUTION METHOD
  168.  2000 CONTINUE
  169.       IF(MITER.GE.1) GO TO 2020
  170.       IWARN=IWARN+1
  171.       WRITE(IOCODE,2010)
  172.  2010 FORMAT(/5X,5('-'),'ALLOWED # OF LSOR ITER LESS THAN 1')
  173.  2020 IF(OMEGA.GE.1.0.AND.OMEGA.LE.2.0) GO TO 2040
  174.       IWARN=IWARN+1
  175.       WRITE(IOCODE,2030)
  176.  2030 FORMAT(/5X,5('-'),'OMEGA LT 1 OR OMEGA GT 2')
  177.  2040 IF(TOL.GT.0) GO TO 2048
  178.       IWARN=IWARN+1
  179.       WRITE(IOCODE,2046)
  180.  2046 FORMAT(/5X,5('-'),'TOL LESS THAN OR EGUAL TO 0')
  181.  2048 IF(TOL1.GE.0.0) GO TO 2060
  182.       IWARN=IWARN+1
  183.       WRITE(IOCODE,2050)
  184.  2050 FORMAT(/5X,5('-'),'TOL1 IS NEGATIVE')
  185.  2060 IF(DSMAX.GT.0.0.AND.DSMAX.LE.1.0) GO TO 2080
  186.       IFATAL=IFATAL+1
  187.       WRITE(IOCODE,2070)
  188.  2070 FORMAT(/5X,5('-'),'DSMAX ERROR')
  189.  2080 IF(DPMAX.GT.0.0) GO TO 3000
  190.       IFATAL=IFATAL+1
  191.       WRITE(IOCODE,2090)
  192.  2090 FORMAT(/5X,5('-'),'DPMAX ERROR')
  193.  3000 CONTINUE
  194.       RETURN
  195.       END
  196. C................................................................GRIDSZ
  197.       SUBROUTINE GRIDSZ(IOCODE,II,JJ,KK)
  198. C      MACHINE DEPENDENT INCLUDE STATEMENT
  199. $INCLUDE:'PARAMS.FOR'
  200. C      READ GRID DESCRIPTION
  201.       REAL KX,KY,KZ
  202.       COMMON /TSTDAT/ IFATAL,IWARN
  203.       COMMON /SPARM/ KX(LP1,LP2,LP3),KY(LP1,LP2,LP3),KZ(LP1,LP2,LP3),
  204.      & EL(LP1,LP2,LP3),TX(LP4,LP2,LP3),TY(LP1,LP5,LP3),TZ(LP1,LP2,LP6),
  205.      & PDAT(LP1,LP2,LP3),PDATUM,GRAD
  206.       COMMON /VECTOR/ DX(LP1,LP2,LP3),DY(LP1,LP2,LP3),DZ(LP1,LP2,LP3),
  207.      & DZNET(LP1,LP2,LP3),IQN1(LP11),IQN2(LP11),IQN3(LP11),IHEDIN(80)
  208.       DIMENSION SUM(LP1,LP2),VAREL(LP1,LP2),RDXL(LP1),RDYL(LP2),
  209.      & RDZL(LP3),RKEL(LP3),RDZNET(LP3)
  210.       READ(20,69)
  211.       WRITE(IOCODE,70)
  212.    70 FORMAT(//T45,'***** INITIALIZATION DATA *****',//)
  213.       READ(20,*) II,JJ,KK
  214.       READ(20,69)
  215. C*****READ INPUT CODES FOR DX,DY,DZ
  216.       READ(20,*)KDX,KDY,KDZ,KDZNET
  217. C*****ESTABLISH GRID BLOCK LENGTH (DX) DISTRIBUTION
  218.       IF(KDX.GE.0)GO TO 180
  219.       READ(20,*)DXC
  220.       DO 175 K=1,KK
  221.       DO 175 J=1,JJ
  222.       DO 175 I=1,II
  223. 175   DX(I,J,K)=DXC
  224.       WRITE(IOCODE,56)
  225.       WRITE(IOCODE,29)DXC
  226.       GO TO 195
  227.   180 IF(KDX.GT.0)GO TO 185
  228.       READ(20,*)(RDXL(I),I=1,II)
  229.       DO 187 K=1,KK
  230.       DO 187 J=1,JJ
  231.       DO 187 I=1,II
  232.   187 DX(I,J,K)=RDXL(I)
  233.       DO 182 I=1,II
  234.   182 WRITE(IOCODE,511)I,RDXL(I)
  235.       GO TO 195
  236. 185   WRITE(IOCODE,43)
  237.       K=1
  238.       WRITE(IOCODE,38)K
  239.       DO 190 J=1,JJ
  240.       READ(20,*)(DX(I,J,K),I=1,II)
  241.   190 WRITE(IOCODE,72)(DX(I,J,K),I=1,II)
  242.       DO 194 K=2,KK
  243.       WRITE(IOCODE,38) K
  244.       DO 194 J=1,JJ
  245.       DO 193 I=1,II
  246.   193 DX(I,J,K)=DX(I,J,1)
  247.   194 WRITE(IOCODE,72) (DX(I,J,K),I=1,II)
  248. 195   CONTINUE
  249.       WRITE(IOCODE,56)
  250. C*****ESTABLISH GRID BLOCK LENGTH (DY) DISTRIBUTION
  251.       IF(KDY.GE.0)GO TO 200
  252.       READ(20,*)DYC
  253.       DO 202 K=1,KK
  254.       DO 202 J=1,JJ
  255.       DO 202 I=1,II
  256. 202   DY(I,J,K)=DYC
  257.       WRITE(IOCODE,56)
  258.       WRITE(IOCODE,33)DYC
  259.       GO TO 220
  260.   200 IF(KDY.GT.0)GO TO 207
  261.       READ(20,*)(RDYL(J),J=1,JJ)
  262.       DO 205 K=1,KK
  263.       DO 205 J=1,JJ
  264.       DO 205 I=1,II
  265.   205 DY(I,J,K)=RDYL(J)
  266.       DO 210 J=1,JJ
  267.   210 WRITE(IOCODE,512)J,RDYL(J)
  268.       GO TO 220
  269. 207   WRITE(IOCODE,47)
  270.       K=1
  271.       WRITE(IOCODE,38)K
  272.       DO 215 J=1,JJ
  273.       READ(20,*)(DY(I,J,K),I=1,II)
  274.   215 WRITE(IOCODE,72)(DY(I,J,K),I=1,II)
  275.       DO 214 K=2,KK
  276.        WRITE(IOCODE,38) K
  277.       DO 214 J=1,JJ
  278.       DO 213 I=1,II
  279.   213 DY(I,J,K)=DY(I,J,1)
  280.   214 WRITE(IOCODE,72) (DY(I,J,K),I=1,II)
  281. 220   CONTINUE
  282.       WRITE(IOCODE,56)
  283. C*****ESTABLISH GRID BLOCK LENGTH (DZ) DISTRIBUTION
  284.       IF(KDZ.GE.0)GO TO 225
  285.       READ(20,*)DZC
  286.       DO 230 K=1,KK
  287.       DO 230 J=1,JJ
  288.       DO 230 I=1,II
  289. 230   DZ(I,J,K)=DZC
  290.       WRITE(IOCODE,56)
  291.       WRITE(IOCODE,36)DZC
  292.       GO TO 245
  293.   225 IF(KDZ.GT.0)GO TO 232
  294.       READ(20,*)(RDZL(K),K=1,KK)
  295.       DO 235 K=1,KK
  296.       DO 235 J=1,JJ
  297.       DO 235 I=1,II
  298.   235 DZ(I,J,K)=RDZL(K)
  299.       DO 237 K=1,KK
  300.   237 WRITE(IOCODE,513)K,RDZL(K)
  301.       GO TO 245
  302. 232   WRITE(IOCODE,48)
  303.       DO 240 K=1,KK
  304.       WRITE(IOCODE,38)K
  305.       DO 242 J=1,JJ
  306.       READ(20,*)(DZ(I,J,K),I=1,II)
  307.   242 WRITE(IOCODE,72)(DZ(I,J,K),I=1,II)
  308. 240   CONTINUE
  309. 245   CONTINUE
  310.       WRITE(IOCODE,56)
  311.       IF(KDZNET.GE.0) GO TO 255
  312.       READ(20,*) DZNETC
  313.       DO 250 K=1,KK
  314.       DO 250 J=1,JJ
  315.       DO 250 I=1,II
  316.   250 DZNET(I,J,K)=DZNETC
  317.       WRITE(IOCODE,56)
  318.       WRITE(IOCODE,252) DZNETC
  319.   252 FORMAT(T15,'GRID BLOCK NET THICKNESS (DZNET) IS INITIALLY',
  320.      & ' SET AT',F10.4,' FOR ALL NODES',//)
  321.       GO TO 274
  322.   255 IF(KDZNET.GT.0) GO TO 268
  323.       READ(20,*) (RDZNET(K),K=1,KK)
  324.       DO 260 K=1,KK
  325.       DO 260 J=1,JJ
  326.       DO 260 I=1,II
  327.   260 DZNET(I,J,K)=RDZNET(K)
  328.       DO 266 K=1,KK
  329.   266 WRITE(IOCODE,267) K,RDZNET(K)
  330.   267 FORMAT(T15,'GRID SIZE (DZNET) IN LAYER ',I5,
  331.      & ' IS INITIALLY SET AT',F8.2,' FOR ALL NODES',/)
  332.       GO TO 274
  333.   268 WRITE(IOCODE,269)
  334.   269 FORMAT(//T15,8('*'),'GRID BLOCK NET THICKNESS ',
  335.      & '(DZNET) DISTRIBUTION',8('*')/)
  336.       DO 273 K=1,KK
  337.        WRITE(IOCODE,38) K
  338.       DO 270 J=1,JJ
  339.       READ(20,*) (DZNET(I,J,K),I=1,II)
  340.   270 WRITE(IOCODE,72) (DZNET(I,J,K),I=1,II)
  341.   273 CONTINUE
  342.   274 CONTINUE
  343.       WRITE(IOCODE,56)
  344. C********GRID BLOCK LENGTH MODIFICATIONS
  345.       READ(20,69)
  346.       READ(20,*) NUMDX,NUMDY,NUMDZ,NUMDZN,IDCODE
  347.       IF(NUMDX.EQ.0) GO TO 8531
  348.       WRITE(IOCODE,31)
  349.       DO 275 L=1,NUMDX
  350.       READ(20,*) I1,I2,J1,J2,K1,K2,REGVAL
  351.       WRITE(IOCODE,32) I1,I2,J1,J2,K1,K2,REGVAL
  352.       DO 275 K=K1,K2
  353.       DO 275 J=J1,J2
  354.       DO 275 I=I1,I2
  355.       DX(I,J,K)=REGVAL
  356.   275 CONTINUE
  357.       IF(IDCODE.NE.1)GO TO 8531
  358.       WRITE(IOCODE,43)
  359.       DO 853 K=1,KK
  360.        WRITE(IOCODE,38)K
  361.       DO 854 J=1,JJ
  362.   854 WRITE(IOCODE,72)(DX(I,J,K),I=1,II)
  363. 853   CONTINUE
  364. 8531  CONTINUE
  365.       IF(NUMDY.EQ.0) GO TO 8551
  366.       WRITE(IOCODE,34)
  367.       DO 276 L=1,NUMDY
  368.       READ(20,*) I1,I2,J1,J2,K1,K2,REGVAL
  369.       WRITE(IOCODE,32) I1,I2,J1,J2,K1,K2,REGVAL
  370.       DO 276 K=K1,K2
  371.       DO 276 J=J1,J2
  372.       DO 276 I=I1,I2
  373.       DY(I,J,K)=REGVAL
  374.   276 CONTINUE
  375.       IF(IDCODE.NE.1) GO TO 8551
  376.       WRITE(IOCODE,47)
  377.       DO 855 K=1,KK
  378.        WRITE(IOCODE,38)K
  379.       DO 856 J=1,JJ
  380.   856 WRITE(IOCODE,72)(DY(I,J,K),I=1,II)
  381. 855   CONTINUE
  382. 8551  CONTINUE
  383.       IF(NUMDZ.EQ.0) GO TO 8571
  384.       WRITE(IOCODE,37)
  385.       DO 277 L=1,NUMDZ
  386.       READ(20,*) I1,I2,J1,J2,K1,K2,REGVAL
  387.       WRITE(IOCODE,32) I1,I2,J1,J2,K1,K2,REGVAL
  388.       DO 277 K=K1,K2
  389.       DO 277 J=J1,J2
  390.       DO 277 I=I1,I2
  391.       DZ(I,J,K)=REGVAL
  392.   277 CONTINUE
  393.       IF(IDCODE.NE.1) GO TO 8571
  394.       WRITE(IOCODE,48)
  395.       DO 857 K=1,KK
  396.        WRITE(IOCODE,38)K
  397.       DO 858 J=1,JJ
  398.   858 WRITE(IOCODE,72)(DZ(I,J,K),I=1,II)
  399. 857   CONTINUE
  400. 8571  CONTINUE
  401.       IF(NUMDZN.EQ.0) GO TO 890
  402.       WRITE(IOCODE,860)
  403.   860 FORMAT(//T15,8('*'),'GRID BLOCK NET THICKNESS',
  404.      & ' (DZNET) NODE MODIFICATIONS',8('*'),//T15,
  405.      & '   I1  I2  J1  J2  K1  K2  NEW DZNET VALUE')
  406.       DO 865 L=1,NUMDZN
  407.       READ(20,*) I1,I2,J1,J2,K1,K2,REGVAL
  408.       WRITE(IOCODE,32) I1,I2,J1,J2,K1,K2,REGVAL
  409.       DO 865 K=K1,K2
  410.       DO 865 J=J1,J2
  411.       DO 865 I=I1,I2
  412.       DZNET(I,J,K)=REGVAL
  413.   865 CONTINUE
  414.       IF(IDCODE.NE.1) GO TO 890
  415.       WRITE(IOCODE,269)
  416.       DO 885 K=1,KK
  417.        WRITE(IOCODE,38) K
  418.       DO 880 J=1,JJ
  419.   880 WRITE(IOCODE,72) (DZNET(I,J,K),I=1,II)
  420.   885 CONTINUE
  421.   890 CONTINUE
  422. C      GRID BLOCK SIZE CHECK
  423.       DO 900 K=1,KK
  424.       DO 900 J=1,JJ
  425.       DO 900 I=1,II
  426.       IF(DX(I,J,K).GE.0.0) GO TO 892
  427.       IFATAL=IFATAL+1
  428.       WRITE(IOCODE,891) I,J,K
  429.   891 FORMAT(/,5X,5('-'),'GRID BLOCK DX ERROR AT IJK =',3I5)
  430.   892 IF(DY(I,J,K).GE.0.0) GO TO 894
  431.       IFATAL=IFATAL+1
  432.       WRITE(IOCODE,893) I,J,K
  433.   893 FORMAT(/,5X,5('-'),'GRID BLOCK DY ERROR AT IJK =',3I5)
  434.   894 IF(DZ(I,J,K).GE.0.0) GO TO 896
  435.       IFATAL=IFATAL+1
  436.       WRITE(IOCODE,895) I,J,K
  437.   895 FORMAT(/,5X,5('-'),'GRID BLOCK DZ ERROR AT IJK =',3I5)
  438.   896 IF(DZNET(I,J,K).GE.0.0) GO TO 900
  439.       IFATAL=IFATAL+1
  440.       WRITE(IOCODE,897) I,J,K
  441.   897 FORMAT(/,5X,5('-'),'GRID BLOCK DZNET ERROR AT IJK =',3I5)
  442.   900 CONTINUE
  443. C*****ESTABLISH NODE MID-POINT ELEVATION.
  444.       READ(20,69)
  445.       READ(20,*)KEL
  446.       IF(KEL.GT.1)GO TO 951
  447.       IF(KEL.EQ.1)GO TO 920
  448.       READ(20,*)ELEV
  449.       DO 910 J=1,JJ
  450.       DO 910 I=1,II
  451.       VAREL(I,J)=ELEV
  452. 910   CONTINUE
  453.   920 IF(KEL.NE.1)GO TO 930
  454.       DO 922 J=1,JJ
  455.       READ(20,*)(VAREL(I,J),I=1,II)
  456.   922 CONTINUE
  457.   930 CONTINUE
  458.       DO 923 I=1,II
  459.       DO 923 J=1,JJ
  460.   923 SUM(I,J)=0.
  461.       DO 926 K=1,KK
  462.       DO 926 J=1,JJ
  463.       DO 926 I=1,II
  464.       IF(K.GE.2) GO TO 924
  465.       EL(I,J,K)=VAREL(I,J)
  466.       GO TO 926
  467.   924 CONTINUE
  468.       DEL=SUM(I,J)+DZ(I,J,K-1)
  469.       EL(I,J,K)=VAREL(I,J)+DEL
  470.       SUM(I,J)=DZ(I,J,K-1)+SUM(I,J)
  471.   926 CONTINUE
  472.       GO TO 601
  473.   951 IF(KEL.EQ.3)GO TO 961
  474.       READ(20,*)(RKEL(K),K=1,KK)
  475.       DO 953 K=1,KK
  476.       DO 953 J=1,JJ
  477.       DO 953 I=1,II
  478.   953 EL(I,J,K)=RKEL(K)
  479.       GO TO 971
  480.   961 DO 963 K=1,KK
  481.       DO 965 J=1,JJ
  482.   965 READ(20,*)(EL(I,J,K),I=1,II)
  483.   963 CONTINUE
  484.   971 CONTINUE
  485.   601 WRITE(IOCODE,390)
  486.       DO 600 K=1,KK
  487.        WRITE(IOCODE,38)K
  488.       DO 600 J=1,JJ
  489.       WRITE(IOCODE,72)(EL(I,J,K),I=1,II)
  490. 600   CONTINUE
  491.       DO 973 K=1,KK
  492.       DO 973 J=1,JJ
  493.       DO 973 I=1,II
  494.   973 EL(I,J,K)=EL(I,J,K)+DZ(I,J,K)*0.5
  495. 69    FORMAT(40A2)
  496. 56    FORMAT(//)
  497.    72 FORMAT(1X,15F8.0)
  498. 29    FORMAT(T15,'GRID BLOCK LENGTH (DX) IS INITIALLY',
  499.      &' SET AT',F10.4,' FOR ALL NODES'//)
  500.    31 FORMAT(//T15,'********GRID BLOCK LENGTH (DX) NODE MODIFICATIONS',
  501.      & '**********',//T15,
  502.      & '   I1  I2  J1  J2  K1  K2  NEW DX VALUE')
  503.    32 FORMAT(15X,6I4,2X,E10.4)
  504. 33    FORMAT(T15,'GRID BLOCK WIDTH  (DY) IS INITIALLY',
  505.      &' SET AT',F10.4,' FOR ALL NODES'//)
  506.    34 FORMAT(//T15,'********GRID BLOCK WIDTH  (DY) NODE MODIFICATIONS',
  507.      & '**********',//T15,
  508.      & '   I1  I2  J1  J2  K1  K2  NEW DY VALUE')
  509. 36    FORMAT(T15,'GRID BLOCK GROSS THICKNESS (DZ) IS INITIALLY',
  510.      &' SET AT',F10.4,' FOR ALL NODES'//)
  511.    37 FORMAT(//T15,'********GRID BLOCK GROSS THICKNESS (DZ) NODE ',
  512.      & 'MODIFICATIONS**********',//T15,
  513.      & '   I1  I2  J1  J2  K1  K2  NEW DZ VALUE')
  514. 38    FORMAT(/1X,'K =',I2/)
  515. 390   FORMAT(//T15,'********** DEPTHS TO GRID BLOCK TOPS **********'/)
  516. 43    FORMAT(//T15,'********GRID BLOCK LENGTH (DX) DISTRIBUTION********'
  517.      &/)
  518. 47    FORMAT(//T15,'********GRID BLOCK WIDTH  (DY) DISTRIBUTION********'
  519.      &/)
  520.    48 FORMAT(//T15,8('*'),'GRID BLOCK GROSS THICKNESS (DZ) ',
  521.      & 'DISTRIBUTION',8('*'),/)
  522.   511 FORMAT(T15,'GRID SIZE (DX) IN COLUMN',I5,' IS INITIALLY SET AT'
  523.      &,F8.2,' FOR ALL NODES',/)
  524.   512 FORMAT(T15,'GRID SIZE (DY) IN ROW   ',I5,' IS INITIALLY SET AT'
  525.      &,F8.2,' FOR ALL NODES',/)
  526.   513 FORMAT(T15,'GRID SIZE (DZ) IN LAYER ',I5,' IS INITIALLY SET AT'
  527.      &,F8.2,' FOR ALL NODES',/)
  528.       RETURN
  529.       END
  530. C.................................................................INTCOM
  531.       SUBROUTINE INTCOM(IREG,X,Y,N,XO,YO)
  532. C      MACHINE DEPENDENT INCLUDE STATEMENT
  533. $INCLUDE:'PARAMS.FOR'
  534.       DIMENSION X(LP10,LP9),Y(LP10,LP9)
  535.       IF(XO.GE.X(IREG,N)) YO=Y(IREG,N)
  536.       IF(XO.GE.X(IREG,N)) RETURN
  537.       DO 10 I=2,N
  538.       IF(XO.GT.X(IREG,I)) GO TO 10
  539.       YO=Y(IREG,I)
  540.       RETURN
  541.   10  CONTINUE
  542.       END
  543. C.................................................................INTERP
  544.       SUBROUTINE INTERP(IREG,X,Y,N,XO,YO)
  545. C      MACHINE DEPENDENT INCLUDE STATEMENT
  546. $INCLUDE:'PARAMS.FOR'
  547.       DIMENSION X(LP10,LP9),Y(LP10,LP9)
  548.       IF(XO.GE.X(IREG,N)) YO=Y(IREG,N)
  549.       IF(XO.GE.X(IREG,N)) RETURN
  550.       DO 10 I=2,N
  551.       IF(XO.GE.X(IREG,I)) GO TO 10
  552.       YO=Y(IREG,I-1) +(XO-X(IREG,I-1))*(Y(IREG,I)-Y(IREG,I-1))
  553.      & /(X(IREG,I)-X(IREG,I-1))
  554.       RETURN
  555.   10  CONTINUE
  556.       END
  557. C.................................................................INTPVT
  558.       SUBROUTINE INTPVT(IREG,BPT,RM,X,Y,N,XO,YO)
  559. C      MACHINE DEPENDENT INCLUDE STATEMENT
  560. $INCLUDE:'PARAMS.FOR'
  561.       DIMENSION X(LP10,LP9),Y(LP10,LP9)
  562.       IF(XO.GT.BPT)GO TO 100
  563.       IF(XO.GE.X(IREG,N))YO=Y(IREG,N)
  564.       IF(XO.GE.X(IREG,N))RETURN
  565.       DO 10 I=2,N
  566.       IF(XO.GE.X(IREG,I))GO TO 10
  567.       YO=Y(IREG,I-1)+(XO-X(IREG,I-1))*(Y(IREG,I)-Y(IREG,I-1))
  568.      & /(X(IREG,I)-X(IREG,I-1))
  569.       RETURN
  570.    10 CONTINUE
  571.   100 CONTINUE
  572.       DO 200 I=2,N
  573.       IF(BPT.GE.X(IREG,I))GO TO 200
  574.       YOBP=Y(IREG,I-1)+(BPT-X(IREG,I-1))*(Y(IREG,I)-Y(IREG,I-1))
  575.      & /(X(IREG,I)-X(IREG,I-1))
  576.       YO=YOBP+RM*(XO-BPT)
  577.       RETURN
  578.   200 CONTINUE
  579.       END
  580. C........................................................................LSORX
  581.       SUBROUTINE LSORX(NX,NY,NZ,OMEGA,TOL,TOL1,MITER,DELT,DELT0,KSN,
  582.      & N,NITER)
  583. C      MACHINE DEPENDENT INCLUDE STATEMENT
  584. $INCLUDE:'PARAMS.FOR'
  585. C      ITERATIVE SOLUTION: X-DIRECTION TRIDIAGONAL ALGORITHM
  586.       COMMON /COEF/ AW(LP1,LP2,LP3),AE(LP1,LP2,LP3),AN(LP1,LP2,LP3),
  587.      & AS(LP1,LP2,LP3),AB(LP1,LP2,LP3),AT(LP1,LP2,LP3),E(LP1,LP2,LP3),
  588.      & B(LP1,LP2,LP3)
  589.       COMMON /SARRAY/ PN(LP1,LP2,LP3),IOCODE,IDMAX,
  590.      & SON(LP1,LP2,LP3),SWN(LP1,LP2,LP3),SGN(LP1,LP2,LP3),
  591.      & A1(LP1,LP2,LP3),A2(LP1,LP2,LP3),A3(LP1,LP2,LP3),
  592.      & SUM(LP1,LP2,LP3),GAM(LP1,LP2,LP3),QS(LP1,LP2,LP3)
  593.       COMMON /SPRTPS/ P(LP1,LP2,LP3),SO(LP1,LP2,LP3),SW(LP1,LP2,LP3),
  594.      & SG(LP1,LP2,LP3)
  595.       COMMON /VECTOR/ DX(LP1,LP2,LP3),DY(LP1,LP2,LP3),DZ(LP1,LP2,LP3),
  596.      & DZNET(LP1,LP2,LP3),IQN1(LP11),IQN2(LP11),IQN3(LP11),IHEDIN(80)
  597.       COMMON /TRIDI/ UM(LP15),AZL(LP15),BZL(LP15),CZL(LP15),DZL(LP15),
  598.      & UZL(LP15),X(LP15),BETA(LP15),GAMMA(LP15),W(LP15)
  599.       DIV=DELT/DELT0
  600.       NITER=0
  601.       DMAX=1.0
  602.       RHO1=0.0
  603.       THETA=0.0
  604.    11 CONTINUE
  605.       TW = 1.0 - OMEGA
  606.       DMAX0=DMAX
  607.       THETA0=THETA
  608.       IF(NITER.GE.MITER)WRITE(IOCODE,30)NITER,TOL,DMAX
  609.       IF(NITER.GE.MITER)RETURN
  610.       NITER = NITER +1
  611.       DMAX=0.0
  612.       DO 20 K=1,NZ
  613.       DO 20 J=1,NY
  614.       DO 15 I=1,NX
  615.       UM(I) = P(I,J,K)
  616.       AZL(I) = AW(I,J,K)
  617.       BZL(I)=E(I,J,K)
  618.       CZL(I)=AE(I,J,K)
  619.       DZL(I)=B(I,J,K)
  620.       IF(NY.EQ.1) GO TO 14
  621.       JM=J-1
  622.       JP=J+1
  623.       IF(J.EQ.1) JM=1
  624.       IF(J.EQ.NY) JP=NY
  625.       DZL(I)=DZL(I) - AS(I,J,K)*P(I,JM,K) - AN(I,J,K)*P(I,JP,K)
  626.    14 IF(NZ.EQ.1) GO TO 15
  627.       KM=K-1
  628.       KP=K+1
  629.       IF(K.EQ.1)KM=1
  630.       IF(K.EQ.NZ) KP=NZ
  631.       DZL(I)=DZL(I) - AT(I,J,K)*P(I,J,KM) - AB(I,J,K)*P(I,J,KP)
  632. 15    CONTINUE
  633. C     THIS ALGORITHM SOLVES THE TRIDIAGONAL SYSTEM
  634. C     GENERATED BY THE SYSTEM OF N EQUATIONS
  635. C     AZL(I)*U(I-1) + BZL(I)*U(I) + CZL(I)*U(I+1) = DZL(I)
  636.       BETA(1)=BZL(1)
  637.       GAMMA(1)=DZL(1)/BZL(1)
  638.       NXM=NX-1
  639. C      COMPUTE FORWARD SOLUTION.
  640.       DO 1010 ITRI=1,NXM
  641.       W(ITRI)=CZL(ITRI)/BETA(ITRI)
  642.       ITRIP=ITRI+1
  643.  1010 BETA(ITRIP)=BZL(ITRIP)-AZL(ITRIP)*W(ITRI)
  644.       DO 1020 ITRI=2,NX
  645.       ITRIM=ITRI-1
  646.  1020 GAMMA(ITRI)=(DZL(ITRI)-AZL(ITRI)*GAMMA(ITRIM))/BETA(ITRI)
  647. C      COMPUTE BACK SOLUTION
  648.       UZL(NX)=GAMMA(NX)
  649.       DO 1030 JTRI=1,NXM
  650.       ITRI=NX-JTRI
  651.       ITRIP=ITRI+1
  652.  1030 UZL(ITRI)=GAMMA(ITRI)-W(ITRI)*UZL(ITRIP)
  653.       DO 16 I=1,NX
  654.       GSLSOR=UZL(I)
  655.       P(I,J,K) = TW*UM(I) + OMEGA*GSLSOR
  656.       ARG=P(I,J,K) - UM(I)
  657.       DM=ABS(ARG)
  658.       IF(DM.GT.DMAX) DMAX = DM
  659. 16    CONTINUE
  660. 20    CONTINUE
  661.       IF(TOL1.EQ.0.0)GO TO 25
  662.       THETA=DMAX/DMAX0
  663.       DELTA=THETA-THETA0
  664.       ARG=DELTA
  665.       ARG=ABS(ARG)
  666.       IF(ARG.GT.TOL1)GO TO 25
  667.       OM=OMEGA-1.0
  668.       RHO1=(THETA+OM)*(THETA+OM)/(THETA*OMEGA*OMEGA)
  669.       IF(RHO1.GE.1.0)GO TO 25
  670.       ARG=1.0-RHO1
  671.       OMEGA=2.0/(1.0+SQRT(ARG))
  672. 25    IF(DMAX.GT.TOL)GO TO 11
  673.       IF(N.NE.KSN) GO TO 300
  674.       WRITE(IOCODE,40)NITER,OMEGA,DMAX,THETA,RHO1
  675.  300  CONTINUE
  676. 40    FORMAT(T5,'CONVERGENCE(LSORX) HAS BEEN REACHED AFTER ',I3,
  677.      &' ITERATIONS',5X,'OMEGA = ',F6.3/
  678.      &T5,'DMAX = ',F10.6,5X,'THETA = ',F10.6,5X,'RHO1 = ',F10.6/)
  679. 30    FORMAT(T15,'CONVERGENCE(LSORX) WAS NOT REACHED IN ',I5,
  680.      &' ITERATIONS'/T15,'TOL = ',F10.7,10X,'DMAX = ',F15.7)
  681.       RETURN
  682.       END
  683. C........................................................................LSORY
  684.       SUBROUTINE LSORY(NX,NY,NZ,OMEGA,TOL,TOL1,MITER,DELT,DELT0,KSN,
  685.      & N,NITER)
  686. C      MACHINE DEPENDENT INCLUDE STATEMENT
  687. $INCLUDE:'PARAMS.FOR'
  688. C      ITERATIVE SOLUTION: Y-DIRECTION TRIDIAGONAL ALGORITHM
  689.       COMMON /COEF/ AW(LP1,LP2,LP3),AE(LP1,LP2,LP3),AN(LP1,LP2,LP3),
  690.      & AS(LP1,LP2,LP3),AB(LP1,LP2,LP3),AT(LP1,LP2,LP3),E(LP1,LP2,LP3),
  691.      & B(LP1,LP2,LP3)
  692.       COMMON /SARRAY/ PN(LP1,LP2,LP3),IOCODE,IDMAX,
  693.      & SON(LP1,LP2,LP3),SWN(LP1,LP2,LP3),SGN(LP1,LP2,LP3),
  694.      & A1(LP1,LP2,LP3),A2(LP1,LP2,LP3),A3(LP1,LP2,LP3),
  695.      & SUM(LP1,LP2,LP3),GAM(LP1,LP2,LP3),QS(LP1,LP2,LP3)
  696.       COMMON /SPRTPS/ P(LP1,LP2,LP3),SO(LP1,LP2,LP3),SW(LP1,LP2,LP3),
  697.      & SG(LP1,LP2,LP3)
  698.       COMMON /VECTOR/ DX(LP1,LP2,LP3),DY(LP1,LP2,LP3),DZ(LP1,LP2,LP3),
  699.      & DZNET(LP1,LP2,LP3),IQN1(LP11),IQN2(LP11),IQN3(LP11),IHEDIN(80)
  700.       COMMON /TRIDI/ UM(LP15),AZL(LP15),BZL(LP15),CZL(LP15),DZL(LP15),
  701.      & UZL(LP15),X(LP15),BETA(LP15),GAMMA(LP15),W(LP15)
  702.       DIV=DELT/DELT0
  703.       NITER=0
  704.       DMAX=1.0
  705.       RHO1=0.0
  706.       THETA=0.0
  707.    11 CONTINUE
  708.       TW = 1.0 - OMEGA
  709.       DMAX0=DMAX
  710.       THETA0=THETA
  711.       IF(NITER.GE.MITER)WRITE(IOCODE,30)NITER,TOL,DMAX
  712.       IF(NITER.GE.MITER)RETURN
  713.       NITER = NITER +1
  714.       DMAX=0.0
  715.       DO 20 K=1,NZ
  716.       DO 20 I=1,NX
  717.       DO 15 J=1,NY
  718.       UM(J) = P(I,J,K)
  719.       AZL(J) = AS(I,J,K)
  720.       BZL(J)=E(I,J,K)
  721.       CZL(J)=AN(I,J,K)
  722.       DZL(J)=B(I,J,K)
  723.       IF(NX.EQ.1) GO TO 14
  724.       IM=I-1
  725.       IP=I+1
  726.       IF(I.EQ.1) IM=1
  727.       IF(I.EQ.NX) IP=NX
  728.       DZL(J)=DZL(J) - AW(I,J,K)*P(IM,J,K) - AE(I,J,K)*P(IP,J,K)
  729.    14 IF(NZ.EQ.1) GO TO 15
  730.       KM=K-1
  731.       KP=K+1
  732.       IF(K.EQ.1)KM=1
  733.       IF(K.EQ.NZ) KP=NZ
  734.       DZL(J)=DZL(J) - AT(I,J,K)*P(I,J,KM) - AB(I,J,K)*P(I,J,KP)
  735. 15    CONTINUE
  736. C     THIS ALGORITHM SOLVES THE TRIDIAGONAL SYSTEM
  737. C     GENERATED BY THE SYSTEM OF N EQUATIONS
  738. C     AZL(J)*U(J-1) + BZL(J)*U(J) + CZL(J)*U(J+1) = DZL(J)
  739.       BETA(1)=BZL(1)
  740.       GAMMA(1)=DZL(1)/BZL(1)
  741.       NYM=NY-1
  742. C      COMPUTE FORWARD SOLUTION.
  743.       DO 1010 ITRI=1,NYM
  744.       W(ITRI)=CZL(ITRI)/BETA(ITRI)
  745.       ITRIP=ITRI+1
  746.  1010 BETA(ITRIP)=BZL(ITRIP)-AZL(ITRIP)*W(ITRI)
  747.       DO 1020 ITRI=2,NY
  748.       ITRIM=ITRI-1
  749.  1020 GAMMA(ITRI)=(DZL(ITRI)-AZL(ITRI)*GAMMA(ITRIM))/BETA(ITRI)
  750. C      COMPUTE BACK SOLUTION
  751.       UZL(NY)=GAMMA(NY)
  752.       DO 1030 JTRI=1,NYM
  753.       ITRI=NY-JTRI
  754.       ITRIP=ITRI+1
  755.  1030 UZL(ITRI)=GAMMA(ITRI)-W(ITRI)*UZL(ITRIP)
  756.       DO 16 J=1,NY
  757.       GSLSOR=UZL(J)
  758.       P(I,J,K) = TW*UM(J) + OMEGA*GSLSOR
  759.       ARG=P(I,J,K) - UM(J)
  760.       DM=ABS(ARG)
  761.       IF(DM.GT.DMAX) DMAX = DM
  762. 16    CONTINUE
  763. 20    CONTINUE
  764.       IF(TOL1.EQ.0.0)GO TO 25
  765.       THETA=DMAX/DMAX0
  766.       DELTA=THETA-THETA0
  767.       ARG=DELTA
  768.       ARG=ABS(ARG)
  769.       IF(ARG.GT.TOL1)GO TO 25
  770.       OM=OMEGA-1.0
  771.       RHO1=(THETA+OM)*(THETA+OM)/(THETA*OMEGA*OMEGA)
  772.       IF(RHO1.GE.1.0)GO TO 25
  773.       ARG=1.0-RHO1
  774.       OMEGA=2.0/(1.0+SQRT(ARG))
  775. 25    IF(DMAX.GT.TOL)GO TO 11
  776.       IF(N.NE.KSN) GO TO 300
  777.       WRITE(IOCODE,40)NITER,OMEGA,DMAX,THETA,RHO1
  778.  300  CONTINUE
  779. 40    FORMAT(T5,'CONVERGENCE(LSORY) HAS BEEN REACHED AFTER ',I3,
  780.      &' ITERATIONS',5X,'OMEGA = ',F6.3/
  781.      &T5,'DMAX = ',F10.6,5X,'THETA = ',F10.6,5X,'RHO1 = ',F10.6/)
  782. 30    FORMAT(T15,'CONVERGENCE(LSORY) WAS NOT REACHED IN ',I5,
  783.      &' ITERATIONS'/T15,'TOL = ',F10.7,10X,'DMAX = ',F15.7)
  784.       RETURN
  785.       END
  786. C.................................................................LSORZ
  787.       SUBROUTINE LSORZ(NX,NY,NZ,OMEGA,TOL,TOL1,MITER,DELT,DELT0,KSN,
  788.      & N,NITER)
  789. C      MACHINE DEPENDENT INCLUDE STATEMENT
  790. $INCLUDE:'PARAMS.FOR'
  791. C      ITERATIVE SOLUTION: Z-DIRECTION TRIDIAGONAL ALGORITHM
  792.       COMMON /COEF/ AW(LP1,LP2,LP3),AE(LP1,LP2,LP3),AN(LP1,LP2,LP3),
  793.      & AS(LP1,LP2,LP3),AB(LP1,LP2,LP3),AT(LP1,LP2,LP3),E(LP1,LP2,LP3),
  794.      & B(LP1,LP2,LP3)
  795.       COMMON /SARRAY/ PN(LP1,LP2,LP3),IOCODE,IDMAX,
  796.      & SON(LP1,LP2,LP3),SWN(LP1,LP2,LP3),SGN(LP1,LP2,LP3),
  797.      & A1(LP1,LP2,LP3),A2(LP1,LP2,LP3),A3(LP1,LP2,LP3),
  798.      & SUM(LP1,LP2,LP3),GAM(LP1,LP2,LP3),QS(LP1,LP2,LP3)
  799.       COMMON /SPRTPS/ P(LP1,LP2,LP3),SO(LP1,LP2,LP3),SW(LP1,LP2,LP3),
  800.      & SG(LP1,LP2,LP3)
  801.       COMMON /VECTOR/ DX(LP1,LP2,LP3),DY(LP1,LP2,LP3),DZ(LP1,LP2,LP3),
  802.      & DZNET(LP1,LP2,LP3),IQN1(LP11),IQN2(LP11),IQN3(LP11),IHEDIN(80)
  803.       COMMON /TRIDI/ UM(LP15),AZL(LP15),BZL(LP15),CZL(LP15),DZL(LP15),
  804.      & UZL(LP15),X(LP15),BETA(LP15),GAMMA(LP15),W(LP15)
  805.       DIV=DELT/DELT0
  806.       NITER=0
  807.       DMAX=1.0
  808.       RHO1=0.0
  809.       THETA=0.0
  810.    11 CONTINUE
  811.       TW = 1.0 - OMEGA
  812.       DMAX0=DMAX
  813.       THETA0=THETA
  814.       IF(NITER.GE.MITER)WRITE(IOCODE,30)NITER,TOL,DMAX
  815.       IF(NITER.GE.MITER)RETURN
  816.       NITER = NITER +1
  817.       DMAX=0.0
  818.       DO 20 I=1,NX
  819.       DO 20 J=1,NY
  820.       DO 15 K=1,NZ
  821.       UM(K) = P(I,J,K)
  822.       AZL(K) = AT(I,J,K)
  823.       BZL(K)=E(I,J,K)
  824.       CZL(K)=AB(I,J,K)
  825.       DZL(K)=B(I,J,K)
  826.       IF(NY.EQ.1) GO TO 14
  827.       JM=J-1
  828.       JP=J+1
  829.       IF(J.EQ.1) JM=1
  830.       IF(J.EQ.NY) JP=NY
  831.       DZL(K)=DZL(K) - AS(I,J,K)*P(I,JM,K) - AN(I,J,K)*P(I,JP,K)
  832.    14 IF(NX.EQ.1) GO TO 15
  833.       IM=I-1
  834.       IP=I+1
  835.       IF(I.EQ.1)IM=1
  836.       IF(I.EQ.NX) IP=NX
  837.       DZL(K)=DZL(K) - AW(I,J,K)*P(IM,J,K) - AE(I,J,K)*P(IP,J,K)
  838.    15 CONTINUE
  839. C     THIS ALGORITHM SOLVES THE TRIDIAGONAL SYSTEM
  840. C     GENERATED BY THE SYSTEM OF N EQUATIONS
  841. C     AZL(K)*U(K-1) + BZL(K)*U(K) + CZL(K)*U(K+1) = DZL(K)
  842.       BETA(1)=BZL(1)
  843.       GAMMA(1)=DZL(1)/BZL(1)
  844.       NZM=NZ-1
  845. C      COMPUTE FORWARD SOLUTION.
  846.       DO 1010 ITRI=1,NZM
  847.       W(ITRI)=CZL(ITRI)/BETA(ITRI)
  848.       ITRIP=ITRI+1
  849.  1010 BETA(ITRIP)=BZL(ITRIP)-AZL(ITRIP)*W(ITRI)
  850.       DO 1020 ITRI=2,NZ
  851.       ITRIM=ITRI-1
  852.  1020 GAMMA(ITRI)=(DZL(ITRI)-AZL(ITRI)*GAMMA(ITRIM))/BETA(ITRI)
  853. C      COMPUTE BACK SOLUTION
  854.       UZL(NZ)=GAMMA(NZ)
  855.       DO 1030 JTRI=1,NZM
  856.       ITRI=NZ-JTRI
  857.       ITRIP=ITRI+1
  858.  1030 UZL(ITRI)=GAMMA(ITRI)-W(ITRI)*UZL(ITRIP)
  859.       DO 16 K=1,NZ
  860.       GSLSOR=UZL(K)
  861.       P(I,J,K) = TW*UM(K) + OMEGA*GSLSOR
  862.       ARG=P(I,J,K) - UM(K)
  863.       DM=ABS(ARG)
  864.       IF(DM.GT.DMAX) DMAX = DM
  865. 16    CONTINUE
  866. 20    CONTINUE
  867.       IF(TOL1.EQ.0.0)GO TO 25
  868.       THETA=DMAX/DMAX0
  869.       DELTA=THETA-THETA0
  870.       ARG=DELTA
  871.       ARG=ABS(ARG)
  872.       IF(ARG.GT.TOL1)GO TO 25
  873.       OM=OMEGA-1.0
  874.       RHO1=(THETA+OM)*(THETA+OM)/(THETA*OMEGA*OMEGA)
  875.       IF(RHO1.GE.1.0)GO TO 25
  876.       ARG=1.0-RHO1
  877.       OMEGA=2.0/(1.0+SQRT(ARG))
  878. 25    IF(DMAX.GT.TOL)GO TO 11
  879.       IF(N.NE.KSN) GO TO 300
  880.       WRITE(IOCODE,40)NITER,OMEGA,DMAX,THETA,RHO1
  881.  300  CONTINUE
  882. 40    FORMAT(T5,'CONVERGENCE(LSORZ) HAS BEEN REACHED AFTER ',I3,
  883.      &' ITERATIONS',5X,'OMEGA = ',F6.3/
  884.      &T5,'DMAX = ',F10.6,5X,'THETA = ',F10.6,5X,'RHO1 = ',F10.6/)
  885. 30    FORMAT(T15,'CONVERGENCE(LSORZ) WAS NOT REACHED IN ',I5,
  886.      &' ITERATIONS'/T15,'TOL = ',F10.7,10X,'DMAX = ',F15.7)
  887.       RETURN
  888.       END
  889.