SUBROUTINE CONOPT(M, N, MM, NN, MN, A, NET, IPRINT, NP, IJ)       CON   10
C THIS SUBPROGRAM OPTIMIZES THE CONTOUR SCAN AND IT RETURNS
C THE OPTIMUM CONTOUR SCANNING PATH STORED IN THE ARRAY IJ.
C                   .....INPUT VARIABLES.....
C M  -NUMBER OF SUBDIVISIONS IN X-DIRECTION. DIMENSION OF THE
C     ARRAY A IN I DIRECTION.
C N  -AS ABOVE BUT IN Y AND J DIRECTIONS.
C A  -WORKING ARRAY WITH NUMERICAL VALUES OF DIRICHLET
C     BOUNDARY CONDITIONS.  A(M,N).
C MM -NUMBER OF SUBDIVISIONS IN X-DIRECTION. DIMENSION OF THE
C     ARRAY NET IN I DIRECTION. MM=M+2.
C NN -AS ABOVE BUT IN Y AND J DIRECTIONS. NN=N+2.
C NET-AN INTEGER ARRAY WITH FLAGS LOCATING THE DIRICHLET
C     AND NEUMANN BOUNDARY CONDITIONS. 0=VARIABLES,
C     1=DIRICHLET B.C., 2=NEUMANN B.C.  NET(MM,NN). THESE
C     FLAGS ARE SHIFTED WITH RESPECT TO THE ARRAY A SUCH
C     THAT THE LOCATION OF A POINT ON THE DIRICHLET BOUNDARY
C     A(I,J) CORRESPONDS TO NET(I+1,J+1).
C MN -DIMENSION OF THE ARRAY IJ. IT DEPENDS ON THE LENGTH
C     OF THE CONTOUR PATH, NR. PREFERABLY MN=M*N.
C IPRINT -IF IPRINT=0, THE PRINTOUT IN CONOPT IS SUPPRESSED.
C                   .....OUTPUT VARIABLES.....
C NP -NUMBER OF POINTS SCANNED BY THE CONTOUR PATH DETERMINED
C     BY THE SUBROUTINE CONTR.
C IJ -AN INTEGER ARRAY CONTAINING COORDINATES I,J OF THE
C     PONTS (I,J)  SCANNED BY THE CONTOUR PATH. THEY ARE
C     DETERMINED BY CONTR. THE ARRAY IJ HAS 3 COLUMNS, THE
C     LAST ONE IS USED FOR COMPUTATIONAL PURPOSES IN CONTR,
C     AND, ON THE OUTPUT, IT CONTAINS THE CORRESPONDING FLAGS
C     FROM THE ARRAY NET.
C                   .....INTERNAL VARIABLES.....
C K  -NUMBER OF STARTING POINTS DETERMINED BY THE SUBROUTINE
C     SPOINT. K IS USED ONLY IN THIS SUBROUTINE.
C IJS-AN INTEGER ARRAY. IT CONTAINS 4 COLUMNS, THE FIRST TWO
C     OF THEM CONTAIN COORDINATES I,J OF THE STARTING POINTS
C     (I,J) DETERMINED BY THE SUBROUTINE SPOINT, THE THIRD
C     ONE CONTAINS CODED VALUES OF DIRICHLET B.C. ADJACENT
C     TO A POINT (I,J), AND THE LAST ONE CONTAINS THE NUMBER
C     OF BREAKS IN THE CONTOUR SCAN WHICH IS INITIATED FROM
C     THE CORRESPONDING POINT (I,J). THE NUMBER OF ROWS IN
C     IJS CAN BE CHANGED, IF NECESSARY.
C KK -NUMBER OF HIGH VALUE DIRICHLET BOUNDARY CONDITIONS,
C     DETERMINED BY THE SUBROUTINE HIVA. KK IS USED IN SPOINT
C AA -AN ARRAY CONTAINING THE HIGH VALUE DIRICHLET B.C.,
C     FOUND BY HIVA. THE LENGTH OF THE ARRAY AA CAN BE
C     CHANGED, IF NECESSARY.
C                   .....NOTES.....
C THE PARAMETERS M,N,MM,NN,A,NET AND IPRINT MUST BE SUPPLIED
C BY THE USER (INPUT). PARAMETERS NP AND IJ ARE RETURNED BY
C THE SUBROUTINE CONOPT (OUTPUT).  THE MAIN PROGRAM COMMUNI-
C CATES WITH THE SUBPROGRAM CONOPT THROUGH THE PARAMETER LIST
C ONLY.  ALL THE LABELED COMMON STORAGE WITHIN THE SUBPROG-
C RAM CONOPT IS FOR THE INTERNAL USE ONLY.
C IF THE CONDITIONS MM= M+2 AND NN= N+2 ARE NOT SATISFIED,
C EXECUTION OF THE SUBPROGRAM IS INTERRUPTED AND A MESSAGE IS
C PRINTED OUT  (AUTOMATIC ADJUSTMENT OF MM AND NN MIGHT CAUSE
C SEVERE ERRORS IN THE MAIN PROGRAM).  IF THERE IS NO DIRICH-
C LET BOUNDARY POINT (NO C1), THE COMPUTER AUTOMATICALLY IN-
C SERTS SUCH A POINT IN THE REGION R1 WITHIN THE BOUNDARY C2
C (THIS MERELY PINS DOWN ONE SOLUTION OF THE NEUMANN PROBLEM).
C IF THERE ARE MORE INTERIOR POINTS THAN ALLOWED BY THE IN-
C PUT PARAMETER MN, EXECUTION IS INTERRUPTED AND A MESSAGE
C IS PRINTED OUT  (MN MAY NOT BE ADJUSTED AUTOMATICALLY BY
C THE COMPUTER).
      DIMENSION A(M,N), NET(MM,NN), IJ(MN,3)
      COMMON /SPT/ K, IJS(20,4) /CS/ NR /HV/ KK, AA(20)
      IF (MM.NE.M+2 .OR. NN.NE.N+2) GO TO 70
      IP = IPRINT
      IF (IP.NE.0) WRITE (6,99999)
