      SUBROUTINE  ASYMBL ( NG, NCOEFF, MODE, IDIAG, NODES, SS, 
     1    CC, M, NE, NDFREE, NITER, LPTEST, LHOMO, 
     2    NHOMO, NULCOL, N, NSPACE, NELFRE, NRB, NQP, NGEOM,
     3    NPARM, NNPFIX, NNPFLO, MISCFX, MISCFL, NLPFIX, NLPFLO, 
     4    LNODE, INDEX, X, DDOLD, COORD, S, C, H, DGH, B, E, EB, 
     5    STRAIN, STRAN0, STRESS, BODY, PT, WT, XYZ, DLH, G, DLG,
     6    AJ, AJINV, HINTG, D, PRTLPT, FLTNP, FLTEL, FLTMIS,
     7    ELPROP, PRTLPT, PRTMAT, MISFIX, MISFLO, NPFIX, LPFIX, 
     8    LPROP, LPPROP, NTAPE1, NTAPE2, NTAPE3, NTAPE4, NTAPE5 )
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C         ASSEMBLE SYSTEM EQUATIONS AND STORE POST
C                 SOLUTION ELEMENT DATA
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
CDP   IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION  CC(NDFREE), SS(NCOEFF), IDIAG(NDFREE), NODES(NE,N)
C     JUST PASSING THROUGH
C         SYSTEM DATA
      DIMENSION  X(M,NSPACE), DDOLD(NDFREE), LNODE(N), INDEX(NELFRE)
C         SYSTEM PROPERTIES
      DIMENSION  PRTLPT(N,0:NNPFLO), FLTNP(M,0:NNPFLO), 
     1           FLTEL(NE,0:NLPFLO), NPFIX(M,0:NNPFIX), 
     2           LPFIX(NE,0:NLPFIX)
C        FOR USE IN ELSQ OR ELCOL:
      DIMENSION  COORD(N,NSPACE), S(NELFRE,NELFRE)
      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)
      DIMENSION  PT(NPARM,NQP), WT(NQP), XYZ(3), DLH(NSPACE,N), 
     1           G(NGEOM), DLG(NPARM,NGEOM), AJ(NSPACE,NSPACE), 
     2           AJINV(NSPACE,NSPACE), HINTG(N)
      DIMENSION  D(NELFRE), FLTMIS(0:MISCFL), ELPROP(0:NLPFLO), 
     1           PRTMAT(0:NLPFLO), MISFIX(0:MISCFX), 
     2           MISFLO(0:MISCFL), LPROP(0:NLPFIX), 
     3           LPPROP(0:NNPFIX)
