$DO66
C.................................................................CODES
      SUBROUTINE CODES(IOCODE,KSM1,KSN1,KCO1,NN,
     &FACT1,FACT2,TMAX,KSOL,MITER,OMEGA,TOL,TOL1,
     & KSN,KSM,KCO,KCOFF,DSMAX,DPMAX,
     & NUMDIS,IRK,THRUIN,
     &WORMAX,GORMAX,PAMIN,PAMAX)
C      MACHINE DEPENDENT INCLUDE STATEMENT
$INCLUDE:'PARAMS.FOR'
      REAL KROT,KRWT,KRGT,KROGT,MUOT,MUWT,MUGT
      COMMON /SPVT/ SAT(LP7,LP9),KROT(LP7,LP9),KRWT(LP7,LP9),
     & BGT(LP7,LP9),
     & KRGT(LP7,LP9),ITHREE(LP7),RSOT(LP7,LP9),BWPT(LP7,LP9),
     & PCOWT(LP7,LP9),PCGOT(LP7,LP9),KROGT(LP7,LP9),SWR(LP7),
     & POT(LP7,LP9),MUOT(LP7,LP9),BOT(LP7,LP9),BOPT(LP7,LP9),
     & RSOPT(LP7,LP9),PWT(LP7,LP9),MUWT(LP7,LP9),BWT(LP7,LP9),
     & RSWT(LP7,LP9),RSWPT(LP7,LP9),PGT(LP7,LP9),MUGT(LP7,LP9),
     & BGPT(LP7,LP9),CRT(LP7,LP9),IPVT(LP1,LP2,LP3),IROCK(LP1,LP2,LP3),
     & NROCK,NPVT,PSIT(LP7,LP9),PRT(LP7,LP9),WOROCK(LP7),GOROCK(LP7)
C      READ DEBUG, MATRIX SOLUTION AND RUN CONTROL INFO
      COMMON /TSTDAT/ IFATAL,IWARN
      READ(20,69)
      READ(20,*) KSN1,KSM1,KCO1,KCOFF
