      SUBROUTINE  ELSQ (N, NSPACE, NELFRE, NRB, NQP, NGEOM,
     1    NPARM, NNPFIX, NNPFLO, MISCFX, MISCFL, NLPFIX,
     2    NLPFLO, COORD, S, C, H, DGH, B, E, EB, STRAIN,
     3    STRAN0, STRESS, BODY, PT, WT, XYZ, DLH, G, DLG,
     4    AJ, AJINV, HINTG, D, PRTLPT, FLTMIS, ELPROP, PRTLPT, 
     5    PRTMAT, MISFIX, LSHAPE, LPROP, LPPROP, NTAPE1, NTAPE2, 
     6    NTAPE3, NTAPE4, NTAPE5, LNODE, NG, IE )
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     GENERATE ELEMENT SQUARE MATRIX, OPTIONAL COLUMN MATRIX
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * *
CDP   IMPLICIT   REAL*8  (A-H,O-Z)
C        ALWAYS USED
      DIMENSION  COORD(N,NSPACE), S(NELFRE,NELFRE)
C        USUALLY USED
      DIMENSION  C(NELFRE), H(N), DGH(NSPACE,N), B(NRB,NELFRE), 
     1           E(NRB,NRB), EB(NRB,NELFRE), STRAIN(NRB+2), 
     2           STRAN0(NRB), STRESS(NRB+2), BODY(NSPACE)
C        OPTIONAL FOR NUMERICAL INTEGRATION
      DIMENSION  PT(NPARM,0:NQP), WT(0:NQP), DLH(NSPACE,N), 
     1           G(NGEOM), DLG(NPARM,NGEOM), AJ(NSPACE,NSPACE), 
     2           AJINV(NSPACE,NSPACE), HINTG(N), LNODE(N),
     3           XYZ(NSPACE)
C        OPTIONAL PROPERTY AND SOLUTION VALUES
      DIMENSION  D(NELFRE), PRTLPT(N,0:NNPFLO), FLTMIS(0:MISCFL), 
     1           ELPROP(0:NLPFLO), PRTLPT(N,0:NNPFLO), 
     2           PRTMAT(0:NLPFLO), MISFIX(0:MISCFX), 
     3                             LPROP(0:NLPFIX), 
     4           LPPROP(0:NNPFIX)
