      SUBROUTINE  BFLUX (FLUX, COORD, LBN, N, NSPACE, NFLUX, 
     1            NG, C, S, IOPT, NQP, NPARM, H, DGH, 
     2            PT, WT, XYZ, DLH, G, DLG, AJ, AJINV )
C     * * * * * * * * * * * * * * * * * * * * * * * * * *
C         PROBLEM DEPENDENT BOUNDARY FLUX CONTRIBUTIONS
C     * * * * * * * * * * * * * * * * * * * * * * * * * *
C        ALWAYS USED
      DIMENSION  COORD(LBN,NSPACE), FLUX(LBN,NG), C(NFLUX),
     1           S(NFLUX,NFLUX)
C        OPTIONAL FOR NUMERICAL INTEGRATION
      DIMENSION   H(N), DGH(NSPACE,N), PT(NPARM,NQP), 
     1            WT(NQP), XYZ(3), DLH(NSPACE,N), G(LBN), 
     2            DLG(NPARM,LBN), AJ(NSPACE,NSPACE), 
     3            AJINV(NSPACE,NSPACE)
C
C     FLUX    = SPECIFIED BOUNDARY FLUX COMPONENTS
C     COORD   = SPATIAL COORDINATES OF SEGMENT NODES
C     LBN     = NO. OF NODES ON AN ELEMENT BOUNDARY SEGMENT
C     NSPACE  = DIMENSION OF SOLUTION SPACE
C     NFLUX   = N*NG = MAXIMUM NUMBER OF FLUX CONTRIBUTIONS
C     C       = BOUNDARY FLUX COLUMN MATRIX CONTRIBUTIONS
C     S       = BOUNDARY FLUX SQUARE MATRIX
C     NG      = NUMBER OF PARAMETERS PER NODE POINT
C     IOPT    = PROBLEM MATRIX REQUIREMENTS, MUST BE SET.
C             = 1, CALCULATE C ONLY
C             = 2, CALCULATE S ONLY
C             = 3, CALCULATE BOTH C AND S
C     AJ      = JACOBIAN
C     AJINV   = JACOBIAN INVERSE
C     DGH     = GLOBAL DERIVATIVES INTERPOLATION FUNCTIONS
C     DLG     = LOCAL DERIVATIVES GEOMETRIC INTERPOLATION
C     DLH     = LOCAL DERIVATIVES INTERPOLATION FUNCTIONS
C     G       = GEOMETRIC INTERPOLATION FUNCTIONS
C     H       = SOLUTION INTERPOLATION FUNCTIONS
C     N       = NUMBER OF SOLUTION NODES, N = LBN USUALLY
C     NQP     = NUMBER OF QUADRATURE POINTS
C     NPARM   = PARAMETRIC GEOMETRY NODES, = NSPACE USUALLY
C     PT      = QUADRATURE COORDINATES
C     WT      = QUADRATURE WEIGHTS
C     XYZ     = SPACE COORDINATES AT A POINT
      IOPT = 0
C     ....................................................
C       ** BFLUX PROBLEM DEPENDENT STATEMENTS FOLLOW **
C     ....................................................
C     LINEAR VARIATION OF 2 COMPONENTS ON STRAIGHT LINE
CDP   SQRT(Z) = DSQRT(Z)
      IOPT = 1
      DX = COORD(2,1) - COORD(1,1)
      DY = COORD(2,2) - COORD(1,2)
      DL = SQRT (DX*DX + DY*DY)
      DLD6 = DL/6.
C     ODD = X FORCE, EVEN = Y FORCE
      C(1) = DLD6*( 2.*FLUX(1,1) + FLUX(2,1) )
      C(2) = DLD6*( 2.*FLUX(1,2) + FLUX(2,2) )
      C(3) = DLD6*( FLUX(1,1) + 2.*FLUX(2,1) )
      C(4) = DLD6*( FLUX(1,2) + 2.*FLUX(2,2) )
      RETURN
      END
