C-------------------------------------------------------
C     THIS SUBROUTINE CALCULATE ELEMENT STIFFNESS MATRIX
C-------------------------------------------------------
      SUBROUTINE ELESTF(SE1,SE2,SE3,YONG,RNU,ETHTA,NEM,NOD)
      IMPLICIT REAL*8(A-H,O-Z)
    DIMENSION DCMATX(3,3),DDMATX(2,3)
    DIMENSION SHAPH(NOD),SHAPHD(NOD),ETHTA(NOD)
    DIMENSION BAMAT(3,2*NOD),BBMAT(3,2*NOD),UNMAT(2,2*NOD)
    DIMENSION AKM(2*NOD,2*NOD),BKM1(2*NOD,2*NOD),BKM2(2*NOD,2*NOD)
    DIMENSION SAKM(2*NOD,2*NOD),SBKM(2*NOD,2*NOD),CKM(2*NOD,2*NOD)
    DIMENSION SE1(2*NOD,2*NOD),SE2(2*NOD,2*NOD),SE3(2*NOD,2*NOD)
    DIMENSION DCBA(3,2*NOD),DCBB(3,2*NOD),DDBA(2,2*NOD),DDBB(2,2*NOD)
    DIMENSION BAT(2*NOD,3),BBT(2*NOD,3),UNT(2*NOD,2)
    DIMENSION YONG(NEM),RNU(NEM)
    COMMON/GAUSSP/GAUSP(17),GAUSW(17)
    COMMON/NG/NGAUSP
      COMMON/VAR/IEM
    WRITE(*,*) IEM
    NEVAB=2*NOD
     DO 40 I=1,NEVAB
    DO 40 J=1,NEVAB
    SE1(I,J)=0.0D0
    SE2(I,J)=0.0D0
    SE3(I,J)=0.0D0
   40 CONTINUE
C
C     EVALUATE THE D-MATRIX FOR THE IEMth ELEMENT 
C     BY GAUSS NUMERICAL INTEGRATION
C
    CALL DMATIX(DCMATX,DDMATX,YONG,RNU,NEM)
C
C     DO-LOOP ON NUMBER OF GAUSS POINTS BEGINS HERE
C
      DO 10 IGAUSP=1,NGAUSP
    ETA=GAUSP(IGAUSP)
    WET=GAUSW(IGAUSP)
C
C     CALL SUBROUTINE SHAPEH TO EVALUATE THE INTERPOLATION 
C     FUNCTIONS AND THEIR GLOBAL DERIVATIVES AT THE GAUSS POINT ETA
C
    CALL SHAPEH(SHAPH,SHAPHD,ETA,NOD)
C
C     COMPUTE COEFFICIENT MATRICES
C
    CALL BAMATX(SHAPH,BAMAT,NOD)
    CALL BBMATX(SHAPH,SHAPHD,BBMAT,ETHTA,NOD)
    CALL UNMATX(SHAPH,UNMAT,NOD)
C
C     CALCULATE THE TRANSPOSE BT-MATRIX OF B-MATRIX
C
      DO 50 I=1,2*NOD
    DO 50 J=1,3
    BAT(I,J)=BAMAT(J,I)
    BBT(I,J)=BBMAT(J,I)
   50 CONTINUE
      DO 60 I=1,2*NOD
     DO 60 J=1,2
     UNT(I,J)=UNMAT(J,I)
   60 CONTINUE
C
C     EVALUATE BT*(D*B)-MATRIX
C
    CALL ABMUT(DCMATX,3,3,BAMAT,2*NOD,DCBA)
    CALL ABMUT(BAT,2*NOD,3,DCBA,2*NOD,AKM)
    CALL ABMUT(DCMATX,3,3,BBMAT,2*NOD,DCBB)
    CALL ABMUT(BBT,2*NOD,3,DCBB,2*NOD,CKM)
    CALL ABMUT(BAT,2*NOD,3,DCBB,2*NOD,BKM1)
    CALL ABMUT(BBT,2*NOD,3,DCBA,2*NOD,BKM2)
    CALL ABMUT(DDMATX,2,3,BAMAT,2*NOD,DDBA)
     CALL ABMUT(UNT,2*NOD,2,DDBA,2*NOD,SAKM)
    CALL ABMUT(DDMATX,2,3,BBMAT,2*NOD,DDBB)
    CALL ABMUT(UNT,2*NOD,2,DDBB,2*NOD,SBKM)