C EVERY KSM1 TH STEP SOLUTION MATRIX WILL BE WRITTEN
C EVERY KSN1 TH STEP LSORX DATA WILL BE WRITTEN
C EVERY KCO1 TH STEP COMPRESSIBILITIES & VOLUME FACTORS WILL BE WRITTEN
      READ(20,69)
      READ(20,*)NN,FACT1,FACT2,TMAX,WORMAX,GORMAX,PAMIN,PAMAX
      WRITE(IOCODE,59)NN,FACT1,FACT2,TMAX,WORMAX,GORMAX,PAMIN,PAMAX
      IF(WORMAX.NE.0.0) GO TO 20
      READ(20,*) (WOROCK(I),I=1,NROCK)
      WRITE(IOCODE,5)
    5 FORMAT(//T15,'ROCK REGION SPECIFIED WOR MAXIMA:')
      DO 10 I=1,NROCK
      WRITE(IOCODE,8) I,WOROCK(I)
    8 FORMAT(T20,'ROCK REGION ',I3,' WOR MAX (STB/STB) = ',F8.0)
   10 CONTINUE
   20 IF(GORMAX.NE.0.0) GO TO 40
      READ(20,*) (GOROCK(I),I=1,NROCK)
      WRITE(IOCODE,25)
   25 FORMAT(//T15,'ROCK REGION SPECIFIED GOR MAXIMA:')
      DO 30 I=1,NROCK
      WRITE(IOCODE,28) I,GOROCK(I)
   28 FORMAT(T20,'ROCK REGION ',I3,' GOR MAX (SCF/STB) = ',F8.0)
   30 CONTINUE
   40 CONTINUE
C NN--MAXIMUM NUMBER OF TIME-STEPS
C FACT1--FACTOR FOR INCREASING TIME-STEP SIZE
C FACT2--FACTOR FOR DECREASING TIME-STEP SIZE
C TMAX--MAXIMUM SIMULATION TIME
      READ(20,69)
      READ(20,*)KSOL,MITER,OMEGA,TOL,TOL1,DSMAX,DPMAX
      IF(KSOL.EQ.1) WRITE(IOCODE,73)
      IF(KSOL.EQ.2) WRITE(IOCODE,75) MITER,OMEGA,TOL,TOL1
      IF(KSOL.EQ.3) WRITE(IOCODE,76) MITER,OMEGA,TOL,TOL1
      IF(KSOL.EQ.4) WRITE(IOCODE,77) MITER,OMEGA,TOL,TOL1
      IF(KSOL.GT.4) THEN
        WRITE(IOCODE,*) 'IMPROPER SOLUTION METHOD SPECIFIED,KSOL=',KSOL
        IFATAL=IFATAL + 1
      ENDIF
      WRITE(IOCODE,99) DSMAX,DPMAX
   73 FORMAT(//T15,'SOLUTION METHOD IS BAND.')
   75 FORMAT(//T15,'SOLUTION METHOD IS LSORX:',
     &/T20,'MAXIMUM NUMBER OF ITERATIONS      (MITER) = ',5X,I5,
     &/T20,'INITIAL ACCELERATION PARAMETER    (OMEGA) = ',F10.4,
     &/T20,'MAXIMUM PRESSURE RESIDUAL           (TOL) = ',F10.4,
     &/T20,'PARAMETER FOR CHANGING OMEGA       (TOL1) = ',F10.4)
   76 FORMAT(//T15,'SOLUTION METHOD IS LSORY:',
     &/T20,'MAXIMUM NUMBER OF ITERATIONS      (MITER) = ',5X,I5,
     &/T20,'INITIAL ACCELERATION PARAMETER    (OMEGA) = ',F10.4,
     &/T20,'MAXIMUM PRESSURE RESIDUAL           (TOL) = ',F10.4,
     &/T20,'PARAMETER FOR CHANGING OMEGA       (TOL1) = ',F10.4)
   77 FORMAT(//T15,'SOLUTION METHOD IS LSORZ:',
     &/T20,'MAXIMUM NUMBER OF ITERATIONS      (MITER) = ',5X,I5,
     &/T20,'INITIAL ACCELERATION PARAMETER    (OMEGA) = ',F10.4,
     &/T20,'MAXIMUM PRESSURE RESIDUAL           (TOL) = ',F10.4,
     &/T20,'PARAMETER FOR CHANGING OMEGA       (TOL1) = ',F10.4)
   99 FORMAT(/T15,'AUTOMATIC TIME STEP CRITERIA:',
     &/T20,'MAXIMUM ALLOWED SATURATION CHANGE (DSMAX) = ',F10.4,
     &/T20,'MAXIMUM ALLOWED PRESSURE CHANGE   (DPMAX) = ',F10.4/)
   59 FORMAT(1H1/T15,'RUN CONTROL PARAMETERS:',/,
     & T20,'MAXIMUM NUMBER OF TIME-STEPS =',I5/
     &T20,'FACTOR FOR INCREASING DELT =',F10.3,3X,
     &'WHEN DSMAX AND DPMAX NOT EXCEEDED.',/,
     &T20,'FACTOR FOR DECREASING DELT =',F10.3,3X,
     &'WHEN DSMAX OR DPMAX IS EXCEEDED.',/,
     &T20,'MAXIMUM SIMULATION TIME =',F11.3/,
     &T20,'MAXIMUM RESERVOIR WOR/TIME-STEP =',F8.0,' STB/STB'/
     &T20,'MAXIMUM RESERVOIR GOR/TIME-STEP =',F8.0,' SCF/STB'/
     &T20,'MINIMUM AVERAGE RESERVOIR PRESSURE/TIME-STEP =',F8.0/
     &T20,'MAXIMUM AVERAGE RESERVOIR PRESSURE/TIME-STEP =',F8.0)
      KSN=KSN1
      KSM=KSM1
      KCO=KCO1
69    FORMAT(40A2)
C      NUM DIS AND IMPES/R-K SOLN CONTROLS
      READ(20,69)
      READ(20,*) NUMDIS,IRK,THRUIN
      IF(IRK.EQ.0.AND.NUMDIS.EQ.0) WRITE(IOCODE,200)
  200 FORMAT(//T15,'IMPES FORMULATION SELECTED; ',
     & ' SINGLE-POINT UPSTREAM WEIGHTING.')
      IF(IRK.EQ.0.AND.NUMDIS.EQ.1) WRITE(IOCODE,220)
  220 FORMAT(//T15,'IMPES FORMULATION SELECTED; ',
     & ' TWO-POINT UPSTREAM WEIGHTING.')
      IF(IRK.GT.0.AND.NUMDIS.EQ.0) WRITE(IOCODE,240)
     & THRUIN
  240 FORMAT(//T15,'STABILISED IMPES FORMULATION:',
     & /T20,'SINGLE-POINT UPSTREAM WEIGHTING',
     & /T20,'USER SPEC THROUGHPUT, FRACTION  (THRUIN) =',5X,F10.4,/)
      IF(IRK.GT.0.AND.NUMDIS.EQ.1) WRITE(IOCODE,260)
     & THRUIN
  260 FORMAT(//T15,'STABILISED IMPES FORMULATION:',
     & /T20,'TWO-POINT UPSTREAM WEIGHTING',
     & /T20,'USER SPEC THROUGHPUT, FRACTION  (THRUIN) =',5X,F10.4,/)
C      CHECK CODE DATA
C      DEBUG
      IF(KSN1.EQ.0) GO TO 1020
      IWARN=IWARN+1
      WRITE(IOCODE,1010)
 1010 FORMAT(/5X,5('-'),'LSOR DEBUG OUTPUT ON')
 1020 IF(KSM1.EQ.0) GO TO 1040
      IWARN=IWARN+1
      WRITE(IOCODE,1030)
 1030 FORMAT(/5X,5('-'),'SOLN METHOD DEBUG OUTPUT ON')
 1040 IF(KCO1.EQ.0) GO TO 1060
      IWARN=IWARN+1
      WRITE(IOCODE,1050)
 1050 FORMAT(/5X,5('-'),'COMP AND FVF DEBUG OUTPUT ON')
 1060 IF(KCOFF.EQ.0) GO TO 1400
      IWARN=IWARN+1
      WRITE(IOCODE,1070)
 1070 FORMAT(/5X,5('-'),'DEN AND SAT DEBUG OUTPUT ON')
C      RUN CONTROL
 1400 CONTINUE
      IF(NN.GE.1) GO TO 1420
      IFATAL=IFATAL+1
      WRITE(IOCODE,1410)
 1410 FORMAT(/5X,5('-'),'MAX # OF TIME STEPS ERROR')
 1420 IF(FACT1.GE.1.0) GO TO 1440
      IFATAL=IFATAL+1
      WRITE(IOCODE,1430)
 1430 FORMAT(/5X,5('-'),'FACT1 ERROR')
 1440 IF(FACT2.GT.0.0.AND.FACT2.LE.1.0) GO TO 1460
      IFATAL=IFATAL+1
      WRITE(IOCODE,1450)
 1450 FORMAT(/5X,5('-'),'FACT2 ERROR')
 1460 IF(TMAX.GT.0.0) GO TO 1480
      IFATAL=IFATAL+1
      WRITE(IOCODE,1470)
 1470 FORMAT(/5X,5('-'),'TMAX ERROR')
 1480 IF(WORMAX.GE.0.0) GO TO 1500
      IFATAL=IFATAL+1
      WRITE(IOCODE,1490)
 1490 FORMAT(/5X,5('-'),'WORMAX ERROR')
 1500 IF(GORMAX.GE.0.0) GO TO 1520
      IFATAL=IFATAL+1
      WRITE(IOCODE,1510)
 1510 FORMAT(/5X,5('-'),'GORMAX ERROR')
 1520 IF(PAMIN.GT.0.0) GO TO 1540
      IFATAL=IFATAL+1
      WRITE(IOCODE,1530)
 1530 FORMAT(/5X,5('-'),'PAMIN ERROR')
 1540 IF(PAMAX.GT.PAMIN) GO TO 2000
      IFATAL=IFATAL+1
      WRITE(IOCODE,1550)
 1550 FORMAT(/5X,5('-'),'PAMAX ERROR')
C      SOLUTION METHOD
 2000 CONTINUE
      IF(MITER.GE.1) GO TO 2020
      IWARN=IWARN+1
      WRITE(IOCODE,2010)
 2010 FORMAT(/5X,5('-'),'ALLOWED # OF LSOR ITER LESS THAN 1')
 2020 IF(OMEGA.GE.1.0.AND.OMEGA.LE.2.0) GO TO 2040
      IWARN=IWARN+1
      WRITE(IOCODE,2030)
 2030 FORMAT(/5X,5('-'),'OMEGA LT 1 OR OMEGA GT 2')
 2040 IF(TOL.GT.0) GO TO 2048
      IWARN=IWARN+1
      WRITE(IOCODE,2046)
 2046 FORMAT(/5X,5('-'),'TOL LESS THAN OR EGUAL TO 0')
 2048 IF(TOL1.GE.0.0) GO TO 2060
      IWARN=IWARN+1
      WRITE(IOCODE,2050)
 2050 FORMAT(/5X,5('-'),'TOL1 IS NEGATIVE')
 2060 IF(DSMAX.GT.0.0.AND.DSMAX.LE.1.0) GO TO 2080
      IFATAL=IFATAL+1
      WRITE(IOCODE,2070)
 2070 FORMAT(/5X,5('-'),'DSMAX ERROR')
 2080 IF(DPMAX.GT.0.0) GO TO 3000
      IFATAL=IFATAL+1
      WRITE(IOCODE,2090)
 2090 FORMAT(/5X,5('-'),'DPMAX ERROR')
 3000 CONTINUE
      RETURN
      END
C................................................................GRIDSZ
      SUBROUTINE GRIDSZ(IOCODE,II,JJ,KK)
C      MACHINE DEPENDENT INCLUDE STATEMENT
$INCLUDE:'PARAMS.FOR'
C      READ GRID DESCRIPTION
      REAL KX,KY,KZ
      COMMON /TSTDAT/ IFATAL,IWARN
      COMMON /SPARM/ KX(LP1,LP2,LP3),KY(LP1,LP2,LP3),KZ(LP1,LP2,LP3),
     & EL(LP1,LP2,LP3),TX(LP4,LP2,LP3),TY(LP1,LP5,LP3),TZ(LP1,LP2,LP6),
     & PDAT(LP1,LP2,LP3),PDATUM,GRAD
      COMMON /VECTOR/ DX(LP1,LP2,LP3),DY(LP1,LP2,LP3),DZ(LP1,LP2,LP3),
     & DZNET(LP1,LP2,LP3),IQN1(LP11),IQN2(LP11),IQN3(LP11),IHEDIN(80)
      DIMENSION SUM(LP1,LP2),VAREL(LP1,LP2),RDXL(LP1),RDYL(LP2),
     & RDZL(LP3),RKEL(LP3),RDZNET(LP3)
      READ(20,69)
      WRITE(IOCODE,70)
   70 FORMAT(//T45,'***** INITIALIZATION DATA *****',//)
      READ(20,*) II,JJ,KK
      READ(20,69)
C*****READ INPUT CODES FOR DX,DY,DZ
      READ(20,*)KDX,KDY,KDZ,KDZNET
C*****ESTABLISH GRID BLOCK LENGTH (DX) DISTRIBUTION
      IF(KDX.GE.0)GO TO 180
      READ(20,*)DXC
      DO 175 K=1,KK
      DO 175 J=1,JJ
      DO 175 I=1,II
175   DX(I,J,K)=DXC
      WRITE(IOCODE,56)
      WRITE(IOCODE,29)DXC
      GO TO 195
  180 IF(KDX.GT.0)GO TO 185
      READ(20,*)(RDXL(I),I=1,II)
      DO 187 K=1,KK
      DO 187 J=1,JJ
      DO 187 I=1,II
  187 DX(I,J,K)=RDXL(I)
      DO 182 I=1,II
  182 WRITE(IOCODE,511)I,RDXL(I)
      GO TO 195
185   WRITE(IOCODE,43)
      K=1
      WRITE(IOCODE,38)K
      DO 190 J=1,JJ
      READ(20,*)(DX(I,J,K),I=1,II)
  190 WRITE(IOCODE,72)(DX(I,J,K),I=1,II)
      DO 194 K=2,KK
      WRITE(IOCODE,38) K
      DO 194 J=1,JJ
      DO 193 I=1,II
  193 DX(I,J,K)=DX(I,J,1)
  194 WRITE(IOCODE,72) (DX(I,J,K),I=1,II)
195   CONTINUE
      WRITE(IOCODE,56)
C*****ESTABLISH GRID BLOCK LENGTH (DY) DISTRIBUTION
      IF(KDY.GE.0)GO TO 200
      READ(20,*)DYC
      DO 202 K=1,KK
      DO 202 J=1,JJ
      DO 202 I=1,II
202   DY(I,J,K)=DYC
      WRITE(IOCODE,56)
      WRITE(IOCODE,33)DYC
      GO TO 220
  200 IF(KDY.GT.0)GO TO 207
      READ(20,*)(RDYL(J),J=1,JJ)
      DO 205 K=1,KK
      DO 205 J=1,JJ
      DO 205 I=1,II
  205 DY(I,J,K)=RDYL(J)
      DO 210 J=1,JJ
  210 WRITE(IOCODE,512)J,RDYL(J)
      GO TO 220
207   WRITE(IOCODE,47)
      K=1
      WRITE(IOCODE,38)K
      DO 215 J=1,JJ
      READ(20,*)(DY(I,J,K),I=1,II)
  215 WRITE(IOCODE,72)(DY(I,J,K),I=1,II)
      DO 214 K=2,KK
       WRITE(IOCODE,38) K
      DO 214 J=1,JJ
      DO 213 I=1,II
  213 DY(I,J,K)=DY(I,J,1)
  214 WRITE(IOCODE,72) (DY(I,J,K),I=1,II)
220   CONTINUE
      WRITE(IOCODE,56)
C*****ESTABLISH GRID BLOCK LENGTH (DZ) DISTRIBUTION
      IF(KDZ.GE.0)GO TO 225
      READ(20,*)DZC
      DO 230 K=1,KK
      DO 230 J=1,JJ
      DO 230 I=1,II
230   DZ(I,J,K)=DZC
      WRITE(IOCODE,56)
      WRITE(IOCODE,36)DZC
      GO TO 245
  225 IF(KDZ.GT.0)GO TO 232
      READ(20,*)(RDZL(K),K=1,KK)
      DO 235 K=1,KK
      DO 235 J=1,JJ
      DO 235 I=1,II
  235 DZ(I,J,K)=RDZL(K)
      DO 237 K=1,KK
  237 WRITE(IOCODE,513)K,RDZL(K)
      GO TO 245
232   WRITE(IOCODE,48)
      DO 240 K=1,KK
      WRITE(IOCODE,38)K
      DO 242 J=1,JJ
      READ(20,*)(DZ(I,J,K),I=1,II)
  242 WRITE(IOCODE,72)(DZ(I,J,K),I=1,II)
240   CONTINUE
245   CONTINUE
      WRITE(IOCODE,56)
      IF(KDZNET.GE.0) GO TO 255
      READ(20,*) DZNETC
      DO 250 K=1,KK
      DO 250 J=1,JJ
      DO 250 I=1,II
  250 DZNET(I,J,K)=DZNETC
      WRITE(IOCODE,56)
      WRITE(IOCODE,252) DZNETC
  252 FORMAT(T15,'GRID BLOCK NET THICKNESS (DZNET) IS INITIALLY',
     & ' SET AT',F10.4,' FOR ALL NODES',//)
      GO TO 274
  255 IF(KDZNET.GT.0) GO TO 268
      READ(20,*) (RDZNET(K),K=1,KK)
      DO 260 K=1,KK
      DO 260 J=1,JJ
      DO 260 I=1,II
  260 DZNET(I,J,K)=RDZNET(K)
      DO 266 K=1,KK
  266 WRITE(IOCODE,267) K,RDZNET(K)
  267 FORMAT(T15,'GRID SIZE (DZNET) IN LAYER ',I5,
     & ' IS INITIALLY SET AT',F8.2,' FOR ALL NODES',/)
      GO TO 274
  268 WRITE(IOCODE,269)
  269 FORMAT(//T15,8('*'),'GRID BLOCK NET THICKNESS ',
     & '(DZNET) DISTRIBUTION',8('*')/)
      DO 273 K=1,KK
       WRITE(IOCODE,38) K
      DO 270 J=1,JJ
      READ(20,*) (DZNET(I,J,K),I=1,II)
  270 WRITE(IOCODE,72) (DZNET(I,J,K),I=1,II)
  273 CONTINUE
  274 CONTINUE
      WRITE(IOCODE,56)
C********GRID BLOCK LENGTH MODIFICATIONS
      READ(20,69)
      READ(20,*) NUMDX,NUMDY,NUMDZ,NUMDZN,IDCODE
      IF(NUMDX.EQ.0) GO TO 8531
      WRITE(IOCODE,31)
      DO 275 L=1,NUMDX
      READ(20,*) I1,I2,J1,J2,K1,K2,REGVAL
      WRITE(IOCODE,32) I1,I2,J1,J2,K1,K2,REGVAL
      DO 275 K=K1,K2
      DO 275 J=J1,J2
      DO 275 I=I1,I2
      DX(I,J,K)=REGVAL
  275 CONTINUE
      IF(IDCODE.NE.1)GO TO 8531
      WRITE(IOCODE,43)
      DO 853 K=1,KK
       WRITE(IOCODE,38)K
      DO 854 J=1,JJ
  854 WRITE(IOCODE,72)(DX(I,J,K),I=1,II)
853   CONTINUE
8531  CONTINUE
      IF(NUMDY.EQ.0) GO TO 8551
      WRITE(IOCODE,34)
      DO 276 L=1,NUMDY
      READ(20,*) I1,I2,J1,J2,K1,K2,REGVAL
      WRITE(IOCODE,32) I1,I2,J1,J2,K1,K2,REGVAL
      DO 276 K=K1,K2
      DO 276 J=J1,J2
      DO 276 I=I1,I2
      DY(I,J,K)=REGVAL
  276 CONTINUE
      IF(IDCODE.NE.1) GO TO 8551
      WRITE(IOCODE,47)
      DO 855 K=1,KK
       WRITE(IOCODE,38)K
      DO 856 J=1,JJ
  856 WRITE(IOCODE,72)(DY(I,J,K),I=1,II)
855   CONTINUE
8551  CONTINUE
      IF(NUMDZ.EQ.0) GO TO 8571
      WRITE(IOCODE,37)
      DO 277 L=1,NUMDZ
      READ(20,*) I1,I2,J1,J2,K1,K2,REGVAL
      WRITE(IOCODE,32) I1,I2,J1,J2,K1,K2,REGVAL
      DO 277 K=K1,K2
      DO 277 J=J1,J2
      DO 277 I=I1,I2
      DZ(I,J,K)=REGVAL
  277 CONTINUE
      IF(IDCODE.NE.1) GO TO 8571
      WRITE(IOCODE,48)
      DO 857 K=1,KK
       WRITE(IOCODE,38)K
      DO 858 J=1,JJ
  858 WRITE(IOCODE,72)(DZ(I,J,K),I=1,II)
857   CONTINUE
8571  CONTINUE
      IF(NUMDZN.EQ.0) GO TO 890
      WRITE(IOCODE,860)
  860 FORMAT(//T15,8('*'),'GRID BLOCK NET THICKNESS',
     & ' (DZNET) NODE MODIFICATIONS',8('*'),//T15,
     & '   I1  I2  J1  J2  K1  K2  NEW DZNET VALUE')
      DO 865 L=1,NUMDZN
      READ(20,*) I1,I2,J1,J2,K1,K2,REGVAL
      WRITE(IOCODE,32) I1,I2,J1,J2,K1,K2,REGVAL
      DO 865 K=K1,K2
      DO 865 J=J1,J2
      DO 865 I=I1,I2
      DZNET(I,J,K)=REGVAL
  865 CONTINUE
      IF(IDCODE.NE.1) GO TO 890
      WRITE(IOCODE,269)
      DO 885 K=1,KK
       WRITE(IOCODE,38) K
      DO 880 J=1,JJ
  880 WRITE(IOCODE,72) (DZNET(I,J,K),I=1,II)
  885 CONTINUE
  890 CONTINUE
C      GRID BLOCK SIZE CHECK
      DO 900 K=1,KK
      DO 900 J=1,JJ
      DO 900 I=1,II
      IF(DX(I,J,K).GE.0.0) GO TO 892
      IFATAL=IFATAL+1
      WRITE(IOCODE,891) I,J,K
  891 FORMAT(/,5X,5('-'),'GRID BLOCK DX ERROR AT IJK =',3I5)
  892 IF(DY(I,J,K).GE.0.0) GO TO 894
      IFATAL=IFATAL+1
      WRITE(IOCODE,893) I,J,K
  893 FORMAT(/,5X,5('-'),'GRID BLOCK DY ERROR AT IJK =',3I5)
  894 IF(DZ(I,J,K).GE.0.0) GO TO 896
      IFATAL=IFATAL+1
      WRITE(IOCODE,895) I,J,K
  895 FORMAT(/,5X,5('-'),'GRID BLOCK DZ ERROR AT IJK =',3I5)
  896 IF(DZNET(I,J,K).GE.0.0) GO TO 900
      IFATAL=IFATAL+1
      WRITE(IOCODE,897) I,J,K
  897 FORMAT(/,5X,5('-'),'GRID BLOCK DZNET ERROR AT IJK =',3I5)
  900 CONTINUE
C*****ESTABLISH NODE MID-POINT ELEVATION.
      READ(20,69)
      READ(20,*)KEL
      IF(KEL.GT.1)GO TO 951
      IF(KEL.EQ.1)GO TO 920
      READ(20,*)ELEV
      DO 910 J=1,JJ
      DO 910 I=1,II
      VAREL(I,J)=ELEV
910   CONTINUE
  920 IF(KEL.NE.1)GO TO 930
      DO 922 J=1,JJ
      READ(20,*)(VAREL(I,J),I=1,II)
  922 CONTINUE
  930 CONTINUE
      DO 923 I=1,II
      DO 923 J=1,JJ
  923 SUM(I,J)=0.
      DO 926 K=1,KK
      DO 926 J=1,JJ
      DO 926 I=1,II
      IF(K.GE.2) GO TO 924
      EL(I,J,K)=VAREL(I,J)
      GO TO 926
  924 CONTINUE
      DEL=SUM(I,J)+DZ(I,J,K-1)
      EL(I,J,K)=VAREL(I,J)+DEL
      SUM(I,J)=DZ(I,J,K-1)+SUM(I,J)
  926 CONTINUE
      GO TO 601
  951 IF(KEL.EQ.3)GO TO 961
      READ(20,*)(RKEL(K),K=1,KK)
      DO 953 K=1,KK
      DO 953 J=1,JJ
      DO 953 I=1,II
  953 EL(I,J,K)=RKEL(K)
      GO TO 971
  961 DO 963 K=1,KK
      DO 965 J=1,JJ
  965 READ(20,*)(EL(I,J,K),I=1,II)
  963 CONTINUE
  971 CONTINUE
  601 WRITE(IOCODE,390)
      DO 600 K=1,KK
       WRITE(IOCODE,38)K
      DO 600 J=1,JJ
      WRITE(IOCODE,72)(EL(I,J,K),I=1,II)
600   CONTINUE
      DO 973 K=1,KK
      DO 973 J=1,JJ
      DO 973 I=1,II
  973 EL(I,J,K)=EL(I,J,K)+DZ(I,J,K)*0.5
69    FORMAT(40A2)
56    FORMAT(//)
   72 FORMAT(1X,15F8.0)
29    FORMAT(T15,'GRID BLOCK LENGTH (DX) IS INITIALLY',
     &' SET AT',F10.4,' FOR ALL NODES'//)
   31 FORMAT(//T15,'********GRID BLOCK LENGTH (DX) NODE MODIFICATIONS',
     & '**********',//T15,
     & '   I1  I2  J1  J2  K1  K2  NEW DX VALUE')
   32 FORMAT(15X,6I4,2X,E10.4)
33    FORMAT(T15,'GRID BLOCK WIDTH  (DY) IS INITIALLY',
     &' SET AT',F10.4,' FOR ALL NODES'//)
   34 FORMAT(//T15,'********GRID BLOCK WIDTH  (DY) NODE MODIFICATIONS',
     & '**********',//T15,
     & '   I1  I2  J1  J2  K1  K2  NEW DY VALUE')
36    FORMAT(T15,'GRID BLOCK GROSS THICKNESS (DZ) IS INITIALLY',
     &' SET AT',F10.4,' FOR ALL NODES'//)
   37 FORMAT(//T15,'********GRID BLOCK GROSS THICKNESS (DZ) NODE ',
     & 'MODIFICATIONS**********',//T15,
     & '   I1  I2  J1  J2  K1  K2  NEW DZ VALUE')
38    FORMAT(/1X,'K =',I2/)
390   FORMAT(//T15,'********** DEPTHS TO GRID BLOCK TOPS **********'/)
43    FORMAT(//T15,'********GRID BLOCK LENGTH (DX) DISTRIBUTION********'
     &/)
47    FORMAT(//T15,'********GRID BLOCK WIDTH  (DY) DISTRIBUTION********'
     &/)
   48 FORMAT(//T15,8('*'),'GRID BLOCK GROSS THICKNESS (DZ) ',
     & 'DISTRIBUTION',8('*'),/)
  511 FORMAT(T15,'GRID SIZE (DX) IN COLUMN',I5,' IS INITIALLY SET AT'
     &,F8.2,' FOR ALL NODES',/)
  512 FORMAT(T15,'GRID SIZE (DY) IN ROW   ',I5,' IS INITIALLY SET AT'
     &,F8.2,' FOR ALL NODES',/)
  513 FORMAT(T15,'GRID SIZE (DZ) IN LAYER ',I5,' IS INITIALLY SET AT'
     &,F8.2,' FOR ALL NODES',/)
      RETURN
      END
C.................................................................INTCOM
      SUBROUTINE INTCOM(IREG,X,Y,N,XO,YO)
C      MACHINE DEPENDENT INCLUDE STATEMENT
$INCLUDE:'PARAMS.FOR'
      DIMENSION X(LP10,LP9),Y(LP10,LP9)
      IF(XO.GE.X(IREG,N)) YO=Y(IREG,N)
      IF(XO.GE.X(IREG,N)) RETURN
      DO 10 I=2,N
      IF(XO.GT.X(IREG,I)) GO TO 10
      YO=Y(IREG,I)
      RETURN
  10  CONTINUE
      END
C.................................................................INTERP
      SUBROUTINE INTERP(IREG,X,Y,N,XO,YO)
C      MACHINE DEPENDENT INCLUDE STATEMENT
$INCLUDE:'PARAMS.FOR'
      DIMENSION X(LP10,LP9),Y(LP10,LP9)
      IF(XO.GE.X(IREG,N)) YO=Y(IREG,N)
      IF(XO.GE.X(IREG,N)) RETURN
      DO 10 I=2,N
      IF(XO.GE.X(IREG,I)) GO TO 10
      YO=Y(IREG,I-1) +(XO-X(IREG,I-1))*(Y(IREG,I)-Y(IREG,I-1))
     & /(X(IREG,I)-X(IREG,I-1))
      RETURN
  10  CONTINUE
      END
C.................................................................INTPVT
      SUBROUTINE INTPVT(IREG,BPT,RM,X,Y,N,XO,YO)
C      MACHINE DEPENDENT INCLUDE STATEMENT
$INCLUDE:'PARAMS.FOR'
      DIMENSION X(LP10,LP9),Y(LP10,LP9)
      IF(XO.GT.BPT)GO TO 100
      IF(XO.GE.X(IREG,N))YO=Y(IREG,N)
      IF(XO.GE.X(IREG,N))RETURN
      DO 10 I=2,N
      IF(XO.GE.X(IREG,I))GO TO 10
      YO=Y(IREG,I-1)+(XO-X(IREG,I-1))*(Y(IREG,I)-Y(IREG,I-1))
     & /(X(IREG,I)-X(IREG,I-1))
      RETURN
   10 CONTINUE
  100 CONTINUE
      DO 200 I=2,N
      IF(BPT.GE.X(IREG,I))GO TO 200
      YOBP=Y(IREG,I-1)+(BPT-X(IREG,I-1))*(Y(IREG,I)-Y(IREG,I-1))
     & /(X(IREG,I)-X(IREG,I-1))
      YO=YOBP+RM*(XO-BPT)
      RETURN
  200 CONTINUE
      END
C........................................................................LSORX
      SUBROUTINE LSORX(NX,NY,NZ,OMEGA,TOL,TOL1,MITER,DELT,DELT0,KSN,
     & N,NITER)
C      MACHINE DEPENDENT INCLUDE STATEMENT
$INCLUDE:'PARAMS.FOR'
C      ITERATIVE SOLUTION: X-DIRECTION TRIDIAGONAL ALGORITHM
      COMMON /COEF/ AW(LP1,LP2,LP3),AE(LP1,LP2,LP3),AN(LP1,LP2,LP3),
     & AS(LP1,LP2,LP3),AB(LP1,LP2,LP3),AT(LP1,LP2,LP3),E(LP1,LP2,LP3),
     & B(LP1,LP2,LP3)
      COMMON /SARRAY/ PN(LP1,LP2,LP3),IOCODE,IDMAX,
     & SON(LP1,LP2,LP3),SWN(LP1,LP2,LP3),SGN(LP1,LP2,LP3),
     & A1(LP1,LP2,LP3),A2(LP1,LP2,LP3),A3(LP1,LP2,LP3),
     & SUM(LP1,LP2,LP3),GAM(LP1,LP2,LP3),QS(LP1,LP2,LP3)
      COMMON /SPRTPS/ P(LP1,LP2,LP3),SO(LP1,LP2,LP3),SW(LP1,LP2,LP3),
     & SG(LP1,LP2,LP3)
      COMMON /VECTOR/ DX(LP1,LP2,LP3),DY(LP1,LP2,LP3),DZ(LP1,LP2,LP3),
     & DZNET(LP1,LP2,LP3),IQN1(LP11),IQN2(LP11),IQN3(LP11),IHEDIN(80)
      COMMON /TRIDI/ UM(LP15),AZL(LP15),BZL(LP15),CZL(LP15),DZL(LP15),
     & UZL(LP15),X(LP15),BETA(LP15),GAMMA(LP15),W(LP15)
      DIV=DELT/DELT0
      NITER=0
      DMAX=1.0
      RHO1=0.0
      THETA=0.0
   11 CONTINUE
      TW = 1.0 - OMEGA
      DMAX0=DMAX
      THETA0=THETA
      IF(NITER.GE.MITER)WRITE(IOCODE,30)NITER,TOL,DMAX
      IF(NITER.GE.MITER)RETURN
      NITER = NITER +1
      DMAX=0.0
      DO 20 K=1,NZ
      DO 20 J=1,NY
      DO 15 I=1,NX
      UM(I) = P(I,J,K)
      AZL(I) = AW(I,J,K)
      BZL(I)=E(I,J,K)
      CZL(I)=AE(I,J,K)
      DZL(I)=B(I,J,K)
      IF(NY.EQ.1) GO TO 14
      JM=J-1
      JP=J+1
      IF(J.EQ.1) JM=1
      IF(J.EQ.NY) JP=NY
      DZL(I)=DZL(I) - AS(I,J,K)*P(I,JM,K) - AN(I,J,K)*P(I,JP,K)
   14 IF(NZ.EQ.1) GO TO 15
      KM=K-1
      KP=K+1
      IF(K.EQ.1)KM=1
      IF(K.EQ.NZ) KP=NZ
      DZL(I)=DZL(I) - AT(I,J,K)*P(I,J,KM) - AB(I,J,K)*P(I,J,KP)
15    CONTINUE
C     THIS ALGORITHM SOLVES THE TRIDIAGONAL SYSTEM
C     GENERATED BY THE SYSTEM OF N EQUATIONS
C     AZL(I)*U(I-1) + BZL(I)*U(I) + CZL(I)*U(I+1) = DZL(I)
      BETA(1)=BZL(1)
      GAMMA(1)=DZL(1)/BZL(1)
      NXM=NX-1
C      COMPUTE FORWARD SOLUTION.
      DO 1010 ITRI=1,NXM
      W(ITRI)=CZL(ITRI)/BETA(ITRI)
      ITRIP=ITRI+1
 1010 BETA(ITRIP)=BZL(ITRIP)-AZL(ITRIP)*W(ITRI)
      DO 1020 ITRI=2,NX
      ITRIM=ITRI-1
 1020 GAMMA(ITRI)=(DZL(ITRI)-AZL(ITRI)*GAMMA(ITRIM))/BETA(ITRI)
C      COMPUTE BACK SOLUTION
      UZL(NX)=GAMMA(NX)
      DO 1030 JTRI=1,NXM
      ITRI=NX-JTRI
      ITRIP=ITRI+1
 1030 UZL(ITRI)=GAMMA(ITRI)-W(ITRI)*UZL(ITRIP)
      DO 16 I=1,NX
      GSLSOR=UZL(I)
      P(I,J,K) = TW*UM(I) + OMEGA*GSLSOR
      ARG=P(I,J,K) - UM(I)
      DM=ABS(ARG)
      IF(DM.GT.DMAX) DMAX = DM
16    CONTINUE
20    CONTINUE
      IF(TOL1.EQ.0.0)GO TO 25
      THETA=DMAX/DMAX0
      DELTA=THETA-THETA0
      ARG=DELTA
      ARG=ABS(ARG)
      IF(ARG.GT.TOL1)GO TO 25
      OM=OMEGA-1.0
      RHO1=(THETA+OM)*(THETA+OM)/(THETA*OMEGA*OMEGA)
      IF(RHO1.GE.1.0)GO TO 25
      ARG=1.0-RHO1
      OMEGA=2.0/(1.0+SQRT(ARG))
25    IF(DMAX.GT.TOL)GO TO 11
      IF(N.NE.KSN) GO TO 300
      WRITE(IOCODE,40)NITER,OMEGA,DMAX,THETA,RHO1
 300  CONTINUE
40    FORMAT(T5,'CONVERGENCE(LSORX) HAS BEEN REACHED AFTER ',I3,
     &' ITERATIONS',5X,'OMEGA = ',F6.3/
     &T5,'DMAX = ',F10.6,5X,'THETA = ',F10.6,5X,'RHO1 = ',F10.6/)
30    FORMAT(T15,'CONVERGENCE(LSORX) WAS NOT REACHED IN ',I5,
     &' ITERATIONS'/T15,'TOL = ',F10.7,10X,'DMAX = ',F15.7)
      RETURN
      END
C........................................................................LSORY
      SUBROUTINE LSORY(NX,NY,NZ,OMEGA,TOL,TOL1,MITER,DELT,DELT0,KSN,
     & N,NITER)
C      MACHINE DEPENDENT INCLUDE STATEMENT
$INCLUDE:'PARAMS.FOR'
C      ITERATIVE SOLUTION: Y-DIRECTION TRIDIAGONAL ALGORITHM
      COMMON /COEF/ AW(LP1,LP2,LP3),AE(LP1,LP2,LP3),AN(LP1,LP2,LP3),
     & AS(LP1,LP2,LP3),AB(LP1,LP2,LP3),AT(LP1,LP2,LP3),E(LP1,LP2,LP3),
     & B(LP1,LP2,LP3)
      COMMON /SARRAY/ PN(LP1,LP2,LP3),IOCODE,IDMAX,
     & SON(LP1,LP2,LP3),SWN(LP1,LP2,LP3),SGN(LP1,LP2,LP3),
     & A1(LP1,LP2,LP3),A2(LP1,LP2,LP3),A3(LP1,LP2,LP3),
     & SUM(LP1,LP2,LP3),GAM(LP1,LP2,LP3),QS(LP1,LP2,LP3)
      COMMON /SPRTPS/ P(LP1,LP2,LP3),SO(LP1,LP2,LP3),SW(LP1,LP2,LP3),
     & SG(LP1,LP2,LP3)
      COMMON /VECTOR/ DX(LP1,LP2,LP3),DY(LP1,LP2,LP3),DZ(LP1,LP2,LP3),
     & DZNET(LP1,LP2,LP3),IQN1(LP11),IQN2(LP11),IQN3(LP11),IHEDIN(80)
      COMMON /TRIDI/ UM(LP15),AZL(LP15),BZL(LP15),CZL(LP15),DZL(LP15),
     & UZL(LP15),X(LP15),BETA(LP15),GAMMA(LP15),W(LP15)
      DIV=DELT/DELT0
      NITER=0
      DMAX=1.0
      RHO1=0.0
      THETA=0.0
   11 CONTINUE
      TW = 1.0 - OMEGA
      DMAX0=DMAX
      THETA0=THETA
      IF(NITER.GE.MITER)WRITE(IOCODE,30)NITER,TOL,DMAX
      IF(NITER.GE.MITER)RETURN
      NITER = NITER +1
      DMAX=0.0
      DO 20 K=1,NZ
      DO 20 I=1,NX
      DO 15 J=1,NY
      UM(J) = P(I,J,K)
      AZL(J) = AS(I,J,K)
      BZL(J)=E(I,J,K)
      CZL(J)=AN(I,J,K)
      DZL(J)=B(I,J,K)
      IF(NX.EQ.1) GO TO 14
      IM=I-1
      IP=I+1
      IF(I.EQ.1) IM=1
      IF(I.EQ.NX) IP=NX
      DZL(J)=DZL(J) - AW(I,J,K)*P(IM,J,K) - AE(I,J,K)*P(IP,J,K)
   14 IF(NZ.EQ.1) GO TO 15
      KM=K-1
      KP=K+1
      IF(K.EQ.1)KM=1
      IF(K.EQ.NZ) KP=NZ
      DZL(J)=DZL(J) - AT(I,J,K)*P(I,J,KM) - AB(I,J,K)*P(I,J,KP)
15    CONTINUE
C     THIS ALGORITHM SOLVES THE TRIDIAGONAL SYSTEM
C     GENERATED BY THE SYSTEM OF N EQUATIONS
C     AZL(J)*U(J-1) + BZL(J)*U(J) + CZL(J)*U(J+1) = DZL(J)
      BETA(1)=BZL(1)
      GAMMA(1)=DZL(1)/BZL(1)
      NYM=NY-1
C      COMPUTE FORWARD SOLUTION.
      DO 1010 ITRI=1,NYM
      W(ITRI)=CZL(ITRI)/BETA(ITRI)
      ITRIP=ITRI+1
 1010 BETA(ITRIP)=BZL(ITRIP)-AZL(ITRIP)*W(ITRI)
      DO 1020 ITRI=2,NY
      ITRIM=ITRI-1
 1020 GAMMA(ITRI)=(DZL(ITRI)-AZL(ITRI)*GAMMA(ITRIM))/BETA(ITRI)
C      COMPUTE BACK SOLUTION
      UZL(NY)=GAMMA(NY)
      DO 1030 JTRI=1,NYM
      ITRI=NY-JTRI
      ITRIP=ITRI+1
 1030 UZL(ITRI)=GAMMA(ITRI)-W(ITRI)*UZL(ITRIP)
      DO 16 J=1,NY
      GSLSOR=UZL(J)
      P(I,J,K) = TW*UM(J) + OMEGA*GSLSOR
      ARG=P(I,J,K) - UM(J)
      DM=ABS(ARG)
      IF(DM.GT.DMAX) DMAX = DM
16    CONTINUE
20    CONTINUE
      IF(TOL1.EQ.0.0)GO TO 25
      THETA=DMAX/DMAX0
      DELTA=THETA-THETA0
      ARG=DELTA
      ARG=ABS(ARG)
      IF(ARG.GT.TOL1)GO TO 25
      OM=OMEGA-1.0
      RHO1=(THETA+OM)*(THETA+OM)/(THETA*OMEGA*OMEGA)
      IF(RHO1.GE.1.0)GO TO 25
      ARG=1.0-RHO1
      OMEGA=2.0/(1.0+SQRT(ARG))
25    IF(DMAX.GT.TOL)GO TO 11
      IF(N.NE.KSN) GO TO 300
      WRITE(IOCODE,40)NITER,OMEGA,DMAX,THETA,RHO1
 300  CONTINUE
40    FORMAT(T5,'CONVERGENCE(LSORY) HAS BEEN REACHED AFTER ',I3,
     &' ITERATIONS',5X,'OMEGA = ',F6.3/
     &T5,'DMAX = ',F10.6,5X,'THETA = ',F10.6,5X,'RHO1 = ',F10.6/)
30    FORMAT(T15,'CONVERGENCE(LSORY) WAS NOT REACHED IN ',I5,
     &' ITERATIONS'/T15,'TOL = ',F10.7,10X,'DMAX = ',F15.7)
      RETURN
      END
C.................................................................LSORZ
      SUBROUTINE LSORZ(NX,NY,NZ,OMEGA,TOL,TOL1,MITER,DELT,DELT0,KSN,
     & N,NITER)
C      MACHINE DEPENDENT INCLUDE STATEMENT
$INCLUDE:'PARAMS.FOR'
C      ITERATIVE SOLUTION: Z-DIRECTION TRIDIAGONAL ALGORITHM
      COMMON /COEF/ AW(LP1,LP2,LP3),AE(LP1,LP2,LP3),AN(LP1,LP2,LP3),
     & AS(LP1,LP2,LP3),AB(LP1,LP2,LP3),AT(LP1,LP2,LP3),E(LP1,LP2,LP3),
     & B(LP1,LP2,LP3)
      COMMON /SARRAY/ PN(LP1,LP2,LP3),IOCODE,IDMAX,
     & SON(LP1,LP2,LP3),SWN(LP1,LP2,LP3),SGN(LP1,LP2,LP3),
     & A1(LP1,LP2,LP3),A2(LP1,LP2,LP3),A3(LP1,LP2,LP3),
     & SUM(LP1,LP2,LP3),GAM(LP1,LP2,LP3),QS(LP1,LP2,LP3)
      COMMON /SPRTPS/ P(LP1,LP2,LP3),SO(LP1,LP2,LP3),SW(LP1,LP2,LP3),
     & SG(LP1,LP2,LP3)
      COMMON /VECTOR/ DX(LP1,LP2,LP3),DY(LP1,LP2,LP3),DZ(LP1,LP2,LP3),
     & DZNET(LP1,LP2,LP3),IQN1(LP11),IQN2(LP11),IQN3(LP11),IHEDIN(80)
      COMMON /TRIDI/ UM(LP15),AZL(LP15),BZL(LP15),CZL(LP15),DZL(LP15),
     & UZL(LP15),X(LP15),BETA(LP15),GAMMA(LP15),W(LP15)
      DIV=DELT/DELT0
      NITER=0
      DMAX=1.0
      RHO1=0.0
      THETA=0.0
   11 CONTINUE
      TW = 1.0 - OMEGA
      DMAX0=DMAX
      THETA0=THETA
      IF(NITER.GE.MITER)WRITE(IOCODE,30)NITER,TOL,DMAX
      IF(NITER.GE.MITER)RETURN
      NITER = NITER +1
      DMAX=0.0
      DO 20 I=1,NX
      DO 20 J=1,NY
      DO 15 K=1,NZ
      UM(K) = P(I,J,K)
      AZL(K) = AT(I,J,K)
      BZL(K)=E(I,J,K)
      CZL(K)=AB(I,J,K)
      DZL(K)=B(I,J,K)
      IF(NY.EQ.1) GO TO 14
      JM=J-1
      JP=J+1
      IF(J.EQ.1) JM=1
      IF(J.EQ.NY) JP=NY
      DZL(K)=DZL(K) - AS(I,J,K)*P(I,JM,K) - AN(I,J,K)*P(I,JP,K)
   14 IF(NX.EQ.1) GO TO 15
      IM=I-1
      IP=I+1
      IF(I.EQ.1)IM=1
      IF(I.EQ.NX) IP=NX
      DZL(K)=DZL(K) - AW(I,J,K)*P(IM,J,K) - AE(I,J,K)*P(IP,J,K)
   15 CONTINUE
C     THIS ALGORITHM SOLVES THE TRIDIAGONAL SYSTEM
C     GENERATED BY THE SYSTEM OF N EQUATIONS
C     AZL(K)*U(K-1) + BZL(K)*U(K) + CZL(K)*U(K+1) = DZL(K)
      BETA(1)=BZL(1)
      GAMMA(1)=DZL(1)/BZL(1)
      NZM=NZ-1
C      COMPUTE FORWARD SOLUTION.
      DO 1010 ITRI=1,NZM
      W(ITRI)=CZL(ITRI)/BETA(ITRI)
      ITRIP=ITRI+1
 1010 BETA(ITRIP)=BZL(ITRIP)-AZL(ITRIP)*W(ITRI)
      DO 1020 ITRI=2,NZ
      ITRIM=ITRI-1
 1020 GAMMA(ITRI)=(DZL(ITRI)-AZL(ITRI)*GAMMA(ITRIM))/BETA(ITRI)
C      COMPUTE BACK SOLUTION
      UZL(NZ)=GAMMA(NZ)
      DO 1030 JTRI=1,NZM
      ITRI=NZ-JTRI
      ITRIP=ITRI+1
 1030 UZL(ITRI)=GAMMA(ITRI)-W(ITRI)*UZL(ITRIP)
      DO 16 K=1,NZ
      GSLSOR=UZL(K)
      P(I,J,K) = TW*UM(K) + OMEGA*GSLSOR
      ARG=P(I,J,K) - UM(K)
      DM=ABS(ARG)
      IF(DM.GT.DMAX) DMAX = DM
16    CONTINUE
20    CONTINUE
      IF(TOL1.EQ.0.0)GO TO 25
      THETA=DMAX/DMAX0
      DELTA=THETA-THETA0
      ARG=DELTA
      ARG=ABS(ARG)
      IF(ARG.GT.TOL1)GO TO 25
      OM=OMEGA-1.0
      RHO1=(THETA+OM)*(THETA+OM)/(THETA*OMEGA*OMEGA)
      IF(RHO1.GE.1.0)GO TO 25
      ARG=1.0-RHO1
      OMEGA=2.0/(1.0+SQRT(ARG))
25    IF(DMAX.GT.TOL)GO TO 11
      IF(N.NE.KSN) GO TO 300
      WRITE(IOCODE,40)NITER,OMEGA,DMAX,THETA,RHO1
 300  CONTINUE
40    FORMAT(T5,'CONVERGENCE(LSORZ) HAS BEEN REACHED AFTER ',I3,
     &' ITERATIONS',5X,'OMEGA = ',F6.3/
     &T5,'DMAX = ',F10.6,5X,'THETA = ',F10.6,5X,'RHO1 = ',F10.6/)
30    FORMAT(T15,'CONVERGENCE(LSORZ) WAS NOT REACHED IN ',I5,
     &' ITERATIONS'/T15,'TOL = ',F10.7,10X,'DMAX = ',F15.7)
      RETURN
      END