C                     VARIABLES:
C     AJ      = JACOBIAN
C     AJINV   = JACOBIAN INVERSE
C     B       = STRAIN-DISPLACEMENT (GRADIENT) MATRIX
C     BODY    = BODY FORCE VECTOR
C     COORD   = SPATIAL COORDINATES OF ELEMENT'S NODES
C     D       = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
C     DGH     = GLOBAL DERIVATIVES INTERPOLATION FUNCTIONS
C     DLG     = LOCAL DERIVATIVES GEOMETRIC INTERPOLATION
C     DLH     = LOCAL DERIVATIVES INTERPOLATION FUNCTIONS
C     E       = CONSTITUTIVE MATRIX
C     EB      = PRODUCT OF E*B
C     ELPROP  = ELEMENT ARRAY OF FLOATING PT PROPERTIES
C     FLTMIS  = SYSTEM STORAGE OF FLOATING PT MISC PROP
C     G       = GEOMETRIC INTERPOLATION FUNCTIONS
C     H       = SOLUTION INTERPOLATION FUNCTIONS
C     HINTG   = INTEGRAL OF INTERPOLATION FUNCTIONS
C     LPPROP  = INTEGER PROPERTIES AT EACH ELEMENT NODE
C     LPROP   = ARRAY INTEGER POINT ELEMENT PROPERTIES
C     MISFIX  = MISCELLANEOUS INTEGER SYSTEM PROPERTIES
C     N       = NUMBER OF NODES PER ELEMENT
C     NELFRE  = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
C     NGEOM   = NUMBER OF GEOMETRY NODES
C     NMAT    = NUMBER OF MATERIAL TYPES
C     NPARM   = DIMENSION OF PARAMWETRIC SPACE
C     NQP     = NUMBER OF QUADRATURE POINTS
C     NRB     = NUMBER OF ROWS IN B AND E MATRICES
C     NSPACE  = DIMENSION OF SPACE
C     NTAPE1  = UNIT FOR POST SOLUTION MATRICES STORAGE
C     NTAPE2,3,4 = OPTIONAL UNITS FOR USER (USED WHEN > 0)
C     PRTLPT  = REAL PROPERTIES AT ELEMENT NODES
C     PRTMAT  = REAL ELEM PROPERTIES BASED ON MATERIAL NUMBER
C     PT      = QUADRATURE COORDINATES
C     S       = ELEMENT SQUARE MATRIX
C     STRAIN  = STRAIN OR GRADIENT VECTOR
C     STRAN0  = INITIAL STRAIN OR GRADIENT VECTOR
C     STRESS  = STRESS VECTOR
C     WT      = QUADRATURE WEIGHTS
C     XYZ     = SPACE COORDINATES AT A POINT
C     .....................................................
C      ***  ELSQ PROBLEM DEPENDENT STATEMENTS FOLLOW ***
C     .....................................................
C     2-D POISSON EQUATION VIA 6 NODE TRIANGULAR ELEMENT
C     N=6, NG=1, NSPACE=2, NLPFLO=3
C     ELPROP(1) = CONDUCTIVITY
C     ELPROP(2) = SOURCE PER UNIT VOLUME
C     ELPROP(3) = THICKNESS
      XYK  = ELPROP(1)
      Q    = ELPROP(2)
      T    = ELPROP(3)
      IF ( T .LE. 0.0 )  T = 1.0
C     ZERO ELEMENT SQUARE MATRIX AND SOURCE VECTOR
      CALL ZEROA (NELFRE*NELFRE,S)
      CALL ZEROA (NELFRE,C)
C      NUMERICAL INTEGRATION LOOP
      IF ( NTAPE1 .GT. 0 )  WRITE (NTAPE1)  NQP
      DO 110  IP = 1, NQP
C        LOCAL DERIVATIVE
        CALL  DER6T (PT(1,IP),PT(2,IP),DLH)
C        JACOBIAN AND INVERSE AND DETERMINANT
        CALL  JACOB (N,NSPACE,DLH,COORD,AJ)
        CALL  I2BY2 (AJ,AJINV,DET)
        IF ( DET .LE. 0.0 )  THEN
          WRITE (6,*) 'NEGATIVE JACOBIAN'
          DET = ABS( DET )
        ENDIF
C        PHYSICAL DERIVATIVES
        CALL  GDERIV (NSPACE,N,AJINV,DLH,DGH)
C        LOCAL SHAPE 
        CALL  SHP6T (PT(1,IP),PT(2,IP),H)
        DETWT = DET*WT(IP)*T
C        ELEMENT MATRICES
        DO 90  I = 1, N
          DO 80  J = 1, N
   80     S(I,J) = S(I,J) + DETWT*(XYK*DGH(1,I)*DGH(1,J)
     1                           + XYK*DGH(2,I)*DGH(2,J))
          C(I) = C(I) + DETWT*Q*H(I)
   90   CONTINUE
C        STORE GRADIENTS FOR POST PROCESSING on NTAPE1
        IF ( NTAPE1 .GT. 0 )  THEN
C          PHYSICAL COORDINATES
          CALL MMULT (H, COORD, XYZ, 1, N, NSPACE)
          WRITE (NTAPE1)  XYZ, DGH
        ENDIF
  110 CONTINUE
C      STORE STIFFNESS, AND SOURCE FOR ELEM REACTIONS
      IF ( NTAPE2 .GT. 0 ) WRITE (NTAPE2) S, C
      RETURN
      END