C
C     EVALUTE THE NUMERICAL INTERATION
C
     DJ=0.5D0*(ETHTA(NOD)-ETHTA(1))
    DO 70 I=1,2*NOD
     DO 70 J=1,2*NOD
    SE1(I,J)=SE1(I,J)+(AKM(I,J)-2.0D0*SAKM(I,J))*DJ*WET
    SE2(I,J)=SE2(I,J)+(BKM1(I,J)+BKM2(I,J)-2.0D0*SAKM(I,J)
     &-2.0D0*SBKM(I,J))*DJ*WET
    SE3(I,J)=SE3(I,J)+(CKM(I,J)-2.0D0*SBKM(I,J))*DJ*WET
   70 CONTINUE
   10 CONTINUE
      RETURN
    END
C------------------------------------------------------------
C     THIS SUBROUTINE EVALUATE D-MATRIX FOR THE IEMth ELEMENT
C------------------------------------------------------------
      SUBROUTINE DMATIX(DCMATX,DDMATX,YONG,RNU,NEM)
    IMPLICIT REAL*8(A-H,O-Z)
    DIMENSION DCMATX(3,3),DDMATX(2,3),CDMATX(2,3)
    DIMENSION YONG(NEM),RNU(NEM)
    COMMON/PROBM/NTYPE
    COMMON/VAR/IEM
    DO 10 I=1,3
    DO 10 J=1,3
    DCMATX(I,J)=0.0D0
   10 CONTINUE
      DO 20 I=1,2
    DO 20 J=1,3
    DDMATX(I,J)=0.0D0
    CDMATX(I,J)=0.0D0
   20 CONTINUE
      CDMATX(1,1)=1.0D0
    CDMATX(2,3)=1.0D0
    IF(NTYPE.EQ.1) THEN
C
C     D-MATRIX FOR PLANE STRESS CASE
C
      CONST=YONG(IEM)/(1.0D0-RNU(IEM)**2)
    DCMATX(1,1)=CONST
    DCMATX(2,2)=CONST
    DCMATX(1,2)=CONST*RNU(IEM)
    DCMATX(2,1)=CONST*RNU(IEM)
    DCMATX(3,3)=(1.0D0-RNU(IEM))*CONST/2.0D0
    ELSE
C
C     D-MATRIX FOR PLANE STRAIN CASE
C
    CONST=YONG(IEM)*(1.0D0-RNU(IEM))/
     &((1.0D0+RNU(IEM))*(1.0D0-2.0D0*RNU(IEM)))
    DCMATX(1,1)=CONST
    DCMATX(2,2)=CONST
    DCMATX(1,2)=CONST*RNU(IEM)/(1.0D0-RNU(IEM))
    DCMATX(2,1)=CONST*RNU(IEM)/(1.0D0-RNU(IEM))
    DCMATX(3,3)=(1.0D0-2.0D0*RNU(IEM))*CONST/(2.0D0*(1.0D0-RNU(IEM)))
    END IF
C
C     EVALUATE DDMATX BY THE CDMATX*DCMATX
C
      CALL ABMUT(CDMATX,2,3,DCMATX,3,DDMATX)
      RETURN
    END
C-----------------------------------------------------------------
C     THIS SUBROUTINE CALCULATE SHAPE FUNCTION AND ITS DERIVATIVES
C-----------------------------------------------------------------
      SUBROUTINE SHAPEH(SHAPH,SHAPHD,ETA,NOD)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION SHAPH(NOD),SHAPHD(NOD)
C
C     CALCULATE THE QUADRATIC SHAPE FUNCTIONS
C
      N=NOD-2
      SHAPH(1)=0.5D0*(1.0D0-ETA)
    DO 10 I=1,N
      SHAPH(I+1)=ETA**(I-1)*(1.0D0-ETA*ETA)
   10 CONTINUE
    SHAPH(NOD)=0.5D0*(ETA+1.0D0)
C
C     CALCULATE THE DERIVATIVES OF THE QUADRATIC SHAPE FUNCTIONS
C
      SHAPHD(1)=-0.5D0
    SHAPHD(2)=-2.0D0*ETA
    SHAPHD(3)=(1.0D0-ETA*ETA)-2.0D0*ETA*ETA
    DO 20 I=3,N
    SHAPHD(I+1)=DFLOAT(I-1)*ETA**(I-2)*(1.0D0-ETA*ETA)
    &-2.0D0*ETA*ETA**(I-1)
   20 CONTINUE
    SHAPHD(NOD)=0.5D0
    RETURN
    END
