      SUBROUTINE  FFLUX (NF,NBFREE,N,C,S,COORD,NSPACE,JOB,FLUX,
     1                  IGEN,H,DLH,NQ)
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C       FLUX RESULTANTS ON PLANAR FACES (OR ELEMENTS)
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
      COMMON /ELARG/ GPT(9),GWT(9),PT(2,12),WT(9),AJ(9),AJINV(9),
     1               DUM(24),DGH(2,12),DUM3(36),DUM4(144),PRTMAT(14),
     2    TWOPI,DUMMY,INJOB,NOBODY,LOADT,IOPT,MATNO,NG,NIP,ISOMAT
      DIMENSION  COORD(N,NSPACE),C(NBFREE),FLUX(NF),H(N,NQ),
     1           DLH(NSPACE,N,NQ),S(NBFREE,NBFREE)
      DATA KALL, NBSQ  / 1, 0 /
      IF ( KALL .EQ. 1 ) THEN
        KALL = 0
        IF ( JOB.NE.4 .AND. JOB.NE.6 )  STOP 'INVALID JOB IN FFLUX'
        NIP = NQ
        IF ( NQ.EQ.4 ) CALL GAUS2D (2,GPT,GWT,NIP,PT,WT)
        IF ( NQ.EQ.3 ) CALL SYMRUL (NQ,PT,WT)
        NBSQ = NBFREE*NBFREE
        DO 5 IQ = 1,NIP
          IF ( N.EQ.3 ) CALL SHP3T (PT(1,IQ),PT(2,IQ),H(1,IQ))
          IF ( N.EQ.4 ) CALL SHP4Q (PT(1,IQ),PT(2,IQ),H(1,IQ))
          IF ( N.EQ.6 ) CALL SHP6T (PT(1,IQ),PT(2,IQ),H(1,IQ))
          IF ( N.EQ.8 ) CALL SHP8Q (PT(1,IQ),PT(2,IQ),H(1,IQ))
          IF ( N.EQ.3 ) CALL DER3T (PT(1,IQ),PT(2,IQ),DLH(1,1,IQ))
          IF ( N.EQ.4 ) CALL DER4Q (PT(1,IQ),PT(2,IQ),DLH(1,1,IQ))
          IF ( N.EQ.6 ) CALL DER6T (PT(1,IQ),PT(2,IQ),DLH(1,1,IQ))
          IF ( N.EQ.8 ) CALL DER8Q (PT(1,IQ),PT(2,IQ),DLH(1,1,IQ))
  5     CONTINUE
      ENDIF    
      F1  = FLUX(1)
      F2  = FLUX(2)
      F12 = F1*F2
      IF ( NF .EQ. 3 ) F3 = FLUX(3)
C      ZERO WORKSPACES
      CALL ZEROA (NBFREE,C)
      IF ( JOB.EQ.4 ) CALL ZEROA (NBSQ,S)
C      COMPUTE MATRICES BY NUMERICAL INTEGRATION
      DO 900 IQ = 1,NIP
        CALL JACOB (N,NSPACE,DLH(1,1,IQ),COORD,AJ)
        CALL I2BY2 (AJ,AJINV,DET)
        CALL GDERIV (NSPACE,N,AJINV,DLH(1,1,IQ),DGH)
        DETWT = ABS(DET)*WT(IQ)
        IF ( JOB .EQ. 6 ) THEN
C          HYBRID PLATE BENDING PRESSURE RESULTANTS
          IGEN = 1
          X    = DOT(N,H(1,IQ),COORD)
          Y    = DOT(N,H(1,IQ),COORD(1,2))
          PRES = F1 + F2*X + F3*Y
          DO 26 I = 1,N
            I4 = 4*I
  26      C(I4) = C(I4) + H(I,IQ)*PRES*DETWT
        ENDIF
        IF ( JOB .EQ. 4 ) THEN
C          PLANAR FACE (THUS ELEMENT) CONVECTION
          IGEN = 3
          DO 24 I = 1,N
            C(I) = C(I) + H(I,IQ)*F12*DETWT
            DO 34 J = 1,N
  34        S(I,J) = S(I,J) + H(I,IQ)*H(J,IQ)*DETWT*F2
  24      CONTINUE
        ENDIF
  900 CONTINUE
      RETURN
      END