C FIND PROPER STARTING POINTS FOR THE CONTOUR SCAN.
      CALL SPOINT(M, N, MM, NN, A, NET, IP)
C FIND A CONTINUOUS CONTOUR PATH, IF POSSIBLE.
      L = 0
   10 L = L + 1
      CALL CONTR(MM, NN, MN, NET, IJ, L)
      NP = NR
      IF (IJS(L,4).EQ.0) GO TO 50
      IF (L.LT.K) GO TO 10
C SINCE THERE IS NO CONTINUOUS PATH FIND THE PATH WITH THE
C LARGEST NUMBER OF BREAKS.
      IG = 0
      DO 20 L=1,K
        IF (IJS(L,4).GT.IG) IG = IJS(L,4)
   20 CONTINUE
C FIND THE PATH WITH THE SMALLEST NUMBER OF BREAKS.
      IS = IG
      DO 30 L=1,K
        IF (IJS(L,4).LT.IS) IS = IJS(L,4)
   30 CONTINUE
C CHOOSE A PATH WITH THE SMALLEST NUMBER OF BREAKS BEING
C CLOSE TO THE HIGHEST VALUE OF DIRICHLET B.C.
      L = 0
   40 L = L + 1
      IF (IJS(L,4).EQ.IS) GO TO 60
      IF (L.LT.K) GO TO 40
   50 I = IJS(L,3)
      IF (IP.NE.0) WRITE (6,99998) L, (IJS(L,J),J=1,2), AA(I)
      RETURN
   60 I = IJS(L,3)
      IF (IP.NE.0) WRITE (6,99997) IS, AA(I), (IJS(L,J),J=1,2)
C DETERMINE THE FINAL CONTOUR PATH.
      CALL CONTR(MM, NN, MN, NET, IJ, L)
      RETURN
   70 WRITE (6,99996)
      STOP