C                       VARIABLES:
C     AJ      = JACOBIAN
C     AJINV   = JACOBIAN INVERSE
C     B       = STRAIN-DISPLACEMENT (GRADIENT) MATRIX
C     BODY    = BODY FORCE VECTOR
C     CC      = SYSTEM EQUATIONS COLUMN MATRIX
C     COORD   = SPATIAL COORDINATES OF ELEMENT'S NODES
C     D       = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
C     DDOLD   = SYSTEM NODAL PARAMETERS FROM LAST ITERATION
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 REAL PROPERTIES
C     FLTEL   = REAL PROPERTIES OF ELEMENTS
C     FLTMIS  = MISCELLANEOUS REAL PROPERTIES OF SYSTEM
C     FLTNP   = REAL PROPERTIES OF SYSTEM NODES
C     G       = GEOMETRIC INTERPOLATION FUNCTIONS
C     H       = SOLUTION INTERPOLATION FUNCTIONS
C     HINTG   = INTEGRAL OF INTERPOLATION FUNCTIONS
C     IDIAG   = DIAGONAL LOCATION IN SKYLINE VECTOR
C     INDEX   = SYSTEM DOF NUMBERS ASSOCIATED WITH ELEMENT
C     LHOMO   = 1, IF ELEMENT PROPERTIES ARE HOMOGENEOUS
C     LNODE   = THE N ELEMENT INCIDENCES OF THE ELEMENT
C     LPFIX   = SYSTEM ARRAY OF INTEGER ELEM PROPERTIES
C     LPPROP  = INTEGER PROPERTIES AT EACH ELEMENT NODE
C     LPROP   = ARRAY INTEGER ELEMENT PROPERTIES
C     LPTEST  > 0, IF ELEMENT PROPERTIES HAVE BEEN DEFINED
C     M       = NUMBER OF SYSTEM NODES
C     MODE    = MODE OF STORAGE, 0-SKYLINE, 1-BANDED
C     MISFIX  = MISCELLANEOUS INTEGER SYSTEM PROPERTIES
C     N       = NUMBER OF NODES PER ELEMENT
C     NCOEFF  = TOTAL NUMBER OF TERMS IN SS
C     NDFREE  = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM
C     NE      = NUMBER OF ELEMENTS
C     NELFRE  = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
C     NG      = NUMBER OF PARAMETERS PER NODE
C     NGEOM   = NUMBER OF GEOMETRY NODES
C     NHOMO   = 1, IF NODAL PROPERTIES ARE HOMOGENEOUS
C     NITER   = NO. OF ITERATIONS TO BE RUN (USUALLY 1)
C     NMAT    = NUMBER OF MATERIAL TYPES
C     NODES   = NODAL INCIDENCES OF ALL ELEMENTS
C     NPARM   = DIMENSION OF PARAMWETRIC SPACE
C     NPFIX   = INTEGER PROPERTIES AT ALL NODES
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     NULCOL  > 0, IF ELEMENT COLUMN MATRIX IS ALWAYS ZERO
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     SS      = SYSTEM EQUATIONS SQUARE MATRIX
C     STRAIN  = STRAIN OR GRADIENT VECTOR
C     STRAN0  = INITIAL STRAIN OR GRADIENT VECTOR
C     STRESS  = STRESS VECTOR
C     WT      = QUADRATURE WEIGHTS
C     X       = COORDINATES OF SYSTEM NODES
C     XYZ     = SPACE COORDINATES AT A POINT
C
C     GENERATE ELEMENT EQUATIONS & POST SOLUTION MATRICES
      DO 10  IE = 1, NE
C-->     EXTRACT ELEMENT NODE NUMBERS
        CALL  LNODES (IE,NE,N,NODES,LNODE)
C-->     CALCULATE DEGREE OF FREEDOM NUMBERS
        CALL  INDXEL (N,NELFRE,NG,LNODE,INDEX)
C-->     GENERATE ELEMENT PROBLEM DEPENDENT MATRICES
        CALL  GENELM ( IE, M, NE, NDFREE, NITER, LPTEST, LHOMO, 
     1    NHOMO, NULCOL, N, NSPACE, NELFRE, NRB, NQP, NGEOM,
     2    NPARM, NNPFIX, NNPFLO, MISCFX, MISCFL, NLPFIX, NLPFLO, 
     3    LNODE, INDEX, X, DDOLD, COORD, S, C, H, DGH, B, E, EB, 
     4    STRAIN, STRAN0, STRESS, BODY, PT, WT, XYZ, DLH, G, DLG,
     5    AJ, AJINV, HINTG, D, PRTLPT, FLTNP, FLTEL, FLTMIS,
     6    ELPROP, PRTLPT, PRTMAT, MISFIX, MISFLO, NPFIX, LPFIX, 
     7    LPROP, LPPROP, NTAPE1, NTAPE2, NTAPE3, NTAPE4, NTAPE5 )
C-->     STORE THE MATRICES IN SYSTEM EQUATIONS
        IF ( MODE .EQ. 0 )  THEN
C          SKYLINE VECTOR STORAGE MODE
          CALL  SKYSTR (NCOEFF,NDFREE,NELFRE,INDEX,IDIAG,S,SS)
        ELSE
C          BANDED STORAGE
          MAXBAN = NCOEFF/NDFREE
          CALL  STORSQ (NDFREE,MAXBAN,NELFRE,INDEX,S,SS)
        ENDIF
        IF ( NULCOL .EQ. 0 )  CALL STORCL (NDFREE,NELFRE,INDEX,C,CC)
 10   CONTINUE
C      ASSEMBLY COMPLETED
      RETURN
      END
