      SUBROUTINE RELAXA(U,F,NX,NY,NZ,H)
C
C     3D RED BLACK RELAXATION    STANDARD IMLEMENTATION
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION U(NX,NY,NZ),F(NX,NY,NZ)
      H2 = H*H
      FACTOR = 1.0D0/6.0D0
C
C     RELAXATION OF THE RED POINTS
C
      IA = 0
      DO 10 K = 2, NZ-1
        DO 20 J = 2, NY-1
          DO 30 I = 2 + IA, NX-1, 2
            U(I,J,K) = FACTOR*(H2*F(I,J,K) +
     *                    (U(I-1,J,K) + U(I+1,J,K) +
     *                     U(I,J-1,K) + U(I,J+1,K) +
     *                     U(I,J,K-1) + U(I,J,K+1)))
30        CONTINUE
        IA = MOD (IA+1,2)
20      CONTINUE
10    CONTINUE
C
C     RELAXATION OF THE BLACK POINTS
C
      IA = 1
      DO 11 K = 2, NZ-1
        DO 21 J = 2, NY-1
          DO 31 I = 2 + IA, NX-1, 2
            U(I,J,K) = FACTOR*(H2*F(I,J,K) +
     *                    (U(I-1,J,K) + U(I+1,J,K) +
     *                     U(I,J-1,K) + U(I,J+1,K) +
     *                     U(I,J,K-1) + U(I,J,K+1)))
31        CONTINUE
        IA = MOD (IA+1,2)
21      CONTINUE
11    CONTINUE
      END
