      SUBROUTINE GOURAU(ICTFAC,NDS,VVV,DIRLUM,PHONG,NF,BSOMB,IOPT)
      DIMENSION VVV(4),DIRLUM(3),PHONG(3,*)
      INTEGER   NF(*)
C
      IF (ICTFAC.EQ.97) THEN
        IF (IOPT.EQ.0) THEN
          DO K=1,NDS
            VVV(K) = ABS( DIRLUM(1)*PHONG(1,NF(K))
     &                  + DIRLUM(2)*PHONG(2,NF(K))
     &                  + DIRLUM(3)*PHONG(3,NF(K)))
          ENDDO
        ELSEIF (NF(5).EQ.0) THEN
          DO K=1,NDS
            VVV(K) = ABS( DIRLUM(1)*PHONG(1,NF(K))
     &                  + DIRLUM(2)*PHONG(2,NF(K))
     &                  + DIRLUM(3)*PHONG(3,NF(K)))
          ENDDO
        ELSE
          DO K=1,NDS
            VVV(K) = BSOMB+.5*ABS( DIRLUM(1)*PHONG(1,NF(K))
     &                           + DIRLUM(2)*PHONG(2,NF(K))
     &                           + DIRLUM(3)*PHONG(3,NF(K)))
          ENDDO
        ENDIF
      ELSE
        IF (IOPT.EQ.0) THEN
          DO K=1,NDS
            TOTO = ABS( DIRLUM(1)*PHONG(1,NF(K))
     &                + DIRLUM(2)*PHONG(2,NF(K))
     &                + DIRLUM(3)*PHONG(3,NF(K)))
            VVV(K) = 2.5*TOTO*TOTO*(TOTO-1.)+TOTO
Cfj            VVV(K) = .1+.9*( DIRLUM(1)*PHONG(1,NF(K))
Cfj     &                     + DIRLUM(2)*PHONG(2,NF(K))
Cfj     &                     + DIRLUM(3)*PHONG(3,NF(K)))**4
          ENDDO
        ELSEIF (NF(5).EQ.0) THEN
          DO K=1,NDS
            TOTO = ABS( DIRLUM(1)*PHONG(1,NF(K))
     &                + DIRLUM(2)*PHONG(2,NF(K))
     &                + DIRLUM(3)*PHONG(3,NF(K)))
            VVV(K) = 2.5*TOTO*TOTO*(TOTO-1.)+TOTO
Cfj            VVV(K) = .1+.9*( DIRLUM(1)*PHONG(1,NF(K))
Cfj     &                     + DIRLUM(2)*PHONG(2,NF(K))
Cfj     &                     + DIRLUM(3)*PHONG(3,NF(K)))**4
          ENDDO
        ELSE
          DO K=1,NDS
            TOTO = ABS( DIRLUM(1)*PHONG(1,NF(K))
     &                + DIRLUM(2)*PHONG(2,NF(K))
     &                + DIRLUM(3)*PHONG(3,NF(K)))
            VVV(K) = BSOMB+.5*(2.5*TOTO*TOTO*(TOTO-1.)+TOTO)
Cfj            VVV(K) = BSOMB+.5*( DIRLUM(1)*PHONG(1,NF(K))
Cfj     &                        + DIRLUM(2)*PHONG(2,NF(K))
Cfj     &                        + DIRLUM(3)*PHONG(3,NF(K)))**4
          ENDDO
        ENDIF
      ENDIF
      IF (NDS.EQ.3) VVV(4) = VVV(3)
      END
