
      SUBROUTINE INIT1 (NP, H, F, G, U, FR)
      IMPLICIT NONE
      INTEGER NP
      DOUBLE PRECISION F, G, H, FR (NP, NP), U (NP, NP)
      EXTERNAL F, G
C
C     COMPUTES INITIAL VALUES ON THE FINEST GRID
C
        INTEGER I, J, N
        DOUBLE PRECISION H2, X, Y

        H2 = H * H
        N  = NP - 1
        DO J = 2, N
          Y = (J - 1) * H
          U (1, J) = G (0.0D0, Y)
          DO I = 2, N
            X = (I - 1) * H
            U  (I, J) = 0.0D0
            FR (I, J) = H2 * F (X, Y)
          ENDDO
          U (NP, J) = G (1.0D0, Y)
        ENDDO
        DO I = 1, NP
          X = (I - 1) * H
          U (I, 1)  = G (X, 0.0D0)
          U (I, NP) = G (X, 1.0D0)
        ENDDO
      END
C
C.......................................................................
C
      SUBROUTINE INT2A (NPC, NPF, UC, UF)
      IMPLICIT NONE
      INTEGER NPC, NPF
      DOUBLE PRECISION UC (NPC, NPC), UF (NPF, NPF)
C
C     COARSE-TO-FINE TRANSFER (BILINEAR INTERPOLATION) AND CORRECTION
C
        INTEGER NC, I, IF, J, JF

        NC = NPC - 1
        DO J= 2, NC
          JF = 2 * J - 1
          DO I = 1, NC
            IF = 2 * I
            UF (IF, JF) = UF (IF, JF) +
     1                    0.5D0 * (UC (I, J) + UC (I + 1, J))
          ENDDO
        ENDDO
        DO J = 1, NC
          JF = 2 * J
          DO I = 1, NC
            IF = 2 * I - 1
            UF (IF, JF) = UF (IF, JF) +
     1                    0.5D0 * (UC (I, J) + UC (I, J + 1))
          ENDDO
        ENDDO
      END
C
C.......................................................................
C
      SUBROUTINE INT4 (NPC, NPF, UC, UF, FRF)
      IMPLICIT NONE
      INTEGER NPC, NPF
      DOUBLE PRECISION FRF (NPF, NPF), UC (NPC, NPC), UF (NPF, NPF)
C
C     4-TH ORDER FMG-INTERPOLATION
C
        INTEGER NC, I, IF, J, JF, JM
        DOUBLE PRECISION HA, HN
C
C       TRANSFER OF COARSE GRID VALUES TO THE FINE GRID
C
        NC = NPC - 1
        JF = 3
        DO J = 2, NC
          IF = 3
          DO I = 2, NC
            UF (IF, JF) = UC (I, J)
            IF = IF + 2
          ENDDO
          JF = JF + 2
        ENDDO
C
C       COMPUTATION OF THE REMAINING EVEN POINTS BY USING THE
C       ROTATED 5-POINT STAR
C
        JF = 2
        JM = 1
        DO J = 2, NPC
          IF = 2
          HA = UC (1, JM) + UC (1, J)
          DO I = 2, NPC
            HN = UC (I, JM) + UC (I, J)
            UF (IF, JF) = 0.25D0 * (FRF (IF, JF) + FRF (IF, JF) +
     1                              HA + HN)
            HA = HN
            IF = IF + 2
          ENDDO
          JM = JM + 1
          JF = JF + 2
        ENDDO
C
C       COMPUTATION OF THE ODD POINTS BY ONE HALF (ODD) RELAXATION STEP
C       **** WHAT?? -VE ITMAX DOESN'T DO THAT !! ****
C
        CALL RELAX (-1, .FALSE., NPF, UF, FRF)
      END
C
C.......................................................................
C
      SUBROUTINE PUTZB (NPC, UC)
      IMPLICIT NONE
      INTEGER NPC
      DOUBLE PRECISION UC (NPC, NPC)