C------------------------------------------------
C     THIS SUBROUTINE CALCULATE THE UNMAT-MATRIX
C------------------------------------------------
      SUBROUTINE UNMATX(SHAPH,UNMAT,NOD)
      IMPLICIT REAL*8(A-H,O-Z)
    DIMENSION UNMAT(2,2*NOD),SHAPH(NOD)
    DO 10 I=1,2
    DO 10 J=1,2*NOD
    UNMAT(I,J)=0.0D0
   10 CONTINUE
      DO 20 I=1,NOD
      UNMAT(1,2*I-1)=SHAPH(I)
        UNMAT(2,2*I)=SHAPH(I)
   20 CONTINUE
        RETURN
        END
C------------------------------------------------
C     THIS SUBROUTINE CALCULATE THE BAMAT-MATRIX
C------------------------------------------------
      SUBROUTINE BAMATX(SHAPH,BAMAT,NOD)
    IMPLICIT REAL*8(A-H,O-Z)
    DIMENSION BAMAT(3,2*NOD),SHAPH(NOD)
    DO 10 I=1,3
    DO 10 J=1,2*NOD
    BAMAT(I,J)=0.0D0
   10 CONTINUE
      DO 20 I=1,NOD
      BAMAT(1,2*I-1)=SHAPH(I)
      BAMAT(3,2*I)=SHAPH(I)
   20 CONTINUE
    RETURN
    END
C------------------------------------------------
C     THIS SUBROUTINE CALCULATE THE BBMAT-MATRIX
C------------------------------------------------
      SUBROUTINE BBMATX(SHAPH,SHAPHD,BBMAT,ETHTA,NOD)
    IMPLICIT REAL*8(A-H,O-Z)
    DIMENSION BBMAT(3,2*NOD),SHAPH(NOD),SHAPHD(NOD),ETHTA(NOD)
    COMMON/VAR/IEM
      DJ=0.5D0*(ETHTA(NOD)-ETHTA(1))
    DO 10 I=1,3
    DO 10 J=1,2*NOD
    BBMAT(I,J)=0.0D0
   10 CONTINUE
      DO 20 I=1,NOD
      BBMAT(1,2*I-1)=SHAPH(I)
    BBMAT(2,2*I-1)=SHAPH(I)
    BBMAT(2,2*I)=SHAPHD(I)/DJ
    BBMAT(3,2*I-1)=SHAPHD(I)/DJ
   20 CONTINUE
    RETURN
    END
C------------------------------------------------------------
C     THIS SUBROUTINE MULTIPLIES THE A-MATRIX BY THE B-MATRIX
C------------------------------------------------------------
      SUBROUTINE ABMUT(AA,N,M,BB,L,CC)
    IMPLICIT REAL*8(A-H,O-Z)
    DIMENSION AA(N,M),BB(M,L),CC(N,L)
    DO 10 I=1,N
    DO 10 J=1,L
    CC(I,J)=0.0D0
    DO 10 K=1,M
    CC(I,J)=CC(I,J)+AA(I,K)*BB(K,J)
   10 CONTINUE
      RETURN
    END
C------------------------------------------------------------    
C     THIS SUBROUTINE CALCULATES THE GAUSS POINT COORDINATES 
C     AND WEIGHTS.
C------------------------------------------------------------
      SUBROUTINE GAUSS(X1,X2,XX,W,N)
    IMPLICIT REAL*8(A-H,O-Z)
    DIMENSION XX(N),W(N)
    PARAMETER (EPS=1.0D-15)
    M=(N+1)/2
    XM=0.5D0*(X1+X2)
    XL=0.5D0*(X2-X1)
    DO 12 I=1,M
    Z=DCOS(3.14159265359D0*(I-0.25D0)/(N+0.5D0))
  1    CONTINUE
    P1=1.0D0
    P2=0.0D0
    DO 11 J=1,N
    P3=P2
    P2=P1
    P1=((2.0D0*J-1.0D0)*Z*P2-(J-1.0D0)*P3)/J
 11    CONTINUE
    PP=N*(Z*P1-P2)/(Z*Z-1.0D0)
    Z1=Z
    Z=Z1-P1/PP
    IF (DABS(Z-Z1).GT.EPS) GO TO 1
    XX(I)=XM-XL*Z
    XX(N+1-I)=XM+XL*Z
    W(I)=2.0D0*XL/((1.0D0-Z*Z)*PP*PP)
    W(N+1-I)=W(I)
 12    CONTINUE
    RETURN
    END