C  Solution of Laplace's equation in 2 dimensions
C  using the 4-color SOR iteration, 9-point difference op.
C  Parameters R, C give array dimensions

C  Four-color scheme is   1 2 1 2 1 2 . . .
C                         4 3 4 3 4 3 . . .
C                         1 2 1 2 1 2 . . .
C                         4 3 4 3 4 3 . . .
C                         . . . . . . . . .

C  fourc9.h  defines parameters and common arrays for
C  managing it all:
C
C     R, C  - arrays are (R, C)
C     I, J  - Coordinate arrays
C      Clr  - Color (1, 2, 3, or 4)
C
C  Parallel arrays are in /FOURC9/

C  Perform n cycles of 4-color SOR on U with parameter W
C  Returns updated function, writes error estimate in mxr
C  Error estimate is max change over n iterations.

      function Sor4(U, n, W, mxr)
      INCLUDE 'fourc9.h'
      real, array(R, C) :: Sor4, U
      integer n
      real W, mxr

      real, array(R, C) ::  W2, W1, Wr
      integer k, kc
      real onemW

C   Wr = * + *

      onemW = 1.0 - W
D     print *, " W, onemW"
D     print *, W, onemW
      W1 = U
      W2 = U
      Wr = U
      do, k = 1, n
       do  kc = 1, 4		! "Colors"
        Wr(:, 2:C-1) = W1(:, 1:C-2) + W1(:, 3:C)
D       print *, " Wr -", kc
D       call pr4(Wr)
        where ( Clr(2:R-1, 2:C-1) == kc )
	 W2(2:R-1, 2:C-1) =
     *	     0.20*(  Wr(2:R-1 ,2:C-1)
     *             + W1(1:R-2, 2:C-1) + W1(3:R, 2:C-1))
     *	   + 0.05*(  Wr(1:R-2, 2:C-1) + Wr(3:R, 2:C-1))
         W1(2:R-1, 2:C-1) =
     *        W*W2(2:R-1, 2:C-1) + onemW*W1(2:R-1, 2:C-1)
        end where
D       print *, " W2 -", kc
D       call pr4(W2)
D       print *, " W1 -", kc+1
D       call pr4(W1)
       end do			! "Colors"
      end do
      mxr = maxval(abs(W1 - U))	! error estimate
      Sor4 = W1
      end function Sor4

C  Print 4 by 4 array section

      subroutine pr4(W)
      INCLUDE 'fourc9.h'
      real, array(R, C) :: W
CMPF  ONDPU W
      print *, W(1, 1:4)
      print *, W(2, 1:4)
      print *, W(3, 1:4)
      print *, W(4, 1:4)
      end subroutine pr4

C  Calculate four-color control array Clr

      subroutine Fourc()
      INCLUDE 'fourc9.h'
      where ( IAND(J, 1) == 1 )
       Clr = 1
      elsewhere
       Clr = 2
      end where
      where ( IAND(I, 1) == 0 )
       Clr = 5 - Clr
      end where
D     print *, " Clr"
D     print *, Clr(1, 1:4)
D     print *, Clr(2, 1:4)
D     print *, Clr(3, 1:4)
D     print *, Clr(4, 1:4)
      end subroutine Fourc


C  Test 4-color SOR

      program TestFourc

      INCLUDE 'fourc9.h'

      INTERFACE
	function Sor4(U, n, W, mxr)
	INCLUDE 'fourc9.h'
	real, array(R, C) :: Sor4, U
CMPF	ONDPU U
	integer n
	real W, mxr
      END INTERFACE

      real, array(R, C) :: U0, Z
      real bsum, W, eps, err, tol
      integer it, in

      I = SPREAD( [1:R], 2, C )   ! Coordinates
      J = SPREAD( [1:C], 1, R )
      Call Fourc()                ! Clr = Colors (1,2,3,4)

C     Z is a solution of the difference equation

      Z = -1.4 + 0.12*I - 0.23*J + 0.001*(I**2 - J**2)

      U0 = Z	! Install boundary values

      bsum = SUM(U0(1,:)) + SUM(U0(R,:)) ! Sum on boundary
     *     + SUM(U0(:,1)) + SUM(U0(:,C))

C     Intialize interior to average on boundary
C     The "-4" prevents counting corners twice

      U0(2:R-1, 2:C-1) = bsum/(2*R + 2*C - 4)

C     Get W
      print *, " W, tol, in >"
      read *, W, tol, in
      eps = 2.0 + tol + tol	! Forces iterations
      print *, " it, eps, err"
      do while ( eps > tol )
D       print *, " U0"
D       call pr4(U0)
	U0 = Sor4(U0, in, W, eps)
        err = maxval(abs(U0 -  Z))
        it = it + in
        print *, it, eps, err
      end do

      end program TestFourc