99999 FORMAT (1H1, 10X, 26HCONTOUR PATH OPTIMIZATION.//)
99998 FORMAT (/1X, 39HTHE CONTINUOUS CONTOUR SCAN HAS BEEN FO,
     * 9HUND AFTER, I3, 15H EVALUATION(S)./17H THE STARTING POI,
     * 8HNT IS A(, I2, 1H,, I2, 4H) AT, E12.4, 5H B.C.//)
99997 FORMAT (/1X, 28HTHE OPTIMUM CONTOUR SCAN HAS, I4, 2H B,
     * 17HREAK(S) (MINIMUM)/30H AND IT STARTS FROM THE HIGHES,
     * 17HT POSSIBLE VALUE,, E12.4/25H THE STARTING POINT IS A(,
     * I2, 1H,, I2, 2H).//)
99996 FORMAT (1H1, 1X, 36HADJUST THE PARAMETERS MM AND NN SUCH,
     * 5H THAT/1X, 8HMM = M+2/1X, 8HNN = N+2//1X, 14HSTOP FROM CONO,
     * 3HPT.)
      END
      SUBROUTINE SPOINT(M, N, MM, NN, A, NET, IP)                       SPO   10
C THE SUBROUTINE DETERMINES THE POSITION OF POINTS IN THE
C VICINITY OF THE HIGHEST AND INTERMEDIATE VALUES OF DIRICHLET
C BOUNDARY CONDITIONS. THE ARRAYS A AND NET SERVE AS AN INPUT.
C ML -AN INTEGER, ML=M+1
C NL -AN INTEGER, NL=N+1
C THE PARAMETERS ARE DEFINED IN CONOPT.
      DIMENSION A(M,N), NET(MM,NN)
      COMMON /COEF/ ML, NL /SPT/ K, IJS(20,4) /HV/ KK, AA(20)
C FIND HIGH VALUE DIRICHLET B.C.
      CALL HIVA(M, N, MM, NN, A, NET, IP)
      K = 0
      DO 100 L=1,KK
        AM = AA(L)
C FIND POINTS SEPARATED FROM THE HIGH VALUE DIRICHLET B.C.
C BY ONE MESH LENGTH. THESE POINTS LIE ON THE ROWS.
        DO 30 I=1,M
          DO 20 J=1,N
            IF (NET(I+1,J+1).NE.0) GO TO 20
            IF (A(I-1,J).EQ.AM .AND. NET(I,J).NE.1) GO TO 10
            IF (A(I+1,J).NE.AM .OR. NET(I+2,J+2).EQ.1) GO TO 20
   10       K = K + 1
            IF (K.GT.20) GO TO 110
            IJS(K,1) = I
            IJS(K,2) = J
            IJS(K,3) = L
   20     CONTINUE
   30   CONTINUE
C FIND POINTS AS ABOVE BUT LYING ON THE COLUMNS.
        DO 60 J=1,N
          DO 50 I=1,M
            IF (NET(I+1,J+1).NE.0) GO TO 50
            IF (A(I,J+1).EQ.AM .AND. NET(I,J+2).NE.1) GO TO 40
            IF (A(I,J-1).NE.AM .OR. NET(I+2,J).EQ.1) GO TO 50
   40       K = K + 1
            IF (K.GT.20) GO TO 110
            IJS(K,1) = I
            IJS(K,2) = J
            IJS(K,3) = L
   50     CONTINUE
   60   CONTINUE
C FIND POINTS ADJACENT TO DIRICHLET B.C., LOCATED AT CORNERS
        DO 90 I=2,ML
          DO 80 J=2,NL
            IF (NET(I,J).NE.0) GO TO 80
            IF (A(I-2,J-1).EQ.AM .AND. NET(I-1,J-1).EQ.1 .AND.
     *       NET(I,J-1).EQ.1) GO TO 70
            IF (A(I-1,J).EQ.AM .AND. NET(I-1,J+1).EQ.1 .AND.
     *       NET(I-1,J).EQ.1) GO TO 70
            IF (A(I,J-1).EQ.AM .AND. NET(I+1,J+1).EQ.1 .AND.
     *       NET(I,J+1).EQ.1) GO TO 70
            IF (A(I-1,J-2).NE.AM .OR. NET(I+1,J-1).NE.1 .OR.
     *       NET(I+1,J).NE.1) GO TO 80
   70       K = K + 1
            IF (K.GT.20) GO TO 110
            IJS(K,1) = I - 1
            IJS(K,2) = J - 1
            IJS(K,3) = L
   80     CONTINUE
   90   CONTINUE
  100 CONTINUE
  110 IF(IP.NE.0) WRITE(6,99999) (I,(IJS(I,J), J=1,2), I=1,K)
      RETURN
99999 FORMAT (/1X, 39HSTARTING POINTS FOR CONTOUR OPTIMIZATIO, 1HN//
     * 5X, 2HNO, 5X, 1HI, 5X, 1HJ/(I7, 2I6))
      END
      SUBROUTINE HIVA(M, N, MM, NN, A, NET, IP)                         HIV   10
C THE SUBROUTINE DETERMINES THE HIGH VALUE DIRICHLET BOUNDA-
C RY CONDITIONS AND IT RETURNS THEM IN THE ARRAY AA.
C PARAMETERS ARE DEFINED IN SPOINT.
      DIMENSION A(M,N), NET(MM,NN)
      COMMON /HV/ KK, AA(20)
C FIND THE MAXIMUM AND MINIMUM VALUES IN THE ARRAY A.
C ARE THERE ANY POINTS OF TYPE S1
      KK = 0
      DO 20 J=1,N
        J1 = J + 1
        DO 10 I=1,M
          IF (NET(I+1,J1).EQ.1) KK = KK + 1
   10   CONTINUE
   20 CONTINUE
C IF KK=0 (NO S1 POINTS), THE SOLUTION WILL BE PINNED DOWN.
      IF (KK) 30, 30, 60
   30 DO 50 J=1,N
        J1 = J + 1
        DO 40 I=1,M
          IF (NET(I+1,J1).NE.0) GO TO 40
          NET(I+1,J1) = 1
          A(I,J) = 50.0
          AA(1) = 50.0
          KK = 1
          GO TO 120
   40   CONTINUE
   50 CONTINUE
      WRITE (6,99999)
      STOP
   60 AM = 0.0
      AS = 100.0
      DO 80 J=1,N
        J1 = J + 1
        DO 70 I=1,M
          IF (NET(I+1,J1).NE.1) GO TO 70
          IF (A(I,J).GT.AM) AM = A(I,J)
          IF (A(I,J).LT.AS) AS = A(I,J)
   70   CONTINUE
   80 CONTINUE
      KK = 1
      AA(1) = AM
C FIND INTERMEDIATE VALUES.
      DO 110 K=1,19
        AM = AS
        DO 100 J=1,N
          J1 = J + 1
          DO 90 I=1,M
            IF (NET(I+1,J1).NE.1) GO TO 90
            IF (A(I,J).GT.AM .AND. A(I,J).LT.AA(K)) AM = A(I,J)
   90     CONTINUE
  100   CONTINUE
        IF (AM.EQ.AS) GO TO 120
        KK = KK + 1
        AA(KK) = AM
  110 CONTINUE
  120 IF (IP.NE.0) WRITE (6,99998) (I,AA(I),I=1,KK)
      RETURN
99999 FORMAT (/1X, 37HTHERE ARE NO FREE POINTS WITHIN A AND,
     * 26H NET. SOLUTION IMPOSSIBLE.//1X, 14HSTOP FROM HIVA,
     * 11H IN CONOPT.)
99998 FORMAT (/1X, 25HHIGH VALUE DIRICHLET B.C./(I7, E15.5))
      END
      SUBROUTINE CONTR(MM, NN, MN, NET, IJ, L1)                         CON   10
C THE SUBROUTINE CONTR TRACES A CONTOUR SCAN WHICH COVERS
C ONLY THE INTERIOR OF THE REGION DEFINED BY BOUNDARY OF AN
C ARBITRARY SHAPE. THE SUBROUTINE REQUIRES A STARTING POINT,
C I,J,AND IT RETURNS THE CONTOUR PATH STORED IN THE ARRAY IJ
C IX -AN ARRAY USED ONLY IN THIS SUBROUTINE FOR SETTING
C     A DIRECTIONAL PATTERN OF I.
C JY -AS ABOVE, BUT FOR J.
C L  -AN INDICATOR OF THE PERFORMED TURN. L=1-TURN LEFT,
C     L=2-GO STRAIGHT, L=3-TURN RIGHT.
C K  -AN INDICATOR OF THE CHOSEN DIRECTION.K IS RELATED TO L
C OTHER PARAMETERS AS DEFINED IN CONOPT.
      DIMENSION NET(MM,NN), IJ(MN,3)
      COMMON /COEF/ ML, NL /TEMP/ IX(5), JY(5)
      COMMON /SPT/ K1, IJS(20,4) /CS/ NR
C SET INITIAL PARAMETERS.
      I = IJS(L1,1) + 1
      J = IJS(L1,2) + 1
      L = 2
      NR = 0
      IJS(L1,4) = 0
      GO TO 30
C FIND INITIAL DIRECTION OF THE SCAN.
   10 IX(5) = IX(1)
      JY(5) = JY(1)
      DO 20 K=1,4
        I = IX(K+1)
        J = JY(K+1)
        I1 = IX(K)
        J1 = JY(K)
        IF (NET(I,J).NE.1 .AND. NET(I1,J1).EQ.1) GO TO 30
   20 CONTINUE
   30 NR = NR + 1
      IF (NR.GT.MN) GO TO 110
C SAVE I,J, AND NET(I,J). PROTECT NET(I,J).
      IJ(NR,1) = I - 1
      IJ(NR,2) = J - 1
      IJ(NR,3) = NET(I,J)
      NET(I,J) = 1
C SET THE DIRECTION PATTERN.
      IX(1) = I - 1
      JY(1) = J
      IX(2) = I
      JY(2) = J + 1
      IX(3) = I + 1
      JY(3) = J
      IX(4) = I
      JY(4) = J - 1
      IF (NR.EQ.1) GO TO 10
C COMPUTE THE RELATIVE DIRECTION.
      K2 = K + L - 2
      K = K2 + 4*(1-((3+K2)/4))
      K2 = K - 1
      IF (K2.EQ.0) GO TO 60
C ROTATE THE PATTERN ACCORDING TO THE RELATIVE DIRECTION.
      DO 50 L2=1,K2
        IX(5) = IX(1)
        JY(5) = JY(1)
        DO 40 L=1,4
          IX(L) = IX(L+1)
          JY(L) = JY(L+1)
   40   CONTINUE
   50 CONTINUE
C CHECK ALL THE POSSIBLE MOVES. IF A MOVE IS POSSIBLE, THE
C VALUE OF L WILL BE PRESERVED.
   60 DO 70 L=1,3
        I = IX(L)
        J = JY(L)
        IF (NET(I,J).NE.1) GO TO 30
   70 CONTINUE
C FIND AND COUNT POSSIBLE BREAKS IN THE SCAN.
      K = 1
      L = 2
      DO 90 I=2,ML
        DO 80 J=2,NL
          IF (NET(I,J).NE.1) IJS(L1,4) = IJS(L1,4) + 1
          IF (NET(I,J).NE.1) GO TO 30
   80   CONTINUE
   90 CONTINUE
C RESET THE ARRAY NET
      DO 100 K=1,NR
        I = IJ(K,1) + 1
        J = IJ(K,2) + 1
        NET(I,J) = IJ(K,3)
  100 CONTINUE
      RETURN
  110 WRITE (6,99999)
      STOP
99999 FORMAT (/1X, 35HTHERE ARE MORE INTERIOR POINTS THAN, 6H ALLOW,
     * 9HED BY MN./1X, 19HINCREASE MN TO M*N./1X, 14HSTOP FROM CONT,
     * 12HR IN CONOPT.)
      END
      SUBROUTINE RASTER(MM, NN, MN, NET, IJ)                            RAS   10
C THE SUBROUTINE RETURNS A RASTER SCAN OF THE INTERIOR OF AN MXN
C ARRAY. ONLY ACTIVE POINTS OF THE SCAN ARE SAVED IN IJ.
      DIMENSION NET(MM,NN), IJ(MN,3)
      COMMON /COEF/ ML, NL /CS/ NR
      WRITE (6,99999)
      NR = 0
      DO 20 I=2,ML
        DO 10 J=2,NL
          IF (NET(I,J).EQ.1) GO TO 10
          NR = NR + 1
          IJ(NR,1) = I - 1
          IJ(NR,2) = J - 1
   10   CONTINUE
   20 CONTINUE
      RETURN
99999 FORMAT (1H1, 10X, 11HRASTER SCAN)
      END
      SUBROUTINE SPIRAL(M, N, MM, NN, MN, NET, IJ)                      SPI   10
C THE SUBROUTINE TRACES AN INWARD SPIRAL PATH WITHIN AN
C M X N ARRAY A. THE SCAN IS RETURNED IN IJ.
C THE PARAMETERS ARE DEFINED IN THE SUBROUTINE CONOPT AND SPOINT.
      DIMENSION NET(MM,NN), IJ(MN,3)
      COMMON /COEF/ ML, NL /CS/ NR
      WRITE (6,99999)
      NR = 0
C FIND THE NUMBER OF RECTANGLES WHICH WILL CONSTITUTE
C THE SPIRAL SCAN.
      IH = MIN0(ML,NL)/2
C BEGIN FROM THE LARGEST RECTANGLE AND SUCCESSIVELY REDUCE
C ITS SIZE.
      DO 90 I1=1,IH
        IS = I1
        K = 0
C LOCATE START AND END POINTS FOR EACH SIDE OF A RECTANGLE.
        DO 80 I2=1,2
          I = IS
          IF (I2.EQ.2) I = IE + 1
          IE = N - IS
          DO 70 I3=1,2
            K = K + 1
C SCAN ONE SIDE OF A RECTANLE AND DISTRIBUTE THE SEQUENCE OF
C POINTS TO I OR TO J.
            DO 60 I4=IS,IE
              GO TO (10, 20, 30, 40), K
   10         J = I4
              GO TO 50
   20         I = I4
              GO TO 50
   30         J = NL - I4
              GO TO 50
   40         I = ML - I4
   50         IF (NET(I+1,J+1).EQ.1) GO TO 60
              NR = NR + 1
              IJ(NR,1) = I
              IJ(NR,2) = J
   60       CONTINUE
            J = IE + 1
            IF (I2.EQ.2) J = IS
            IE = M - IS
   70     CONTINUE
   80   CONTINUE
   90 CONTINUE
      RETURN
99999 FORMAT (1H1, 10X, 11HSPIRAL SCAN)
      END
      SUBROUTINE SORCAR(M, N, MN, A, IJ, C, IP, ER)                     SOR   10
C THIS SUBROUTINE OPTIMIZES THE OVERRELAXATION FACTOR USING CARRE*S
C TECHNIQUE.
      DIMENSION A(M,N), IJ(MN,3), C(MN), GX(3)
      COMMON /LOG/ Q, QQ
      LOGICAL Q, QQ, QQQ
      WRITE (6,99999)
C THE FIRST ITERATION USES BETA=1.0
      BETA = 1.
C THE FIRST INTERVAL OF INT ITERATIONS USES BETA=1.375
      BOPT = 1.375
C SET THE INTERVAL OF ITERATIONS
      INT = 12
      Q = .TRUE.
      QQQ = .FALSE.
      IN = 0
      SUM = 0.0
      JOB = 0
   10 JOB = JOB + 1
      IN = IN + 1
      QQ = IN.LT.INT
      BIG = 0.0
      S = SUM
      SUM = 0.0
      CALL SOR(M, N, MN, A, IJ, C, BETA, SUM, BIG)
      QQ = (.NOT.QQ) .AND. QQQ
C Q=.FALSE. IF BETA OPTIMUM IS FOUND.
      IF (.NOT.Q) GO TO 60
C QQ=.TRUE. IF JOB.GT.INT+1 AND IN=INT
      IF (QQ) GO TO 30
C QQQ=.TRUE. IF JOB.GT.INT+1
      IF (QQQ) GO TO 10
      IF (JOB.EQ.1) BETA = BOPT
      IF (JOB.EQ.1) IN = 0
      IF (IN.LE.INT-3) GO TO 10
C COMPUTE THREE SUCCESSIVE RATIOS OF THE DISPLACEMENT VECTOR NORMS
C AT THE END OF THE FIRST INTERVAL
      L = IN - (INT-3)
      GX(L) = SUM/S
      IF (IN.LT.INT) GO TO 10
      QQQ = .TRUE.
      D1 = GX(1) - GX(2)
      D2 = GX(2) - GX(3)
C CRITERION FOR DETERMINING THE FEASIBILITY OF AITKEN*S
C EXTRAPOLATION, THE DIFFERENCES BETWEEN THE RATIOS MUST DECREASE
C AND HAVE THE SAME SIG
      IF (ABS(D1).LE.ABS(D2) .OR. (D1/D2).LT.0.0) GO TO 20
      GAM = GX(1) - ((GX(2)-GX(1))**2/(GX(1)-2.*GX(2)+GX(3)))
      GO TO 40
   20 WRITE (6,99998)
   30 GAM = SUM/S
   40 IN = 0
      BL = BOPT
C COMPUTE A NEW ESTIMATE OF BETA OPTIMUM
      BOPT = 2./(1.+SQRT(ABS(1.-(GAM+BETA-1.)**2/(GAM*BETA**2))))
C THE FOLLOWING BOPT REDUCTION PREVENTS OVERESTIMATION OF BETA
      BETA = BOPT - (2.-BOPT)/4.
C CRITERION FOR STOPPING THE PROCESS OF IMPROVING BETA.
      IF (ABS(BOPT-BL)/(2.-BOPT).GT.0.05) GO TO 50
C BETA OPTIMUM HAS BEEN FOUND
      Q = .FALSE.
      IN = INT
      BETA = BOPT
      GAM = BETA - 1.
      GR = GAM/(2.-BETA)
      WRITE (6,99997) BETA, JOB
   50 IF (GAM.EQ.1.) GO TO 10
      ERM = BIG*GAM/(1.-GAM)
      GO TO 70
   60 ERM = BIG*GR
   70 IF (ABS(ERM).LE.ER) GO TO 80
      IF (JOB.LT.IP) GO TO 10
   80 WRITE (6,99996) JOB, ERM
C IF BETA OPTIMUM HAS NOT BEEN FOUND YET, BUT ERM.LE.ER, PRINT BETA.
      IF (Q) WRITE (6,99995) BETA
      RETURN
99999 FORMAT (//1X, 30HSOLUTION BY SOR (CARRE) METHOD)
99998 FORMAT (/5X, 45HAITKEN*S EXTRAPOLATION NOT FEASIBLE FOR BETA ,
     * 13HOPTIMIZATION.)
99997 FORMAT (5X, 14HBETA OPTIMUM =, E14.7, 6H AFTER, I4, 7H ITERAT,
     * 5HIONS.)
99996 FORMAT (/10X, 22HNUMBER OF ITERATIONS =, I4/10X, 9HMAXIMUM E,
     * 4HRROR, 8X, 1H=, E10.3, 8H PERCENT)
99995 FORMAT (10X, 22HBETA (NOT OPTIMUM)   =, E14.7)
      END
      SUBROUTINE SOR(M, N, MN, A, IJ, C, BETA, SUM, BIG)                SOR   10
C THIS IS THE SUCCESSIVE OVERRELAXATION ALGORITHM.
C THE SUBROUTINE PERFORMS ONE COMPLETE SWEEP OF THE INTERIOR, AND
C COMPUTES THE ERROR CRITERION COMPONENTS.
      DIMENSION A(M,N), IJ(MN,3), C(MN)
      COMMON /CS/ NR /LOG/ Q, QQ
      EXTERNAL AV
      LOGICAL Q, QQ
      DO 10 K=1,NR
        I = IJ(K,1)
        J = IJ(K,2)
        L = IJ(K,3)
        AL = A(I,J)
        AN = AV(M,N,MN,A,C,I,J,K,L)
        AN = AL + BETA*(AN-AL)
        A(I,J) = AN
        RES = ABS(AN-AL)
        IF (Q) SUM = SUM + RES
        IF (QQ) GO TO 10
        IF (RES.GT.BIG) BIG = RES
   10 CONTINUE
      RETURN
      END
      FUNCTION AV(M, N, MN, A, C, I, J, K, L)                           AV    10
      DIMENSION A(M,N), C(MN)
      GO TO (10, 20, 30, 40, 50), L
C NEUMANN B.C. ARE PARALLEL TO J-AXIS
   10 AV = (2.*A(I+1,J)+C(K)+C(K)+A(I,J+1)+A(I,J-1))/4.
      RETURN
C NEUMANN B.C. ARE PARALLEL TO I-AXIS
   20 AV = (2.*A(I,J+1)+C(K)+C(K)+A(I+1,J)+A(I-1,J))/4.
      RETURN
C NEUMANN B.C. AS IN 1, BUT FOR I.GT.IN 1
   30 AV = (2.*A(I-1,J)+C(K)+C(K)+A(I,J+1)+A(I,J-1))/4.
      RETURN
C NEUMANN B.C. AS IN 2, BUT FOR J.GT. IN 2
   40 AV = (2.*A(I,J-1)+C(K)+C(K)+A(I+1,J)+A(I-1,J))/4.
      RETURN
C ORDINARY 5-POINT OPERATOR FOR THE INTERIOR.
   50 AV = (A(I,J+1)+A(I+1,J)+A(I-1,J)+A(I,J-1))/4.
      RETURN
      END
      SUBROUTINE BOUND(M, N, MM, NN, A, NET)                            BOU   10
C THE SUBROUTINE READS THE BOUNDARY CONDITIONS SPECIFIED IN THE
C DATA CARDS AND SETS THEM INTO THE ARRAYS A AND NET. THE BOUNDARY
C CONDITIONS ARE DETERMINED BY SEGMENTS P1P2 WHERE POINTS P1(I1,J1)
C AND P2(I2,J2) COINCIDE WITH THE MESH POINTS IMPOSED ON THE REGION
C A(M,N). THE POINTS MUST SATISFY THE FOLLOWING CONDITION P1.LE.P2.
C THE VALUES AB ASSOCIATED WITH THE POINTS P1 AND P2 CAN BE AB1.LE.AB2
C OR AB1.GE.AB2. DIRICHLET B.C. ARE DISTINGUISHED FROM NEWMANN B.C.
C BY THE INDICATOR IND WHERE IND=1 FOR DIRICHLET OR IND=2 FOR NEUMANN
C B.C.  EACH DATA CARD CONTAINS ONLY ONE SEGMENT P1P2 (OR A POINT).
C THE INDEX NEXT.NE.0 TERMINATES THE READING PROCESS FOR ONE PROBLEM.
      DIMENSION A(M,N), NET(MM,NN)
      COMMON /COEF/ ML, NL /TEMP/ ID(5), IY(5)
C CLEAR THE ARRAYS.
      DO 20 J=1,N
        DO 10 I=1,M
          A(I,J) = 0.0
          NET(I+1,J+1) = 0
   10   CONTINUE
   20 CONTINUE
C SET BORDERS IN THE ARRAYS A AND NET.
      DO 30 I=2,ML
        NET(I,2) = 1
        NET(I,NL) = 1
   30 CONTINUE
      DO 40 J=2,NL
        NET(2,J) = 1
        NET(ML,J) = 1
   40 CONTINUE
      DO 50 I=1,MM
        NET(I,1) = 1
        NET(I,NN) = 1
   50 CONTINUE
      DO 60 J=1,NN
        NET(1,J) = 1
        NET(MM,J) = 1
   60 CONTINUE
   70 READ (5,99999) I1, J1, AB1, I2, J2, AB2, IND, NEXT
      IF (I1.GT.I2 .OR. J1.GT.J2) GO TO 150
C PLACE BOUNDARY PARALLEL TO I-AXIS INTO A.
      IF (I1.EQ.I2) GO TO 90
      AB = (AB2-AB1)/FLOAT(I2-I1)
      DO 80 I=I1,I2
        A(I,J1) = AB1 + AB*FLOAT(I-I1)
C THE FOLLOWING STATEMENT TRUNCATES ABOVE 5TH DECIMAL DIGIT.
        IF (AB.NE.0.0) A(I,J1) = AINT(A(I,J1)*1.E5)*1.E-5
        NET(I+1,J1+1) = IND
   80 CONTINUE
      GO TO 120
C PLACE BOUNDARY PARALLEL TO J-AXIS INTO A.
   90 IF (J1.EQ.J2) GO TO 110
      AB = (AB2-AB1)/FLOAT(J2-J1)
      DO 100 J=J1,J2
        A(I1,J) = AB1 + AB*FLOAT(J-J1)
        IF (AB.NE.0.0) A(I1,J) = AINT(A(I1,J)*1.E5)*1.E-5
        NET(I1+1,J+1) = IND
  100 CONTINUE
      GO TO 120
C PLACE THE POINT (I1,J1) INTO A.
  110 A(I1,J1) = AB1
      NET(I1+1,J1+1) = IND
  120 IF (NEXT.EQ.0) GO TO 70
C PRINT THE ARRAY NET.
      WRITE (6,99998)
      ND = 1 + (N-1)/10
      DO 130 J=1,ND
        ID(J) = J
  130 CONTINUE
      WRITE (6,99997) (ID(J),J=1,ND)
      DO 140 I=1,M
        WRITE (6,99996) I, (NET(I+1,J),J=2,NL)
  140 CONTINUE
      RETURN
  150 WRITE (6,99995)
      STOP
99999 FORMAT (2(2I5, F10.4), 2I5)
99998 FORMAT (1H1, 10X, 40HTHE FOLLOWING MATRIX REPRESENTS THE NET,,
     * 5X, 11H0=VARIABLES/56X, 26H1=DIRICHLET BOUNDARY COND./56X,
     * 24H2=NEUMANN BOUNDARY COND./)
99997 FORMAT (1H , 35X, 5(I1, 19X)/)
99996 FORMAT (11X, I2, 4X, 50I2)
99995 FORMAT (1H1, 1X, 41HINCORRECT INPUT DATA, I1.GT.I2 OR J1.GT.J,
     * 2H2.)
      END
      SUBROUTINE DINE(M, N, MM, NN, MN, A, NET, IJ, C)                  DIN   10
C THE SUBROUTINE RETURNS REFERENCE MARKS LOCATED IN THE ARRAY IJ(I,3).
C THE MARKS WILL ACCOMPANY THE SCAN AND HENCE SIMPLIFY THE LOGIC OF
C THE ITERATION PROCESS. BOUNDARY CONDITIONS IN THE ARRAY NET ARE
C IDENTIFIED BY 1= DIRICHLET, 2= NEUMANN, AND 0= INTERIOR OF THE
C REGION.
      DIMENSION A(M,N), NET(MM,NN), IJ(MN,3), C(MN)
      COMMON /CS/ NR
C SET MARKS OF REFERENCE, AND SAVE NEUMANN B.C. FROM A INTO C.
C 1-4 = NEUMANN B.C. LOCATED AT (I-1,J),(I,J-1),(I+1,J), AND (I,J+1)
C WITH RESPECT TO THE INTERIOR OF THE REGION, RESPECTIVELY.
C 5-INTERNAL (REGULAR) POINTS.  NOTE THAT THE MARK CHANGE FROM 0 TO 5
C IS INTERNAL TO THE PROGRAM AND IS INTRODUCED TO SIMPLIFY THE
C SELECTION OF EQUATIONS PERTINENT TO EITHER THE NEUMANN POINTS (1-4)
C OR THE REGULAR POINTS (5) /SEE FUNCTION AV/.
      DO 30 K=1,NR
        I = IJ(K,1) + 1
        J = IJ(K,2) + 1
        IF (NET(I,J)-1) 10, 40, 20
   10   IJ(K,3) = 5
        C(K) = 0.0
        GO TO 30
   20   IF (NET(I+1,J).EQ.0) IJ(K,3) = 1
        IF (NET(I,J+1).EQ.0) IJ(K,3) = 2
        IF (NET(I-1,J).EQ.0) IJ(K,3) = 3
        IF (NET(I,J-1).EQ.0) IJ(K,3) = 4
        C(K) = A(I-1,J-1)
   30 CONTINUE
      WRITE (6,99999) NR
      RETURN
   40 WRITE (6,99998) I, J
      STOP
99999 FORMAT (/10X, 30HNUMBER OF SCANNED POINTS, NR= , I5)
99998 FORMAT (/1X, 36HINCORRECT SCANNING SEQUENCE AT NET (, I2, 1H,,
     * I2, 2H).)
      END
      SUBROUTINE INVA(M, N, MM, NN, A, NET)                             INV   10
C THE SUBROUTINE READS A CONSTANT INITIAL VALUE, AND INTRODUCES IT
C INTO THE ARRAY A. THE INITIAL VALUE DATUM FOLLOWS THE BOUNDARY
C CONDITION DATA.
      DIMENSION A(M,N), NET(MM,NN)
      READ (5,99999) V
      DO 20 J=1,N
        DO 10 I=1,M
          IF (NET(I+1,J+1).NE.1) A(I,J) = V
   10   CONTINUE
   20 CONTINUE
      WRITE (6,99998) V
      RETURN
99999 FORMAT (F10.0)
99998 FORMAT (10X, 29HCONSTANT INITIAL VALUE,    V=, F7.2)
      END
      SUBROUTINE PRES(M, N, A)                                          PRE   10
C THE SUBROUTINE PRINTS NUMERICAL RESULTS FROM THE ARRAY A.
C EACH PRINTED PAGE WILL CONTAIN 10 COLUMNS OF A.
      DIMENSION A(M,N)
      COMMON /TEMP/ LA(10)
      DO 30 L=1,N,10
        JS = L - 1
        JE = L + 9
        DO 10 I=1,10
          LA(I) = JS + I
   10   CONTINUE
        WRITE (6,99999) LA
        DO 20 I=1,M
          WRITE (6,99998) I, (A(I,J),J=L,JE)
   20   CONTINUE
   30 CONTINUE
      RETURN
99999 FORMAT (1H1, 15X, 10(I3, 9X)/)
99998 FORMAT (1H , 5H   A(, I2, 1H,, 3X, 10(F10.5, 2X))
      END
      SUBROUTINE MAP(M, N, MM, NN, A, NET)                              MAP   10
C THE SUBROUTINE PRINTS THE EQUIPOTENTIAL MAP.
      DIMENSION A(M,N), NET(MM,NN), CA(10), SP(10)
      COMMON /TEMP/ IZ(10)
      DATA CA(1), CA(2), CA(3), CA(4), CA(5), CA(6), CA(7), CA(8),
     * CA(9), CA(10) /1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9/,
     * ZERO, VAL, BLANK /1HZ,1H$,1H /
C READ EQUIPOTENTIALS TO BE PRINTED OUT.
      READ (5,99999) (SP(I),I=1,10), TOL
C PREPARE THE MAP.
      DO 50 J=1,N
        DO 40 I=1,M
          IF (NET(I+1,J+1).NE.1) GO TO 20
          IF (A(I,J).EQ.0.0) GO TO 10
          A(I,J) = VAL
          GO TO 40
   10     A(I,J) = ZERO
          GO TO 40
   20     DO 30 K=1,10
            IF (A(I,J).LT.SP(K)-TOL .OR. A(I,J).GT.SP(K)+TOL) GO TO
     *       30
            A(I,J) = CA(K)
            GO TO 40
   30     CONTINUE
          A(I,J) = BLANK
   40   CONTINUE
   50 CONTINUE
C WRITE THE MAP.
      WRITE (6,99998)
      DO 60 I=1,M
        WRITE (6,99997) (A(I,J),J=1,N)
   60 CONTINUE
C WRITE THE LIMITS FOR THE EQUIPOTENTIAL STRIPS IN THE PRINTED MAP.
      WRITE (6,99996)
      DO 70 I=1,10
        IZ(I) = I - 1
   70 CONTINUE
      DO 80 I=1,5
        P1 = SP(I) - TOL
        P2 = SP(I) + TOL
        P3 = SP(I+5) - TOL
        P4 = SP(I+5) + TOL
        IF (P1.LT.0.) P1 = 0.
        WRITE (6,99995) IZ(I), P1, P2, IZ(I+5), P3, P4
   80 CONTINUE
      WRITE (6,99994)
      RETURN
C IF THE ARRAY A IS LARGER THAN 30X30, CHANGE THE FORMAT 99997
C TO (15X,50(A1,1X)).
99999 FORMAT (10F7.3, F10.3)
99998 FORMAT (1H1, 7X, 41HTHE FOLLOWING IS A CONTOUR MAP OF THE EQU,
     * 31HIPOTENTIALS WITHIN GIVEN LIMITS/)
99997 FORMAT (15X, 20(A1, 2X)/)
99996 FORMAT (/13X, 8HCONTOURS)
99995 FORMAT (13X, I1, 1H=, E11.4, 4H TO , E11.4, 6H  **  , I1, 1H=,
     * E11.4, 4H TO , E11.4)
99994 FORMAT (13X, 10HBOUNDARIES/13X, 22HZ= ZERO DIRICHLET B.C.,
     * 8X, 4H**  , 25H$= NONZERO DIRICHLET B.C.)
      END

.