C
C     PUTS ZERO TO BOUNDARY OF COARSER GRIDS
C
        INTEGER NC, I, J

        NC = NPC - 1
        DO J = 2, NC
          UC (1, J)   = 0.0D0
          UC (NPC, J) = 0.0D0
        ENDDO
        DO I =1, NPC
          UC (I, 1)   = 0.0D0
          UC (I, NPC) = 0.0D0
        ENDDO
      END
C
C.......................................................................
C
      SUBROUTINE RELAX (ITMAX, IZ, NPF, UF, FF)
      IMPLICIT NONE
      INTEGER ITMAX, NPF
      LOGICAL IZ
      DOUBLE PRECISION FF (NPF, NPF), UF (NPF, NPF)
C
C     RED BLACK RELAXATION ON GRID NPF
C
        INTEGER I, IS, IT, J, NXF, NYF, NXPF, NYPF, NXF2, NYF2

        NXPF = NPF
        NYPF = NPF
        NXF  = NXPF - 1
        NYF  = NYPF - 1
        NXF2 = NXF - 1
        NYF2 = NYF - 1
        DO IT = 1, ITMAX
          IS = 2
          IF (IZ .AND. (IT .LE. 1)) THEN
C
C           ZERO STARTING VALUES
C           RELAXATION OF RED POINTS
C
            DO J = 2, NYF, 2
              DO I = 2, NXF, 2
                UF (I, J) = FF (I, J) * 0.25
              ENDDO
            ENDDO
            DO J = 3, NYF, 2
              DO I = 3, NXF, 2
                UF (I, J) = FF (I, J) * 0.25
              ENDDO
            ENDDO

          ELSE
C
C           NON-ZERO STARTING VALUES
C           RELAXATION OF RED POINTS
C
            DO J = 2, NYF2, 2
              DO I = IS, NXF2, 2
                UF (I, J) = 0.25 * (FF (I, J) +
     1                   (UF (I - 1, J) + (UF (I, J - 1) +
     2                   (UF (I + 1, J) +  UF (I, J + 1)))))
                UF (I + 1, J + 1) = 0.25 * (FF (I + 1, J + 1) +
     1                   (UF (I + 2, J + 1) + (UF (I + 1, J + 2) +
     2                   (UF (I + 1, J)     +  UF (I, J + 1)))))
              ENDDO
            ENDDO
            DO I = 2, NXF, 2
              UF (I, NYF) = 0.25 * (FF (I, NYF) +
     1                 (UF (I - 1, NYF) + (UF (I + 1, NYF) +
     2                 (UF (I, NYF - 1) +  UF (I, NYF + 1)))))
            ENDDO
            DO J = 2, NYF, 2
              UF (NXF, J) = 0.25 * (FF (NXF, J) + 
     1                 (UF (NXF, J - 1) + (UF (NXF, J + 1) +
     2                 (UF (NXF - 1, J) +  UF (NXF + 1, J)))))
            ENDDO
C
          ENDIF
          IS = 3
C
C         RELAXATION OF BLACK POINTS
C
          DO J = 2, NYF2, 2
            DO I = IS, NXF2, 2
              UF (I, J) = 0.25 * (FF (I, J) + 
     1                 (UF (I, J - 1) + (UF (I + 1, J) +
     2                 (UF (I - 1, J) +  UF (I, J + 1)))))
              UF (I - 1, J + 1) = 0.25 * (FF (I - 1, J + 1) +
     1                 (UF (I - 2, J + 1) + (UF (I - 1, J + 2) +
     2                 (UF (I - 1, J)     +  UF (I, J + 1)))))
            ENDDO
          ENDDO
          DO I = 3, NXF, 2
              UF (I, NYF) = 0.25 * (FF (I, NYF) +
     1                 (UF (I - 1, NYF) + (UF (I + 1, NYF) +
     2                 (UF (I, NYF - 1) +  UF (I, NYF + 1)))))
          ENDDO
          DO J = 3, NYF, 2
              UF (NXF, J) = 0.25 * (FF (NXF, J) +
     1                 (UF (NXF, J - 1) + (UF (NXF, J + 1) +
     2                 (UF (NXF - 1, J) +  UF (NXF + 1, J)))))
          ENDDO
        ENDDO
C
      END
C
C.......................................................................
C
      SUBROUTINE RESTR (NPC, NPF, FC, UF, FF)
      IMPLICIT NONE
      INTEGER NPC, NPF
      DOUBLE PRECISION FC (NPC, NPC), FF (NPF, NPF), UF (NPF, NPF)
C
C     COMPUTATION OF THE DEFECT AND FINE-TO-COARSE TRANSFER
C     (HALF-INJECTION)
C
        INTEGER NC, I, IF, J, JF

        NC = NPC - 1
        DO J = 2, NC
          JF = 2 * J - 1
          DO I = 2, NC
            IF = 2 * I - 1
            FC (I, J) = 2 * (FF (IF, JF) -
     1          4.0D0 * UF (IF, JF) + UF (IF, JF - 1) + UF (IF - 1, JF)
     2                          + UF (IF + 1, JF) + UF (IF, JF + 1))
          ENDDO
        ENDDO
      END
C
C.......................................................................
C
      SUBROUTINE TRANS (NPC, NPF, UC, FC, UF, FF)
      IMPLICIT NONE
      INTEGER NPC, NPF
      DOUBLE PRECISION FC (NPC, NPC), FF (NPF, NPF),
     1                 UC (NPC, NPC), UF (NPF, NPF)
C
C     TRANSFER OF F AND U FROM GRID NPF TO NPC
C
        INTEGER NC, I, IF, J, JF

        NC = NPC - 1
        DO J = 2, NC
          JF = 2 * J - 1
          UC (1, J) = UF (1, JF)
          DO I = 2, NC
            IF = 2 * I - 1
            UC (I, J) = UF (IF, JF)
            FC (I, J) = 4.0D0 * FF (IF, JF)
          ENDDO
          UC (NPC, J) = UF (NPF, JF)
        ENDDO
        DO I = 1, NPC
          IF = 2 * I - 1
          UC (I, 1)   = UF (IF, 1)
          UC (I, NPC) = UF (IF, NPF)
        ENDDO
      END
C
C.......................................................................
C
      DOUBLE PRECISION FUNCTION DEFMX (UC, FC, NP)
      IMPLICIT NONE
      INTEGER NP
      DOUBLE PRECISION UC (NP, NP), FC (NP, NP)
C
C     COMPUTES THE MAXIMUM NORM OF THE DEFECT
C
        INTEGER N, I, J
        DOUBLE PRECISION D, HSQR

        N = NP - 1
        HSQR  = DBLE (FLOAT (N * N))
        DEFMX = 0.0D0
        DO J = 2, N
          DO I = 2, N
            D = (FC (I, J) - 4.0D0 * UC (I, J) +
     1           UC (I - 1, J) + UC (I, J - 1) +
     2           UC (I, J + 1) + UC (I + 1, J)) * HSQR
            DEFMX = DMAX1 (DEFMX, DABS (D))
          ENDDO
        ENDDO
      END
C
C.......................................................................
C
      DOUBLE PRECISION FUNCTION DIFMX (UC, NP, SOL)
      IMPLICIT NONE
      INTEGER NP
      DOUBLE PRECISION UC (NP, NP), SOL
      EXTERNAL SOL
C
C     COMPUTES THE MAXIMUM NORM OF THE DIFFERENCE BETWEEN
C     SOL (=SOLUTION OF THE BVP) AND THE VALUES IN UC
C
        INTEGER I, J
        DOUBLE PRECISION H, X, Y

        H     = 1.0D0 / (NP - 1)
        DIFMX = 0.0D0
        DO J = 1, NP
          Y = (J - 1) * H
          DO I = 1, NP
            X = (I - 1) * H
            DIFMX = DMAX1 (DIFMX, DABS (SOL (X, Y) - UC (I, J)))
          ENDDO
        ENDDO
      END

