文章标题

CCC                               
      PROGRAM   FEATP
C====================================================================
C        THIS PROGRAM CAN SOLVE THE ELASTIC PROBLEM OF PLANE STRESS,
C            PLANE STRAIN,AXISYMMETRIC SOLID AND MINDLIN PLATE
C====================================================================
C        7 SUBROUTINES(ALLOCAT,INPUT,ASSEM,STATIC_SOLVE,STRESS,DYNAM,
C            EIGEN)ARE CALLED BY THE MAIN PROGRAM.
C        BESIDES,ANOTHER 18 SUBROUTINES ARE CALLED BY ABOVE 7 SUBROUTINES.
C=======================================================================
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER*5 FILENAME
    DIMENSION IZ(300000),AR(15000000)
C                !IZ-MAXIMUN SPACE FOR DYNAMIC INTEGE ARRAY
C                !AR-MAXIMUM SPACE FOR DYNAMIC REAL ARRAY
      COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI
      COMMON/COMN/NFIX,NPC,GRAV
      COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA
      COMMON/COM3/MND2,NUMPT2
      COMMON/ELEM/NODE,INTX,INTY
      COMMON/DYN/HUV,OMEGA,CC1,CC2,TT,DT,ALPHA,DELTA   
      print*,'===> Please Input the Title of the Problem to be Solved!'
      READ(*,'(a5)')FILENAME
    OPEN(5,FILE=FILENAME//'.inp',STATUS='OLD')
      OPEN(6,FILE=FILENAME//'.dat',STATUS='UNKNOWN')
      OPEN(10,FILE=FILENAME//'.mkp',STATUS='UNKNOWN')
      OPEN(14,FILE=FILENAME//'.dis',STATUS='UNKNOWN')
      OPEN(15,FILE=FILENAME//'.str',STATUS='UNKNOWN')
C*********************************************************************************************************************************
C        ALLOCATE STORAGE SPACE FOR THE ARRAYS OF FE MODEL
C*****************************************************************************************************
      CALL ALLOCAT(M1,M2,M3,M4,M5,M8,M11,M12,N1,N2,N3,N4,N5,N6,N7,
     $       N8,N9,N10,N11,N12,N13,N14,N15,N16,N17,N18,N19,N20,N21,N22,
     $       N23,N24,N25,N26,N27,N31,N32,N33,N34,N35,N36,N37,N38)
C***************************************************************************************************
C          INPUT ALL THE DATA OF FE MODEL AT APPOINTED ADDRESS
C****************************************************************************************************
      CALL INPUT(IZ(M2),IZ(M3),IZ(M4),AR(N1),AR(N2),AR(N3),AR(N4))
C****************************************************************************************************
C          ASSEMBLE THE ELEMENT MATRIXES TO FORM GLOBAL MATRIXES:[M],[K],{P}
C*******************************************************************************************************
      CALL ASSEM(AR(N1),IZ(M2),AR(N2),IZ(M3),AR(N3),IZ(M4),
     $               AR(N4),AR(N11),AR(N12),AR(N13),AR(N14),
     $               IZ(M11),AR(N15),AR(N16),AR(N17))
C********************************************************************************************************
C          GET NODAL DISPLACEMENTS BY SOLVING EQUATION[K]{U}={P}
C                    FOR STATIC PROBLEM(MSOLVE=1)
C********************************************************************************************************
      IF(MSOLV.EQ.1)THEN
      CALL STATIC_SOLVE(AR(N12),AR(N13),AR(N14))
C*********************************************************************************************************
C            COMPUTE STRESSES AT THE INTEGRATION POINTS AND NODES 
C***********************************************************************************************************
      CALL STRESS(AR(N1),IZ(M2),AR(N2),AR(N14),AR(N15),IZ(M11),
     $       AR(N18),AR(N9),AR(N10),IZ(M8),AR(N8))
      ENDIF
C****************************************************************************************************
C       SOLVE DYNAMIC RESPONSE PROBLEM BY CENTRAL-DIFFERENCE METHOD 
C                                 (MSOLVE=2) AND BY NEWMARK METHOD(MSOLV=3)
C*******************************************************************************************************
      IF(MSOLV.EQ.2.OR.MSOLV.EQ.3)
     $CALL DYNAM(AR(N11),AR(N12),AR(N13),AR(N21),AR(N22),
     $         AR(N23),AR(N24),AR(N25),AR(N26))
C************************************************************************************************************
C       SOLVE DYNAMIC CHARACTER PROBLEM BY THE INVERSE ITERATION METHOD 
C             (MSOLV=4) AND BY THE SUBSPACE ITERATION METHOD(MSOLV=5)
C************************************************************************************************************ 
      IF(MSOLV.EQ.4.OR.MSOLV.EQ.5)
     $CALL EIGEN(AR(N11),AR(N12),AR(N31),AR(N32),AR(N33),
     $         AR(N34),AR(N35),AR(N36),AR(N37))
      STOP
      END
C=========================================================================
C=======================SUB:1=============================================
      SUBROUTINE ALLOCAT(M1,M2,M3,M4,M5,M8,M11,M12,N1,N2,N3,N4,N5,
     $  N6,N7,N8,N9,N10,N11,N12,N13,N14,N15,N16,N17,N18,N19,N20,
     $  N21,N22,N23,N24,N25,N26,N27,N31,N32,N33,N34,N35,N36,N37,N38)       

C***********************************************************************************************************
C   INPUT BASIC PAPRMETERS FROM FILE'IN_DAT'
C   ALLOCAT DYANMICAL STORAGE SPACE
C***********************************************************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI
      COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA
      COMMON/COMN/NFIX,NPC,GRAV
      COMMON/COM3/MND2,NUMPT2
C                      !MND-MAXIMAL NODE NUMBER IN ALL ELEMENTS
C                      !NUMEL-NEMBER OF GLOBAL ELEMENTS 
C                      !NUMPT-NUMBER OF GLOBAL NODES
C                      !MBAND-HALF BANDWIDTH(INCLUDING DIAGONAL ELEMENT)
      READ(5,*)MND,NUMEL,NUMPT,MBAND
      READ(5,*)NFIX,NPC,MPROB,MSOLV
C                      !NFIX-NUMBER OF NODES SUBJECTER TO CONSTRIANT
C                      !NPC-NUBMER OF NODES SUBJECTED TO EQUIVALENT LOAD
C==========================================================
C       MPROB=1-PLANE STRESS PROBLEM,  MPROB=2-PLANE STRAIN PROBLEM
C       MPROB=3-AXISYMMETRIC PROBLEM,   MPROB=4-MINDLIN PLATE PROBLEM
C------------------------------------------------------------------------------------------------------------------------------
C       MSOLV=1-STATIC ANALYSIS
C       MSOLV=2-DYNAMIC RESPONSE ANALYSIS BY CENTRAL DIFFERENCE METHOD
C       MSOLV=3-DYNAMIC RESPONSE ANALYSIS BY NEWMARK METHOD
C       MSOLV=4-DYNAMIC CHARACTER ANALYSIS BY INVERSE ITERATION METHOD
C       MSOLV=5-DYNAMIC CHARACTER ANALYSIS BY SUBSPACE ITERATION METHOD
C==================================================================
      IF(MSOLV.NE.4.OR.MSOLV.NE.5)READ(5,*)NMATI,GRAV,MTYPE
      IF(MSOLV.EQ.4.OR.MSOLV.EQ.5)READ(5,*)NMATI,GRAV,MTYPE,NVA
C                               !NMATI-KIND OF MATERIALS
C                               !GRAV-GRAVITY ACCELERATION
C                               !NVA-NUMBER OF EIGENVALUES
C----------------------------------------------------------------------------------------------------------------------------------------
C            MTYPE-CONTROL KEY FOR OUTPUT RESULTS 
C            MTYPE=0-OUTPUT RESULTS INGLUDE GLOBAL MATRIXES AND STRESS AT 
C                         GAUSS POINTS
C            MTYPE=1-OUTPUT RESULTS INCLUDE GLOBAL MATRIXES 
C            MTYPE=2-OUTPUT RESULTS INCLUDE STRESSES AT GAUSS POINTS
C-------------------------------------------------------------------------------------------------------------------------------------------
      IF(MPROB.EQ.1.OR.MPROB.EQ.2)   THEN
      NF=2         !NF-NUMBER OF NODAL FREEDOMS(DISPLACEMENT COMPONENTS)
      NFSTR=3      !NFSTR-NUMBER OF STRESS COMPONENTS
      ENDIF
      IF(MPROB.EQ.3)THEN
      NF=2
      NFSTR=4
      ENDIF
      IF(MPROB.EQ.4)THEN
      NF=3
      NFSTR=5
      ENDIF
C*************************************************************************************************************************
C               ALLOCATE STORAGE SPACE FOR INPUT DATA
C**********************************************************************************************************************
      M1=1
      M2=M1
      M3=M2+NUMEL*14       !IELEM(NUMEL,14)-CODE OF MATERIAL AND NODES
      M4=M3+NFIX*4         !IFIXD(NFIX,4)-CODE OF NODES HAMING CONSTRIANT
      M5=M4+NPC*4          !ILOAD(NPC,4)-CODE OF NODES HAVING LOAD
      N1=1
      N2=N1+NMATI*4        !VAMATI(NMAT1,4)-PARAMETERS OF MATERIALS
      N3=N2+NUMPT*2        !VCOOD(NUMPT,2)-GLOBAL NODE COORDINATES
      N4=N3+NFIX*3         !VFIXD(NFIX,3)-CONSTRIANED DISPLACEMENTS
      N5=N4+NPC*3          !VLOAD(NPC,3)-VALUES OF EQUIVALENT LOAD
C****************************************************************************************************************************
C                  ALLOCATE STORAGE SPACE FOR ELEMENTAL MATRIXES AND GLOBAL MATRIXES
C******************************************************************************************************************************
      MND2=MND*NF          !MND2-NUMBER OF FREEDOMS IN A ELEMENT
      NUMPT2=NUMPT*NF      !NUMPT2-NUMBER OF GLOBAL FREEDOMS
      M8=M5
      M11=M8+NUMPT         !IADD(NUMPT)-USED FOR NODAL STRESS
      M12=M11+MND          !IEL(MND)-NODE CODE IN A ELEMENT
      N6=N5
      N7=N6+NFSTR*MND2     !VSG(NFSTR,MND2)-ELEMENT STRESS MATRIX AT 
                         !INTEGRATION POINTS
      N8=N7+NFSTR*MND2     !VSN(NFSTR,MND2)-ELEMENT STRESS AT NODES 
      N9=N8+NUMPT*NFSTR    !SSN(NUMPT,NFSTR)-STRESSES AT NODES 
      N10=N9+9*NFSTR       !SSS(9,NFSTR)-STRESSES AT INTERGARION POINTS
    N11=N10+4*NFSTR      !VSS(4,NFSTR)-STRESSES AT DESIGNATED
                                !INTEGRATION POINTS`
      N12=N11+NUMPT2*MBAND      !GMM(NUMPT2,MBAND)-GLOBAL MASS MATRIX
      N13=N12+NUMPT2*MBAND      !GKK(NUMPT2,MBAND)-GLOBAL STIFENESS MATRIX
      N14=N13+NUMPT2                !GP(NUMPT2)-GLOBAL LOAD VECTOR
      N15=N14+NUMPT2                !GU(NUMPT2)-GLOBAL DISPLACEMENT VECTOR
      N16=N15+MND*2                 !VXY(MND,2)-NODE COORDINATE IN ELEMENT
      N17=N16+MND2*MND2             !VMM(MND2,MND2)-ELEMENT MASS MATRIX
      N18=N17+MND2*MND2             !VKK(MND2,MND2)-ELEMENT STIFFNESS MATRIX
      N19=N18+MND2                  !VU(MND2)-NODE DISPLACEMENTS IN A ELEMENT
      N20=N19+4                     !VMAE(4)-MATERIAL PARAMETERS IN A ELEMENT
C******************************************************************************************************************************
C                      ALLOCAT STORAGE SPACE FOR DYNAMIC RESPONSE ANALYSIS
C*********************************************************************************************************************************
      IF(MSOLV.EQ.2.OR.MSOLV.EQ.3)THEN
      N21=N20
      N22=N21+NUMPT2*MBAND           !DAMP(NUMPT2,MBAND)-DAMPLING MATRIX
      N23=N22+NUMPT2                 !U0(NUMPT2)-INITIAL DISPLACEMENT VECTOR
      N24=N23+NUMPT2                 !V0(NUMPT2)-INITIAL VELOCITY VECTOR
      N25=N24+NUMPT2                 !A(NUMPT2)-INITIAL ACCELERATION VECTOR
      N26=N25+NUMPT2*MBAND           !AW(NUMPT2,MBAND)-WORKING ARRAY
      N27=N26+NUMPT2                 !B(NUMPT2)-WORKING ARRAY
      ENDIF
C*************************************************************************************************************************************
C                      ALLOCAT STORAGE SPACE FOR DYNAMIC CHARACTER ANALYSIS
C***************************************************************************************************************************************
      IF(MSOLV.EQ.4.OR.MSOLV.EQ.5)THEN
      IF(NVA.EQ.0)THEN
      WRITE(*,11) 
 11   FORMAT('PLEASE READ THE NUMBERS OF EIGENVALUE-NVA=')
      READ(*,*)NVA
      ENDIF
      N31=N20
      N32=N31+NUMPT2*NVA         !AA(NUMPT2,NVA)-INITIAL ITERATION VECTOR
      N33=N32+NUMPT2*NVA         !BB(NUMPT2,NVA)-WORKING ARRAY
      N34=N33+NVA*NVA            !GM(NVA,NVA)-MASS MATRIX IN SUBSPACE
      N35=N34+NVA*NVA            !GK(NVA,NVA)-STIFFNESS MATRIX IN SUBSPACE
      N36=N35+NVA*NVA            !V(NVA,NVA)-EIGENVECTORS IN SUBSPACE
      N37=N36+NVA                !W1(NVA)-WORKING ARRAY
      N38=N37+NVA                !W2(NVA)-EIGENVALUES IN SUBSPACE
      ENDIF
      RETURN 
      END
C=======================SUB:2=========================================================
      SUBROUTINE INPUT(IELEM,IFIXD,ILOAD,VMATI,VCOOD,VFIXD,VLOAD)
C*****************************************************************************************************************************
C               INPUT ALL THE INFORMATION OF FE MODEL
C********************************************************************************************************************************
      IMPLICIT  REAL*8(A-H,O-Z)
      COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI
      COMMON/COMN/NFIX,NPC,GRAV
      COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA
      COMMON/DYN/HUV,OMEGA,CC1,CC2,TT,DT,ALPHA,DELTA
      DIMENSION  IELEM(NUMEL,14),IFIXD(NFIX,4),ILOAD(NPC,4),
     $           VCOOD(NUMPT,2),VFIXD(NFIX,3),VLOAD(NPC,3),
     $           VMATI(NMATI,5)
      WRITE(*,101) 
101   FORMAT(/6X,'## IMPUT DATA FROM FILE<IN_DAT>TO MEMORY #')
C========================================================================================
                       !INPUT NODAL COORDINATES
      READ(5,*)(II,(VCOOD(I,J),J=1,2),I=1,NUMPT)
                         !INPUT ELEMENT INFORMATION
      READ(5,*)(II,(IELEM(I,J),J=1,4+MND),I=1,NUMEL)
                        !IMPUT CONSTRAIN INFORMATION
      READ(5,*)(II,(IFIXD(I,J),J=1,NF+1),(VFIXD(I,J),J=1,NF),I=1,NFIX)

      IF(NPC.GT.0)THEN
                      !INPUT  EQUIVALENT LOAD AT NODES
      READ(5,*)(II,(ILOAD(I,J),J=1,NF+1),(VLOAD(I,J),J=1,NF),I=1,NPC)

      ENDIF
                      !INPUT MATERIAL PARAMETERS
      READ(5,*)(II,(VMATI(I,J),J=1,3),I=1,NMATI)
      IF(MSOLV.EQ.2.OR.MSOLV.EQ.3) THEN
      READ(5,*) 
      READ(5,*)       
      READ(5,*) HUV,FREQ,CC1,CC2,TT,DT,ALPHA,DELTA
C                            !HUV-INPUT CONTROL OF INITIAL DISPLACEMENT AND VELOCITY
C                            (HUV=1:INPUT DISPLACEMENT;HUV=2:INPUT VELOCITY;
C                             HUV=3:INPUT BOTH OF THEM)AD
C                            !OMEGA-FREQUENCE OF LO 
C                            !CC1,CC2-CONSTANTS TO COMPUTE DAMPLING MATRIXS
C                            !TT-TOTAL TIME
C                            !DT-TIME STEP LENGTH
C                            !ALPHA,DELTA-CONSTANTS USED IN THE NEWMARK METHOD
      ENDIF
C***************************************************************************************
C                         OUTPUT ABOVE INPUT DATA
C****************************************************************************
      WRITE(*,102)
102   FORMAT(/5X,'%% OUTPUT INPUT-DATA TO<OUT_DAT> %%')
      WRITE(6,10)
10    FORMAT(/8X,'MAXIMAL-ELEM-NODES,ELEMENTS,NODES,BANDWIDTH')
      WRITE(6,11)    MND,NUMEL,NUMPT,MBAND
11    FORMAT(10X,I10,I11,I9,I9)
      WRITE(6,12)
12    FORMAT(/8X,'FIXED-NODES,EQUIVALENT-LOADS,MATERIAL-KINDS,',
     $                  2X,'GRAVITY')
      WRITE(6,13) NFIX,NPC,NMATI,GRAV
13    FORMAT(11X,I9,2(3X,I9),2X,F16.5)
      WRITE(6,14)
14    FORMAT(/8X,'PROBLEM-KIND,SOLVING-KIND,OUTPUT-KEY')
      WRITE(6,15) MPROB,MSOLV,MTYPE
15    FORMAT(11X,I9,2(3X,I9))
      WRITE(6,16)
16    FORMAT(/8X,'NODAL-FREEDOMS,STRESS-COMPONENTS')
      WRITE(6,17) NF,NFSTR
17    FORMAT(11X,I9,3X,I10)
      WRITE(6,18)                      ! OUTPUT NODAL COORDIATES
18    FORMAT(/8X,'NODAL COORDINATES'/9X'NO.',15X,'X-',11X,'Y-')
      DO 19 I=1,NUMPT
19    WRITE(6,20)I,(VCOOD(I,J),J=1,2)
20    FORMAT(5X,I5,5X,3F15.6)
      WRITE(6,21)(I,I=1,9)             ! OUTPUT ELEMENT INFORMATION      
21    FORMAT(/8X,'ELEMENT  INFORMATION'/3X,'NP.',1X,'NODES',1X,
     $           'MATERIAL',1X,'INTX',1X,'INTY',1Z,9(2X,2HN-,I1))
      DO 22 I=1,NUMEL
22    WRITE(6,23) I,(IELEM(I,J),J=1,MND+4)
23    FORMAT(1X,I5,2I5,3X,2I5,3X,9(1X,I4))
      WRITE(6,24)                      ! OUTPUT XONSTRAIN INFORMATION
24    FORMAT(/8X,'CONSRAINT INFORMATION ON NODES'/
     $            11X,'NO.',2X,'NODE NO.',2X,'X-',1X,'Y-',1X,'Z-',
     $            5X,'X-VALUE',2X,'Y-VALUE',2X,'Z-VALUE')
      DO 25 I=1,NFIX
25    WRITE(6,26) I,(IFIXD(I,J),J=1,4),(VFIXD(I,J),J=1,3)
26    FORMAT(10X,I4,4X,I3,3X,3I3,3X,3E10.3)
      IF(NPC.NE.0) THEN
      WRITE(6,27)                      !  OUTPUT EQUIVALENT LOAD AT NODES
27    FORMAT(/8X,'EQUIVALENT LOAD ON NODES'/12X,'NO.',1X,
     $        'NODE NO.',2X,'X-',1X,'Y-',1X,'Z-',5X,'X-VALUE',2X,
     $        'Y-VALUE',2X,'Z-VALUE')
      DO 28 I=1,NPC
28    WRITE(6,29) I,(ILOAD(I,J),J=1,4),(VLOAD(I,J),J=1,3)
29    FORMAT(10X,I4,4X,I3,4X,3I3,3X,3E10.3)
      ENDIF
      WRITE(6,33)
33    FORMAT(/8X,'MATERIAL PARAMETERS'/5X,'NO.'3X,'E(MODULUS)',2X,
     $      'V(POISSON RATIO)',2X,'DENS(DENSITY)',2X,'TH(THICKNESS)')
      DO 34 I=1,NMATI
34    WRITE(6,35) I,(VMATI(I,J),J=1,4)
35    FORMAT(5X,I2,4E14.3)
C*******
C********OUTPUT THE INPUT PARAMETERS FOR DYNAMIC REPONSE COMPUTATION
C**********
      IF(MSOLV.EQ.2.OR.MSOLV.EQ.3)
     $WRITE(6,40) HUV,OMEGA,CC1,CC2,TT,DT,ALPHA,DELTA
40    FORMAT(/8X,'PARAMETERS FOR DYNAMIC RESPONSE COMPUTATION:'/
     $     5X,'INPUT CONTROL FOR INITIAL DISPLACE AND VELOCITY-HUV=',I3/
     $     5X,'FREQUENCE OF LOAD-OMEGA=',E10.3/
     $     5X,'DAMPING COEFFICIENTS-CC1=',E10.3,8X,'CC2=',E10.3/
     $     5X,'TOTAL TIME-TT=',E12.6,10X,'STEP LENGTH-DT=',E12.6/
     $     5X,'PARAMETERS OF NEWMARK-ALPHA=',E10.3,8X,'DELTA=',E10.3)
C**********
C*******OUTPUT THE INPUT PARAMETER FOR CHARACTER VALUE COMPUTATION
C*********
      IF(MSOLV.EQ.4.OR.MSOLV.EQ.5)  WRITE(6,45) NVA
45    FORMAT(/8X,'PARAMERER FOR DYNAMIC CHARACTER COMPUTATION:'/
     $                         10X,'NUMBERS OF EIGENVALUES-NVA=',I3)
      WRITE(6,301)
301   FORMAT(/'C************************************************',
     $                '******************************',/)
      RETURN
      END
C=========================SUB:3=====================================
      SUBROUTINE ASSEM(VMATI,IELEM,VCOOD,IFIXD,VFIXD,
     $             ILOAD,VLOAD,GMM,GKK,GP,GU,IEL,VXY,VMM,VKK)
C*******************************************************************
C                    ASSEMBLE THE GLOBAL MATRIXES:[M],[K],AND{P}
C                          CALL ELEMENT-MATRIX
C*****************************************************************
      IMPLICIT  REAL*8(A-H,O-Z)
      COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI
      COMMON/COMN/NFIX,NPC,GRAV
      COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA
      COMMON/COM3/MND2,NUMPT2
      COMMON/ELEM/NODE,INTX,INTY
      DIMENSION   IELEM(NUMEL,MND+4),VCOOD(NUMPT,2),IFIXD(NFIX,4),
     $            VFIXD(NFIX,3),ILOAD(NPC,4),VLOAD(NPC,3),
     $            GMM(NUMPT2,MBAND),GKK(NUMPT2,MBAND),GP(NUMPT2),
     $            GU(NUMPT2),VXY(MND,2),IEL(MND),VMATI(NMATI,4),
     $            VMAE(4),VMM(MND2,MND2),VKK(MND2,MND2)
      WRITE(*,101)
101   FORMAT(/5X,'## ASSEMBLE GLOBAL MATRIX[GKK],[GMM]##')
C*********
C**********GLEAR THE MEMORY FOR GLOBAL MATRIX:[M]AND[K]
C********
      DO 10 I=1,NUMPT2
      DO 10 J=1,MBAND
      GKK(I,J)=0.0
      GMM(I,J)=0.0
10    CONTINUE
C*******
C*********LOOP OVER EACH ELEMENT
C*********
      DO 320 IE=1,NUMEL          !NUMEL-NUMBER OF ELEMENTS
C********
C*********FORM INFORMATION OF EACH ELEMENT FROM THE INPUT DATA
C*******
      DO 11 I=1,MND              !MND-MAXIMUM NODE NUMBER IN ALL ELEMENTS
      IEL(I)=IELEM(IE,I+4)
      DO 11 J=1,2
      VXY(I,J)=0.0
      IF(IEL(I).GT.0) VXY(I,J)=VCOOD(IEL(I),J)
11    CONTINUE
      NODE=IELEM(IE,1)
      INTX=IELEM(IE,3)
      INTY=IELEM(IE,4)
      IMATI=IELEM(IE,2)
      DO 13 J=1,4
      VMAE(J)=VMATI(IMATI,J)
13    CONTINUE
C****************************************************************************
C             COMPUTE ELEMENTAL MASS MATRIX[VMM] AND STIFFNESS MATRIX[VKK]
C                       FOR THE ELEMENTS WITH 3-6,4-8,9 NODES
C********************************************************************************
C 
      CALL ELEMENT_MATRIX(IE,VXY,IEL,VMM,VKK,VMAE)
C
C*************************************************************************
C             ASSEMBLE ELEMENT MATRIXES TO FORM THE GLOBAL MASS MATRIX[GMM]
C                       AND GLOBAL STIFFNESS MATRIX[GKK]
C*****************************************************************************
      DO 20 I=1,MND
      DO 20 J=1,MND
      DO 25 II=1,NF
      DO 25 JJ=1,NF
      IH=NF*(I-1)+II
      IV=NF*(J-1)+JJ
      IHH=NF*(IEL(I)-1)+II
      IVV=NF*(IEL(J)-1)+JJ
      IVV=IVV-IHH+1
      IF(IHH.GT.0.AND.IVV.GT.0)
     $               GKK(IHH,IVV)=GKK(IHH,IVV)+VKK(IH,IV)
      IF(IHH.GT.0.AND.IVV.GT.0)
     $               GMM(IHH,IVV)=GMM(IHH,IVV)+VMM(IH,IV)
25    CONTINUE
20    CONTINUE
320   CONTINUE
      WRITE(*,102)
102   FORMAT(/6X,'## DRAW BOUNDARY COMNDITIONS TO FORM [GP]##')
C******************************************************************************
C                         FORM THE GRAVITY LOADING [P]
C*******************************************************************************
      DO 30 I=1,NUMPT
      DO 30 J=1,NF
      GP(NF*(I-1)+J)=0.0
      GU(NF*(I-1)+J)=0.0
      IF(J.EQ.NF)      GU(NF*(I-1)+J)=GRAV          !GRAV-GRAVITY ACCELERATION
30    CONTINUE
      DO 31 I=1,NUMPT2
      GP(I)=GMM(I,1)*GU(I)
      DO 31 K=I+1,I+MBAND-1
      IF(K.LE.NUMPT2) GP(I)=GP(I)+GMM(I,K-I+1)*GU(K)
      IF(2*I-K.GE.1) GP(I)=GP(I)+GMM(2*I-K,K-I+1)*GU(2*I-K)
31    CONTINUE
      IF(GRAV.NE.0.0) THEN
      WRITE(6,*)'LOAD UNDER GRAVITY:'
      DO 32 I=1,NUMPT
32    WRITE(6,74) I,(GP(NF*(I-1)+J),J=1,NF)
      WRITE(6,301)
      ENDIF
C*********************************************************************************
C                         ADD THE NODAL LOAD VECTOR[P]
C********************************************************************************
      DO 40 I=1,NPC
      DO 40 J=1,NF
      IF(ILOAD(I,J+1).NE.0) THEN
      II=NF*(ILOAD(I,1)-1)+J
      GP(II)=GP(II)+VLOAD(I,J)
      ENDIF
40    CONTINUE
C************************************************************************************
C                         DRAW THE NODAL CONSTRAINT
C***********************************************************************************
      IF(MSOLV.EQ.1) THEN
      DO 42 I=1,NFIX
      DO 42 J=1,NF
      IF(IFIXD(I,J+1).NE.0) THEN
      II=NF*(IFIXD(I,1)-1)+J
      GU(II)=VFIXD(I,J)*1E30
      GKK(II,1)=GKK(II,1)*1E30
      IF(GKK(II,1).LE.1E-10) GKK(II,1)=0.
      ENDIF
42    CONTINUE
      ENDIF
      IF(MSOLV.GT.1) THEN
      DO 50 I=1,NFIX                                                      
      DO 50 J=1,NF
      IF(IFIXD(I,J+1).NE.0) THEN
      II=NF*(IFIXD(I,1)-1)+J
      GU(II)=0.0
      GKK(II,1)=1.0
      GMM(II,1)=1.0
      IF(MSOLV.EQ.4.OR.MSOLV.EQ.5) GMM(II,1)=0.0
      DO 60 K=2,MBAND
      GKK(II,K)=0.0
      GMM(II,K)=0.0
      IF(II.LT.K) GOTO 60
      GKK(II-K+1,K)=0.0
      GMM(II-K+1,K)=0.0
60    CONTINUE
      ENDIF
50    CONTINUE
      ENDIF
C****************************************************************************************
C               OUTPUT THE MATRIXES:[M],[K] AND{P} TO FILE 'OUT_MKP'
C                                (WHEN MTYPE=0 OR MTYPE=1)
C***************************************************************************************
      CLOSE(10,STATUS='DELETE')
      IF(MTYPE.EQ.0.OR.MTYPE.EQ.1) THEN
      WRITE(*,105)
105   FORMAT(/6X,'%% OUTPUT GLOBAL MATRIX INTO<OUT_MKP>%%')
C*******
C****************OUTPUT THE GLOBAL MASS MATRIX:[M]
C********
      OPEN(10,FILE='OUT_MKP',STATUS='UNKNOWN')
      WRITE(10,301)
      WRITE(10,*)'GLOBAL MASS MATRIX(GMM):'
      DO 70 I=1,NUMPT
      DO 70 J=1,NF
      II=NF*(I-1)+J
      WRITE(10,71) I,II,(GMM(II,K),K=1,MBAND)                            
70    CONTINUE
71    FORMAT(2I5,20(3X,2E15.6))
C***********
C******************OUTPUT THE GLOBAL STIFFNESS MATRIX[K]
C*********
      WRITE(10,301)
      WRITE(10,*) 'GLOBAL STIFFNESS MATRIX (GKK):'
      DO 72 I=1,NUMPT
      DO 72 J=1,NF
      II=NF*(I-1)+J
      WRITE(10,71) I,II,(GKK(II,K),K=1,MBAND)
72    CONTINUE
C*********
C******************OUTPUT THE GLOBAL LOADING{P}
C********
      WRITE(10,301)
      WRITE(10,*)'GLOBAL LOADING(GP):'
      DO 73 I=1,NUMPT
      WRITE(10,74) I,(GP(NF*(I-1)+J),J=1,NF)
73    CONTINUE
74    FORMAT(5X,I6,4X,3E15.6)
      WRITE(10,301)
301   FORMAT(/3X,'********************',
     $               '***************************',/)
      CLOSE(10)
      ENDIF
      RETURN
      END
C=========================SUB:3-1==========================================================
      SUBROUTINE ELEMENT_MATRIX(IE,VXY,IEL,VMM,VKK,VMAE)
C***************************************************************************************
C          FORM THE ELEMENT MATRIXES:[M],[K] AND[S]
C          CALL SUBROUTINES OF ELEMENT_VD,GAUSS_INTEGRAATION OR
C               HAMMER_INTEGATION,SHAPE_QUADRANGLE_8 OR SHAPE_QUADRANGLE_9
C               OR SHAPE_TRIANGLE,ELEMENT_VB,ELEMENT_JOCABI
C*******************************************************************************************
      IMPLICIT  REAL*8(A-H,O-Z)
      COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI
      COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA
      COMMON/COM3/MND2,NUMPT2
      COMMON/ELEM/NODE,INTX,INTY
      DIMENSION  VXY(MND,2),IEL(MND),VMAE(4),VMM(MND2,MND2),
     $           VKK(MND2,MND2),VSG(NFSTR,MND2)
      DIMENSION  VN(9),VDN(3,9),VD0(3,9),VD(5,5),VB(5,27)
C*******
C**********FORM THE[D] MATRIX ACCORDING TO TYPE OF THE PROBLEM
C********
      CALL ELEMENT_VD(MPROB,VD,VMAE)      !VD-ELASTIC MATRIX
C********
C**********CLEAR THE MEMORY FOR ELEMENT MATRIXES[M]AND[K]
C********
      DO 10 I=1,MND2
      DO 10 J=1,MND2
      VMM(I,J)=0.0
      VKK(I,J)=0.0
10    CONTINUE
C**********************************************************************************************
C          COMPUTE SHAPE FUNCTION VN AND ITS LOCAL DERIVATIVE VDN AT EACH
C                    INTEGRATOPM POINTS IN ELEMENTS
C******************************************************************************************
      IF(NODE.EQ.3.OR.NODE.EQ.6) INTY=1
      DO 302 I=1,INTX     !INTX-INTEGRATION POINTS IN X DIRECTION
      DO 302 J=1,INTY    !INTY-INTEGRATION POINTS IN Y DIRECTION
      IF(NODE.EQ.4.OR.NODE.EQ.8) THEN
      CALL GAUSS_INTEGRATION(INTX,INTY,I,J,X,Y,WXY)
      CALL SHAPE_QUADRANGLE_8(NODE,X,Y,IEL,VN,VDN)
C                             !VN-SHAPE FUNCTION
C                             !VDN-LOCAL DERIVATE OF VN
      ENDIF
      IF(NODE.EQ.9) THEN
      CALL GAUSS_INTEGRATION(INTX,INTY,I,J,X,Y,WXY)
      CALL SHAPE_QUADRANGLE_9(NODE,X,Y,IEL,VN,VDN)
      ENDIF
      IF(NODE.EQ.3.OR.NODE.EQ.6) THEN
      CALL HAMMER_INTEGRATION(INTX,I,X,Y,Z,WXY)
      CALL SHAPE_TRIANGLE(NODE,X,Y,Z,IEL,VN,VDN)
      ENDIF
C******
C*********************FORM THE JACOBI MAXTRIX[J] AT INTEGRATION POINTS
C******
      CALL ELEMENT_JACOBI(MND,VXY,VDN,SJ,VD0)
C                                      !SJ=|J|
C                                      !VD0-GLOBAL DERIVATIVE OF VN
C                                      !VDN-LOCAL DERIVATE OF VN
      IF(SJ.LE.0.0) THEN
      WRITE(*,99)IE,I,J,SJ
99    FORMAT(/3X,'***SJ.LE.0.0 IN ELEMENT=',I4,3X,'INTX=',I2,
     $                   3X,'INTY=',I2,3X,'SJ=',E10.4)
      ENDIF
C***
C**********FORM THE[B]MATRIX AT INTEGRATION POINTS
C***
      CALL ELEMENT_VB(MPROB,MND,VXY,VN,VD0,VB,SR)
                                             !VB-ELEMENT STRAIN MATRIX
                                             !SR-PARAMETER OF INTEGRATION
C***
C*****************FORM THE ELEMENT STRESS MATRIX[S]=[D]*[B]
C***
      DO 20 II=1,NFSTR
      DO 20 JJ=1,MND2
      VSG(II,JJ)=0.0
C                         !VSG-ELEMENT STRESS MATRIX AT INTEGRATION POINT
      DO 20 KK=1,NFSTR
      VSG(II,JJ)=VSG(II,JJ)+VD(II,KK)*VB(KK,JJ)
20    CONTINUE
C***
C*****************FORM ELEMENT STIFFNESS MATRIX:[K]=[B]*[S]
C***
      DO 22 II=1,MND2
      DO 22 JJ=1,MND2
      DO 22 KK=1,NFSTR
      VKK(II,JJ)=VKK(II,JJ)+VB(KK,II)*VSG(KK,JJ)*WXY*SJ*SR
C                          !VKK-ELEMENT STIFFNESS MATRIX
22    CONTINUE                            
C***
C******************FORM THE ELEMENT MASS MATRIX[M]
C***
      DENS=VMAE(3)
      DO 25 II=1,MND
      DO 25 JJ=1,MND
      DO 25 KK=1,NF
      II1=NF*(II-1)+KK
      JJ1=NF*(JJ-1)+KK
      VMM(II1,JJ1)=VMM(II1,JJ1)+DENS*VN(II)*VN(JJ)*WXY*SJ*SR
C                                              !VMM-ELEMENT MASS MATRIX
25    CONTINUE
302   CONTINUE
      RETURN
      END
C================================SUB3-1-1==================================
      SUBROUTINE GAUSS_INTEGRATION(INTX,INTY,I,J,X,Y,WXY)
C*******************************************************************
C                 GET THE INFORMATION OF GAUSS INTEGRATION POINT
C**************************************************************8
      IMPLICIT   REAL*8(A-H,O-Z)
      DIMENSION   GXY(3,3),WG(3,3)
C***
C********************GAUSS INTEGRATION CONSTANTS FOR 1,2 AND 3POINTS
C***
      GXY(1,1)=0.0
      WG(1,1)=2.0
      GXY(1,2)=-0.577350269189626
      GXY(2,2)=0.577350269189626
      WG(1,2)=1.0
      WG(2,2)=1.0
      GXY(1,3)=-0.774596669241483
      GXY(2,3)=0.0
      GXY(3,3)=0.774596669241483
      WG(1,3)=0.555555555555556
      WG(2,3)=0.888888888888889
      WG(3,3)=0.555555555555556
C***
C*******************GET PARAMETERS OF INTEGRATION POINT
C***
      X=GXY(I,INTX)
      Y=GXY(J,INTY)
      WXY=WG(I,INTX)*WG(J,INTY)
      RETURN
      END
C=======================SUB:3-1-2======================================================
      SUBROUTINE SHAPE_QUADRANGLE_8(NODE,X,Y,IEL,VN,VDN)
C***********************************************************************
C         COMPUTE SHAPE FUNCTION OF QUADRANGLE ELEMENT AT INTEGRATION POINT
C*******************************************************************************
      IMPLICIT   REAL*8(A-H,O-Z)
      DIMENSION   IEL(NODE),VN(9),VDN(3,9)
C***
C******************CLEAR THE MEMORY FOR SHAPE FUNCTION AND ITS DERIVATIVES
C***
      DO 10 I=1,9
      VN(I)=0.0           !VN-SHAPE FUNCTION
      DO 10 J=1,3
      VDN(J,I)=0.0      !VDN-LOCAL DERIVATIVE OF VN
10    CONTINUE
C***
C******************SET FUNCTION VALUE FOR QUADRANDGLE ELEMENT OF 4 NODES
C***
      VN(1)=(1-X)*(1-Y)/4
      VN(2)=(1+X)*(1-Y)/4
      VN(3)=(1+X)*(1+Y)/4
      VN(4)=(1-X)*(1+Y)/4
      VDN(1,1)=-(1-Y)/4
      VDN(1,2)=(1-Y)/4
      VDN(1,3)=(1+Y)/4
      VDN(1,4)=-(1+Y)/4
      VDN(2,1)=-(1-X)/4
      VDN(2,2)=-(1+X)/4
      VDN(2,3)=(1+X)/4
      VDN(2,4)=(1-X)/4
C***
C*****************SET FUNCTION VALUE FOR QUADRANGLE ELEMENT OF THE 5-8 NODES
C***
      IF(NODE.EQ.8) THEN
      VN(5)=(1-X*X)*(1-Y)/2
      VN(6)=(1-Y*Y )*(1+X)/2
      VN(7)=(1-X*X)*(1+Y)/2
      VN(8)=(1-Y*Y)*(1-X)/2
      VDN(1,5)=(-2*X)*(1-Y)/2
      VDN(1,6)=(1-Y*Y)*(+1)/2
      VDN(1,7)=(-2*X)*(1+Y)/2
      VDN(1,8)=(1-Y*Y)*(-1)/2
      VDN(2,5)=(1-X*X)*(-1)/2
      VDN(2,6)=(-2*Y)*(1+X)/2
      VDN(2,7)=(1-X*X)*(+1)/2
      VDN(2,8)=(-2*Y)*(1-X)/2
      DO 30 I=1,4
      IF(IEL(4+I).EQ.0) VN(4+I)=0.0           !IEL-NODE CODE IN A ELEMENT
      IF(IEL(4+I).EQ.0) VDN(1,4+I)=0.0
      IF(IEL(4+I).EQ.0) VDN(2,4+I)=0.0
30    CONTINUE
      VN(1)=VN(1)-(VN(5)+VN(8))/2
      VN(2)=VN(2)-(VN(5)+VN(6))/2 
    VN(3)=VN(3)-(VN(6)+VN(7))/2
      VN(4)=VN(4)-(VN(7)+VN(8))/2                     
      DO 40 I=1,2
      VDN(I,1)=VDN(I,1)-(VDN(I,5)+VDN(I,8))/2
      VDN(I,2)=VDN(I,2)-(VDN(I,5)+VDN(I,6))/2
      VDN(I,3)=VDN(I,3)-(VDN(I,6)+VDN(I,7))/2
      VDN(I,4)=VDN(I,4)-(VDN(I,7)+VDN(I,8))/2
40    CONTINUE
      ENDIF
      RETURN
      END
C=========================SUB:3-1-3=================================================================
      SUBROUTINE SHAPE_QUADRANGLE_9(NODE,X,Y,IEL,VN,VDN)
C***********************************************************************************
C      COMPUTE SHAPE FUNCTION OF QUADRNGLE ELEMENT ON INTEGRATION POINT
C************************************************************************************
      IMPLICIT  REAL*8(A-H,O-Z)
      DIMENSION IEL(NODE),VN(9),VDN(3,9)
      DIMENSION IX(9),IY(9),VL0X(3),VL1X(3),VL0Y(3),VL1Y(3)
      DATA IX/1,-1,-1,1,0,-1,0,1,0/
      DATA IY/1,1,-1,-1,1,0,-1,0,0/
C***
C*****************CLEAR THE MEMORY FOR SHAPE FUNCTION AND ITS DERIVATIVES
C***
      DO 10 I=1,9
      VN(I)=0.0                 !VN-SHAPE FUNCTION  
      DO 10 J=1,3
      VDN(J,I)=0.0             !VDN-LOCAL DERIVATIVE OF VN
10    CONTINUE
C***
C******************SET THE VALUE OF LAGRANGEFUNCTION AND ITS DERIVATIVES
C***
      VL0X(1)=0.5*X*(X-1)
      VL0X(2)=1-X*X
      VL0X(3)=0.5*X*(X+1)
      VL1X(1)=X-0.5
      VL1X(2)=-2*X
      VL1X(3)=X+0.5
      VL0Y(1)=0.5*Y*(Y-1)
      VL0Y(2)=1-Y*Y
      VL0Y(3)=0.5*Y*(Y+1)
      VL1Y(1)=Y-0.5
      VL1Y(2)=-2*Y
      VL1Y(3)=Y+0.5
C***
C*******************SET SHAPE FUNCTION FOR QUADRANGLE ELEMENT OF 9 NODES
C***
      DO 20 I=1,9
      IIX=IX(I)+2
      IIY=IY(I)+2
      VN(I)=VL0X(IIX)*VL0Y(IIY)
      VDN(1,I)=VL1X(IIX)*VL0Y(IIY)
      VDN(2,I)=VL0X(IIX)*VL1Y(IIY)
20    CONTINUE
      IF(IEL(1).EQ.0) RETURN
      END
C============================SUB:3-1-4================================================
      SUBROUTINE HAMMER_INTEGRATION(INTX,I,X,Y,Z,WX)
C***************************************************************************
C                  GET THE INFORMATION OF HAMMER INTEGRATION POINT
C***************************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  HXY(3,5),WH(5),INDEX1(7)
C                                !HXY-COORDINATE OF HAMMER INTEGRATION POINT
C                                !WH-WEIGHT OF HAMMER INTEGRATION POINT
      INDEX1(1)=1
      INDEX1(3)=2
      INDEX1(4)=3
      INDEX1(7)=4
C***
C***********INTEGRATION CONSTANTS OF ONE POINT
C***
      HXY(1,1)=1.0/3.0
      HXY(2,1)=1.0/3.0
      HXY(3,1)=1.0/3.0
      WH(1)=1.0
C***
C***********INTEGRATION CONSTANTS OF 3 POINTS
C***
      HXY(1,2)=2.0/3.0
      HXY(2,2)=1.0/6.0
      HXY(3,2)=1.0/6.0
      WH(2)=1.0/3.0
C***
C**********INTEGRATION CONSTANTS OF 4 POINTS
C***
      HXY(1,3)=0.6
      HXY(2,3)=0.2
      HXY(3,3)=0.2
      WH(3)=25.0/48.0
C***
C**********INTEGRATION CONSTANTS OF 7 POINTS 
C***
      A1=0.0597158717
      B1=0.4701420641
      A2=0.7974269853
      B2=0.1012865073
      HXY(1,4)=A1
      HXY(2,4)=B1
      HXY(3,4)=B1
      WH(4)=0.1323941527
      HXY(1,5)=A2
      HXY(2,5)=B2
      HXY(3,5)=B2
      WH(5)=0.1259391805
C***
C**********GET PARAMETERS OF INTEGRATION POINTS
C***
      X=HXY(I+0-(I-1)/3*3,INDEX1(INTX))
      Y=HXY(I+2-(I+1)/3*3,INDEX1(INTX))
      Z=HXY(I+1-(I+0)/3*3,INDEX1(INTX))
      WX=WH(INDEX1(INTX))/2.0
CCCCC
      IF(INTX.EQ.7.AND.I.GE.4) THEN
      J=I-3
      X=HXY(J+0-(J-1)/3*3,5)
      Y=HXY(J+2-(J+1)/3*3,5)
      Z=HXY(J+1-(J+0)/3*3,5)
      WX=WH(5)/2.0
      ENDIF
CCCCC
      IF(INTX.EQ.4.AND.I.EQ.4) THEN
      X=HXY(I+0-(I-1)/3*3,1)
      Y=HXY(I+1-(I+0)/3*3,1)
      Z=HXY(I+2-(I+1)/3*3,1)
      WX=-27.0/48.0/2.0
      ENDIF
CCCCC
      IF(INTX.EQ.7.AND.I.EQ.7) THEN
      X=HXY(1,1)
      Y=HXY(2,1)
      Z=HXY(3,1)
      WX=0.9/8.0
      ENDIF
      RETURN
      END
C==========================SUB:3-1-5==========================================
      SUBROUTINE SHAPE_TRIANGLE(NODE,X,Y,Z,IEL,VN,VDN)
C*********************************************************************
C                COMPUTE SHAPE FUNCTION OF TRIANGLE ELEMENT AT INTEGRATION POINT
C*********************************************************************************
      IMPLICIT   REAL*8(A-H,O-Z)
      DIMENSION   IEL(NODE),VN(9),VDN(3,9)
C***
C******************CLEAR THE MEMORY FOR SHAPE FUNCTION AND ITS DERIVATIVES
C***
      DO 10 I=1,9
      VN(I)=0.0               !VN-SHAPE FUNCTION
      DO 10 J=1,3
      VDN(J,I)=0.0            !VDN-LOCAL DERIVATIVE OF VN
10    CONTINUE
C***
C******************SET THE FUNCTION VALUE FOR 3-NODE TRIANGLE ELEMENT
C***
      VN(1)=X
      VN(2)=Y
      VN(3)=Z
      VDN(1,1)=1.0
      VDN(1,2)=0.0
      VDN(1,3)=0.0
      VDN(2,1)=0.0
      VDN(2,2)=1.0
      VDN(2,3)=0.0
      VDN(3,1)=0.0
      VDN(3,2)=0.0
      VDN(3,3)=1.0
C***
C******************SET THE FUNCTION VALUE FOR TRIANGLE ELEMENT OF 4-6NODES
C***
      IF(NODE.EQ.6) THEN
      VN(4)=4*X*Y
      VN(5)=4*Y*Z
      VN(6)=4*Z*X
      VDN(1,4)=4*Y
      VDN(1,5)=0.0
      VDN(1,6)=4*Z
      VDN(2,4)=4*X
      VDN(2,5)=4*Z
      VDN(2,6)=0.0
      VDN(3,4)=0.0
      VDN(3,5)=4*Y
      VDN(3,6)=4*X
      DO 20 I=1,3
      IF(IEL(3+I).EQ.0) VN(3+I)=0.0
      IF(IEL(3+I).EQ.0) VDN(1,3+I)=0.0
      IF(IEL(3+I).EQ.0) VDN(2,3+I)=0.0
      IF(IEL(3+I).EQ.0) VDN(3,3+I)=0.0
20    CONTINUE
      VN(1)=VN(1)-(VN(4)+VN(6))/2    
      VN(2)=VN(2)-(VN(4)+VN(5))/2           
      VN(3)=VN(3)-(VN(5)+VN(6))/2
      DO 30 I=1,3
      VDN(I,1)=VDN(I,1)-(VDN(I,4)+VDN(I,6))/2    
      VDN(I,2)=VDN(I,2)-(VDN(I,4)+VDN(I,5))/2      
      VDN(I,3)=VDN(I,3)-(VDN(I,5)+VDN(I,6))/2
30    CONTINUE
      ENDIF
      DO 40 I=1,NODE
      VDN(1,I)=VDN(1,I)-VDN(3,I)
      VDN(2,I)=VDN(2,I)-VDN(3,I)
40    CONTINUE
      RETURN
      END
C==============================SUB:3-1-6=================================================
      SUBROUTINE ELEMENT_VD(MPROB,VD,VMAE)
C************************************************************************************
C**************FORM ELEMENT ELASTIC MATRIX[D] ACCORDING TO TYPE OF PROBLEM
C***********************************************************************************
      IMPLICIT  REAL*8(A-H,O-Z)
      DIMENSION   VD(5,5),VMAE(4)
      E=VMAE(1)                                        !E-ELASTIC MODULUS
      V=VMAE(2)                                        !V-POSSION'S RATIO
C***
C**************CLEAR MEMORY FOR THE MATRIX:[D]
C***
      DO 30 I=1,5
      DO 30 J=1,5
      VD(I,J)=0.0                                         !VD-ELASTICMATRIX
30    CONTINUE
C***
C**************COMPUTE[D] MATRIX FOR PLANE STRESS OR PLANE STRAIN PROBLEM
C***
      IF(MPROB.EQ.1.OR.MPROB.EQ.2) THEN
      IF(MPROB.EQ.2)  E=E/(1-V*V)
      IF(MPROB.EQ.2)  V=V/(1-V)
      D0=E/(1-V*V)
      VD(1,1)=D0                                       !MPROB=1--PLANE STRESS
      VD(2,2)=D0                                           
      VD(3,3)=D0*(1-V)/2                               !MPROB=2--PLANE STRAIN
      VD(1,2)=D0*V
      VD(2,1)=D0*V
      ENDIF
C***
C***************COMPUTE [D]MATRIX FOR AXISYMMETRIC PROBLEM
C***
      IF(MPROB.EQ.3) THEN
      D0=E*(1-V)/(1+V)/(1-2*V)
      VD(1,1)=D0
      VD(2,2)=D0
      VD(3,3)=D0*(1-2*V)/2/(1-V) 
      VD(4,4)=D0
      VD(2,1)=D0*V/(1-V)                             !MPROB=3-AXISYMMETRIC
      VD(1,2)=D0*V/(1-V)
      VD(4,1)=D0*V/(1-V)
      VD(1,4)=D0*V/(1-V)
      VD(4,2)=D0*V/(1-V)
      VD(2,4)=D0*V/(1-V)
      ENDIF
C***
C****************COMPUTE[D]MATRIX FOR MINDLIN PLATE PROBLEM
C***
      IF(MPROB.EQ.4) THEN
      TH=VMAE(4)
      D0=E*TH*TH*TH/12/(1-V*V)
      VD(1,1)=D0
      VD(2,2)=D0                           !MPROB=4-MINDLIN PLATE
      VD(3,3)=D0*(1-V)/2
      VD(1,2)=D0*V
      VD(2,1)=D0*V
      VD(4,4)=E/2/(1+V)*TH/(6.0/5.0)
      VD(5,5)=E/2/(1+V)*TH/(6.0/5.0)
      ENDIF
      RETURN
      END
C============= ==================SUB:3-1-7===============================================
      SUBROUTINE ELEMENT_JACOBI(MND,VXY,VDN,SJ,VD0)
C*************************************************************************
C                GET THE DETERMINANT OF  JACOBI MATRIX AND CARTESIAN DERIVATIVES
C********************************************************************************
      IMPLICIT  REAL*8(A-H,O-Z)
      DIMENSION    VXY(MND,2),VDN(3,9),VD0(3,9)
      DIMENSION    VJJ(2,2),VJ1(2,2)
C******
C*******************FORM JACOBI MATRIX
C******
      DO 11 II=1,2
      DO 11 JJ=1,2
      VJJ(II,JJ)=0.0                !VJJ-JACOBI MATRIX[J]
      DO 11 KK=1,MND
      VJJ(II,JJ)=VJJ(II,JJ)+VDN(II,KK)*VXY(KK,JJ)
11    CONTINUE
C***
C********************FORM THE DETERMINANT OF JACOBI MATRIX
C***
      SJ=VJJ(1,1)*VJJ(2,2)-VJJ(1,2)*VJJ(2,1)         !SJ-|J|
C***
C*********************COMPUTE THE INVERSE OF JACOBI MATRIX
C***
      VJ1(1,1)=+VJJ(2,2)/SJ
      VJ1(1,2)=-VJJ(1,2)/SJ           !VJ1-INVERSE OF [J]
      VJ1(2,1)=-VJJ(2,1)/SJ
      VJ1(2,2)=+VJJ(1,1)/SJ
C***
C*********************COMPUTE THE CARTESIAN DERIVATIVES OF SHAPE FUNCTION
C***
      DO 51 II=1,2
      DO 51 JJ=1,9
      VD0(II,JJ)=0.0                 !VD0-GLOBAL DERIVATIVE OF VN
      DO 51 KK=1,2
      VD0(II,JJ)=VD0(II,JJ)+VJ1(II,KK)*VDN(KK,JJ)
51    CONTINUE
      RETURN
      END
C===========================SUB:3-1-8=======================================
      SUBROUTINE ELEMENT_VB(MPROB,MND,VXY,VN,VD0,VB,SR)
C**********************************************************************
C                    FORM ELEMENT STRAIN MATRIX[B]ACCODING TO TYPE OF PROBLEM
C***********************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  VXY(MND,2),VN(9),VD0(3,9),VB(5,27)
                                    !VXY-单元节点坐标
                                    !VN-形函数
                                    !VD0-整体形函数导数
                                    !VB-应变矩阵
C***
C**********************CLEAR MEMORY FOR THE MATRIX[B]
C***

      DO 30 JJ=1,5
      DO 30 II=1,27
      VB(JJ,II)=0.0
30    CONTINUE
C***
C***********************COMPUTE[B] MATRIX FOR PLANE STRESS OR PLANE STRAIN PROBLEM
C***
      IF(MPROB.EQ.1.OR.MPROB.EQ.2) THEN
      SR=1.0
      DO 40 II=1,9
      VB(1,(II-1)*2+1)=VD0(1,II)                      !MPROB=1-PLANE STRESS
      VB(2,(II-1)*2+2)=VD0(2,II) 
      VB(3,(II-1)*2+1)=VD0(2,II)                      !MPROB=2-PLANE STRAIN
      VB(3,(II-1)*2+2)=VD0(1,II) 
40    CONTINUE
      ENDIF
C***
C************************COMPUTE[B]MATRIX FOR AXISYMMETRIX PROBLEM
C***
      IF(MPROB.EQ.3) THEN
    SR=0.0
    DO 45 II=1,MND
    SR=SR+VXY(II,1)*VN(II)
45  CONTINUE                                       !MPROB=3-AXISYMMETRIC
    DO 50 II=1,9
      VB(1,(II-1)*2+1)=VD0(1,II)                    
      VB(2,(II-1)*2+2)=VD0(2,II) 
      VB(3,(II-1)*2+1)=VD0(2,II)                      
      VB(3,(II-1)*2+2)=VD0(1,II) 
    IF(ABS(SR).LT.1E-20) VB(4,(II-1)*2+1)=0.0
    IF(ABS(SR).GE.1E-20) VB(4,(II-1)*2+1)=VN(II)/SR
50  CONTINUE
    SR=SR*2*3.14159265
    ENDIF
C***
C***********************COMPUTE[B]MATRIC FOR THE MINDLIN PLATE PROBLEM
C***
      IF(MPROB.EQ.4) THEN
      SR=1.0                                          !MPROB=4-MINDLIN PLATE
      DO 70 II=1,9
      VB(1,(II-1)*3+1)=-VD0(1,II)
      VB(2,(II-1)*3+2)=-VD0(2,II)
      VB(3,(II-1)*3+1)=-VD0(2,II)
      VB(3,(II-1)*3+2)=-VD0(1,II)
      VB(4,(II-1)*3+1)=-VN(II)
      VB(4,(II-1)*3+3)=-VD0(1,II)
      VB(5,(II-1)*3+2)=-VN(II)
      VB(5,(II-1)*3+3)=-VD0(2,II)
70    CONTINUE
      ENDIF
      RETURN
      END
C=============================SUB:4============================================================
C********************************************************************************************
C       SOLVE STATIC PROBLEM TO GET THE NODE DISPLACEMENTS
C         CALL  SUBROUTINES :DECOMPOS AND BACKSUBS
C***************************************************************************************************
      SUBROUTINE STATIC_SOLVE(GKK,GP,GU)
      IMPLICIT   REAL*8(A-H,O-Z)
      COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI
      COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA
      COMMON/COM3/MND2,NUMPT2
      COMMON/ELEM/NODE,INTX,INTY
      DIMENSION GKK(NUMPT2,MBAND),GP(NUMPT2),GU(NUMPT2)
      WRITE(*,101)
101   FORMAT(/5X,'## SOLVE EQUATION [GKK]{GU}={GP}##')
C***
C*******************TRIANGLE DECOMPOSITION OF THE MATRIX:[GKK]
C***

C***
C********************SUBSTITUTE BACK TO GET THE VECTOR {GU}
C***
      CALL GXIAO(NUMPT2,MBAND,GKK,GP)
      DO 30 I=1,NUMPT2
      GU(I)=GP(I)
30    CONTINUE
C***************************************************************
C                      OUTPUT THE NODAL DISPLACEMENT
C****************************************************************
      WRITE(*,103)
103   FORMAT(/6X,'#OUTPUT NODAL DISPLACEMENT TO FILE<OUT_DIS>#')
      IF(MPROB.EQ.4) THEN
      WRITE(14,40)
      ELSE
      WRITE(14,39)
      ENDIF
39    FORMAT(10X,'NODAL DISPLACEMENTS'/2X,'NO.OF NODES',10X,
     $                   'X-',15X,'Y-')
40    FORMAT(10X,'NODAL DISPLANCEMENTS'/2X,'NO.OF NODES',5X,
     $                  'THETA-X',7X,'THETA-Y',10X,'W-Z')
      DO 41 I=1,NUMPT
      IF(NF.EQ.2) WRITE(14,42) I,(GU(NF*(I-1)+J),J=1,NF) 
      IF(NF.EQ.3) WRITE(14,43) I,(GU(NF*(I-1)+J),J=1,NF) 
41    CONTINUE
42    FORMAT(2X,I5,4X,2E18.8)
43    FORMAT(2X,I5,4X,3E16.8)
      CLOSE(14)
      RETURN
      END
C=============================SUB:4-1==============================
      SUBROUTINE GXIAO(NUMPT2,MBAND,ARRAY1,ARRAY2)
    IMPLICIT REAL*8(A-H,O-Z)
    DIMENSION ARRAY1(NUMPT2,MBAND), ARRAY2(NUMPT2)
    DO 20 K=1,NUMPT2-1
    IF(NUMPT2.GT.K+MBAND-1)THEN
    IM=K+MBAND-1
    ELSE
    IM=NUMPT2
    END IF
    DO 20 I=K+1,IM
    IL=I-K+1
    C=ARRAY1(K,IL)/ARRAY1(K,1)
    DO 10 J=1,MBAND-IL+1
    M=J+I-K
10  ARRAY1(I,J)=ARRAY1(I,J)-C*ARRAY1(K,M)
    ARRAY2(I)=ARRAY2(I)-C*ARRAY2(K)
20  CONTINUE     
    ARRAY2(NUMPT2)=ARRAY2(NUMPT2)/ARRAY1(NUMPT2,1)
    DO 40 I=NUMPT2-1,1,-1
    IF(MBAND.GT.NUMPT2-I+1)THEN
    JO=NUMPT2-I+1
    ELSE
    JO=MBAND
    END IF
    DO 30 J=2,JO
    IH=J+I-1
30  ARRAY2(I)=ARRAY2(I)-ARRAY1(I,J)*ARRAY2(IH)
40  ARRAY2(I)=ARRAY2(I)/ARRAY1(I,1)
50  CONTINUE
    RETURN
    END
C=============================SUB:4-2==============================
      SUBROUTINE DECOMPOS(NUMPT2,MBAND,ARRAY)
C******************************************************************************************
C                   TRIANGLE DECOMPOSITION OF MATRIX
C*********************************************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
    DIMENSION ARRAY(NUMPT2,MBAND)
    DO 14 M=1,NUMPT2-1
    MM=M+MBAND-1
    IF(MM.GT.NUMPT2)MM=NUMPT2
    DO 13 I=M+1,MM
    MI=I+MBAND-1
    IF(MI.GT.NUMPT2)MI=NUMPT2
    DO 10 J=I,MI
    ARRAY(I,J-I+1)=ARRAY(I,J-I+1)-ARRAY(M,I-M+1)*ARRAY(M,J-M+1)
     $ /ARRAY(M,1)
10  CONTINUE
13    CONTINUE
14    CONTINUE
      RETURN
    END
C=============================SUB:4-3===========================================================
      SUBROUTINE BACKSUBS(NUMPT2,MBAND,ARRAY1,ARRAY2)
C***********************************************************************************************
C                    BACKSUBSTITUTION OF TRIANGLE DECOMPOSITION
C**********************************************************************************************
    IMPLICIT REAL*8(A-H,O-Z)
    DIMENSION ARRAY1(NUMPT2,MBAND),ARRAY2(NUMPT2)
    DO 11 M=1,NUMPT2-1
    MM=M+MBAND-1
    IF(MM.GT.NUMPT2)MM=NUMPT2
    DO 11 I=M+1,MM
    MI=I+MBAND-1
    IF(MI.GT.NUMPT2)MI=NUMPT2
    DO 11 J=1,MI
    ARRAY2(I)=ARRAY2(I)-ARRAY1(M,I-M+1)*ARRAY2(M)/ARRAY1(M,1)
11  CONTINUE
    DO 20 I=NUMPT2,1,-1
    MI=I+MBAND-1
    IF(MI.GT.NUMPT2)MI=NUMPT2
    DO 21 J=I+1,MI
    ARRAY2(I)=(ARRAY2(I)-ARRAY1(I,J-I+1)*ARRAY2(J))/ARRAY1(I,1)
21  CONTINUE
    ARRAY2(I)=ARRAY2(I)/ARRAY1(I,1)
20  CONTINUE
    RETURN
    END
C===========================SUB:5==============================================================
      SUBROUTINE STRESS(VMATI,IELEM,VCOOD,GU,VXY,IEL,VU,SSS,
     $                        VSS,IADD,SSN)
C*******************************************************************************
C                    COMPUTE STRESSES AT THE INTEGRATION POINTS AND NODES 
C                     CALL SUBROUTINE STRESS_MATRIX
C***********************************************************************************
      IMPLICIT  REAL*8(A-H,O-Z)
      COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI
      COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA
      COMMON/COM3/MND2,NUMPT2
      COMMON/ELEM/NODE,INTX,INTY
      DIMENSION IELEM(NUMEL,MND+4),VCOOD(NUMPT,2),GU(NUMPT2),
     $             VXY(MND,2),IEL(MND),VU(MND2),VSS(4,NFSTR),
     $             SSS(9,NFSTR),IADD(NUMPT),VMATI(NMATI,4),
     $             SSN(NUMPT,NFSTR),VSNN(9,NFSTR),VMAE(4)
C
      DO 20 I=1,NUMPT
      IADD(I)=0
      DO 20 J=1,NFSTR
      SSN(I,J)=0.0
20    CONTINUE
      WRITE(*,101)
101   FORMAT(/6X,'# OUTPUT ELEMENT STRESS TO FILE<OUT_STR>#')
      IF(MTYPE.EQ.0.OR.MTYPE.EQ.2) THEN
      IF(MPROB.EQ.1.OR.MPROB.EQ.2) WRITE(15,21)
      IF(MPROB.EQ.3) WRITE(15,22)
      IF(MPROB.EQ.4) WRITE(15,23)
      ENDIF
21    FORMAT(10X,'STRESS ON EACH INTEGRATION POINT'/'ELEM-NO.',1X,
     $          'INTEG-NO.',5X,'SINGMA-X',7X,'SIGMA-Y',9X,'TAO-XY')
22    FORMAT(10X,'STRESS ON EACH INTEGRATION POINT'/'ELEM-NO.',1X,
     $          'INTEG-NO.',3X,'SIGMA-X',8X,'SIGMA-Y',8X,'TAO-XY',10X,
     $          'ROTATION')
23    FORMAT(10X,'STRESS ON EACH INTEGRATION POINT'/'ELEM-NO.',1X,
     $          'INTEG.',5X,'M-X',10X,'M-Y',9X,'M-XY',7X,'TAO-XZ',
     $          7X,'TAO-YZ')
C***
C*********FORM ALL THE ELEMENT INFORMATION FROM THE INPUT DATA
      DO 320 IE=1,NUMEL
    DO 30 I=I,MND
      IEL(I)=IELEM(IE,I+4)
      DO 30 J=1,2
      VXY(I,J)=0.0
      IF(IEL(I).GT.0) VXY(I,J)=VCOOD(IEL(I),J)
30    CONTINUE
      NODE=IELEM(IE,1)            !NODE-NUMBER OF NODEES IN THE ELEMENT
      INTX=IELEM(IE,3)            !INTX-INTEGERATION POINTS IN X-DIRECEION
      INTY=IELEM(IE,4)                !INTY-INTEGERATION POINTS IN Y-DIRECEION
      IF(NODE.EQ.3.OR.NODE.EQ.6) INTY=1
      INTXY=INTX*INTY
      IMATI=IELEM(IE,2)
      DO 31 J=1,4
31    VMAE(J)=VMATI(IMATI,J)
      DO 40 I=1,MND
      DO 40 J=1,NF
      IF (IEL(I).GT.0) VU(NF*(I-1)+J)=GU(NF*(IELEM(IE,I+4)-1)+J)
40  CONTINUE
C***
C******************COMPUTE STRESS AT THE INTEGRATION POINTS
C***
      CALL STRESS_MATRIX(IE,VXY,IEL,VMAE,VU,SSS,VSS)
C***
C******************OUTPUT STRESS AT THE INTEGRATION POINT
C***
      IF(MTYPE.EQ.0.OR.MTYPE.EQ.2)THEN
      DO 70 I=1,INTXY
      IF(MPROB.EQ.1.OR.MPROB.EQ.2)
     $               WRITE(15,71) IE,I,(SSS(I,J),J=1,NFSTR)
      IF(MPROB.EQ.3) WRITE(15,72) IE,I,(SSS(I,J),J=1,NFSTR)
      IF(MPROB.EQ.4) WRITE(15,73) IE,I,(SSS(I,J),J=1,NFSTR)
70    CONTINUE
      ENDIF
71    FORMAT(1X,I5,1X,I5,1X,3(2X,E15.6))
72    FORMAT(1X,I5,1X,I5,1X,4(1X,E15.6))
73    FORMAT(1X,I5,1X,I4,1X,5E13.6)
C***
C*****************COMPUTE NODAL STRESS FROM STRESS MATRIX<VSN>BY ELEMENT SMOOTHING
C***
      DO 83 J=1,NFSTR
      IF(NODE.EQ.3.OR.NODE.EQ.6) THEN
      A1=5.0/3.0
      B1=1.0/3.0
      VSNN(1,J)=VSS(1,J)*A1-VSS(2,J)*B1-VSS(3,J)*B1
      VSNN(2,J)=-VSS(1,J)*B1+VSS(2,J)*A1-VSS(3,J)*B1
      VSNN(3,J)=-VSS(1,J)*B1-VSS(2,J)*B1+VSS(3,J)*A1
      VSNN(4,J)=(VSNN(1,J)+VSNN(2,J))/2.0
      VSNN(5,J)=(VSNN(2,J)+VSNN(3,J))/2.0
      VSNN(6,J)=(VSNN(3,J)+VSNN(1,J))/2.0
      ELSE
      A=1.8660254
      B=-0.5
      C=0.1339746
      VSNN(1,J)=VSS(1,J)*A+VSS(2,J)*B+VSS(3,J)*C+VSS(4,J)*B
      VSNN(2,J)=VSS(1,J)*B+VSS(2,J)*A+VSS(3,J)*B+VSS(4,J)*C
      VSNN(3,J)=VSS(1,J)*C+VSS(2,J)*B+VSS(3,J)*A+VSS(4,J)*B
      VSNN(4,J)=VSS(1,J)*B+VSS(2,J)*C+VSS(3,J)*B+VSS(4,J)*A
      VSNN(5,J)=(VSNN(1,J)+VSNN(2,J))/2
      VSNN(6,J)=(VSNN(2,J)+VSNN(3,J))/2
      VSNN(7,J)=(VSNN(3,J)+VSNN(4,J))/2
      VSNN(8,J)=(VSNN(4,J)+VSNN(1,J))/2
      VSNN(9,J)=(VSNN(5,J)+VSNN(7,J))/2
      ENDIF
83    CONTINUE
C***
C*****************COMPUTE THE AVERAGE VALUE OF NODAL STRESSES
C***
    DO 89 I=1,MND
    IEL(I)=IELEM(IE,I+4)
    IH=IEL(I)
      IF(IH.GT.0) THEN
      IADD(IH)=IADD(IH)+1
      DO 90 J=1,NFSTR
      SSN(IH,J)=(SSN(IH,J)*(IADD(IH)-1)+VSNN(I,J))/IADD(IH)                                                              
90    CONTINUE
      ENDIF
89    CONTINUE
320   CONTINUE
C***
C******************OUTPUT STRESS AT THE NODES
C***
      IF(MPROB.EQ.1.OR.MPROB.EQ.2) WRITE(15,91)
      IF(MPROB.EQ.3) WRITE(15,92)
      IF(MPROB.EQ.4) WRITE(15.93)
91    FORMAT(/10X,'STRESSES AT EACH NODE'/2X,'NODE-NO.',8X,
     $            'SIGMA-X',10X,'SIGMA-Y',10X,'TAO-XY',10X,'ROTATION')
92    FORMAT(/10X,'STRESSES AT EACH NODE'/2X,'NODE-NO.',8X,
     $            'SIGMA-X',10X,'SIGMA-Y',10X,'TAO-XY',10X,'ROTATION')
93    FORMAT(/10X,'STRESSES AT EACH NODES'/2X,'NODE-NO.',5X,'M-X',
     $            10X,'M-Y',10X,'M-XY',8X,'TAO-XZ',9X,'TAO-YZ')
      DO 95 I=1,NUMPT 
      IF(MPROB.EQ.1.OR.MPROB.EQ.2)
     $               WRITE(15,96) I,(SSN(I,J),J=1,NFSTR)
      IF(MPROB.EQ.3) WRITE(15,97) I,(SSN(I,J),J=1,NFSTR)
      IF(MPROB.EQ.4) WRITE(15,98) I,(SSN(I,J),J=1,NFSTR)
95    CONTINUE
96    FORMAT(2X,I5,2X,3(2X,E15.6))
97    FORMAT(2X,I5,2X,4(2X,E15.6))
98    FORMAT(1X,I5,1X,5E14.6)
      CLOSE(15)
      RETURN
      END
C==============================SUB:5-1=============================================
      SUBROUTINE STRESS_MATRIX(IE,VXY,IEL,VMAE,VU,SSS,VSS,GU)
C***********************************************************************************
C               FORM THE ELEMENT STRESS MATRIXES:[VSG]AND[VSN]
C*******************************************************************************
      IMPLICIT  REAL*8(A-H,O-Z)
      COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI
      COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA
      COMMON/COM3/MND2,NUMPT2
      COMMON/ELEM/NODE,INTX,INTY
      DIMENSION VXY(MND,2),IEL(MND),VMAE(4),VU(MND2),VSS(4,NFSTR),
     $         SSS(9,NFSTR),VSG(NFSTR,MND2),VSN(NFSTR,MND2),GU(NUMPT2)
      DIMENSION VN(9),VDN(3,9),VD0(3,9),VD(5,5),VB(5,27)
C***
C*******************FORM THE[D]MATRIX ACCORDING TO TYPE OF THE PROBLEM
C***
      CALL ELEMENT_VD(MPROB,VD,VMAE)     !VD-ELASTIC MATRIX
C************************************************************************************
C                      LOOP OVER THE NUMERICAL INTEGRATION POINTS(INTX,INTY)
C*************************************************************************************
      IF(NODE.EQ.3.OR.NODE.EQ.6) INTY=1
      DO 32 I=1,INTX
      DO 32 J=1,INTY
      IF(NODE.EQ.4.OR.NODE.EQ.8) THEN
      CALL GAUSS_INTEGRATION(INTX,INTY,I,J,X,Y,WXY)
      CALL SHAPE_QUADRANGLE_8(NODE,X,Y,IEL,VN,VDN)
      ENDIF                        !VN-SHAPE FUNCTION
      IF(NODE.EQ.9) THEN           !VU-单元节点位移
      CALL GAUSS_INTEGRATION (INTX,INTY,I,J,X,Y,WXY)
      CALL SHAPE_QUADRANGLE_9(NODE,X,Y,IEL,VN,VDN)
      ENDIF                         !VDN-LOCAL DERIVATIVE OF[VN]
      IF(NODE.EQ.3.OR.NODE.EQ.6)THEN
    CALL HAMMER_INTEGRATION(INTX,I,X,Y,Z,WXY)
    CALL SHAPE_TRIANGLE(NODE,X,Y,Z,IEL,VN,VDN)
    ENDIF
    CALL ELEMENT_JACOBI(MND,VXY,VDN,SJ,VD0)    !SJ-|J|
                                    !VD0-GLOBAL DERIVATIVE OF[VN]
      CALL ELEMENT_VB(MPROB,MND,VXY,VN,VD0,VB,SR)
                                    !VB-ELEMENT STRAIN MATRIX
                                    !SR-PARAMETER OF INTEGRATION
C***
C********************FORM ELEMENT STRESS MATRIX[VSG]AT INTEGRATION POINTS
C***                                   ([VSG]=[D]*[B])
      DO 20 II=1,NFSTR
      DO 20 JJ=1,MND2
      VSG(II,JJ)=0.0
      DO 22 KK=1,NFSTR
      VSG(II,JJ)=VSG(II,JJ)+VD(II,KK)*VB(KK,JJ)
22    CONTINUE
20  CONTINUE 
      k=(I-1)*INTY+J
      DO 30 II=1,NFSTR
    Z=0.0
      IF(INTX.EQ.3.AND.INTY.EQ.1) VSS(K,II)=0.0
      IF(INTX.EQ.2.AND.INTY.EQ.2) VSS(K,II)=0.0
    DO 30 JJ=1,MND2
      Z=Z+VSG(II,JJ)*VU(JJ)
    SSS(K,II)=Z
      IF(INTX.EQ.3.AND.INTY.EQ.1) VSS(K,II)=SSS(K,II)
      IF(INTX.EQ.2.AND.INTY.EQ.2) VSS(K,II)=SSS(K,II)
30  CONTINUE
32    CONTINUE
C*********************************************************************************************
C             FORM ELEMENT STESS MATRIX[VSN] AT OPTIMAL STRESS POINTS
C                              ( [VSN]=[D]*[B])
C**************************************************************************************************
      IF(NODE.EQ.3.OR.NODE.EQ.6) THEN
      INTR=3
      INTS=1
      IF(INTX.EQ.3.AND.INTY.EQ.1) GOTO 46
      ELSE
      INTR=2
      INTS=2
      IF(INTX.EQ.2.AND.INTY.EQ.2) GOTO 46
      ENDIF
      DO 45 I=1,INTR
      DO 45 J=1,INTS
      IF(NODE.EQ.4.OR.NODE.EQ.8) THEN
      CALL GAUSS_INTEGRATION(INTR,INTS,I,J,X,Y,WXY)    
      CALL SHAPE_QUADRANGLE_8(NODE,X,Y,IEL,VN,VDN)
      ENDIF
      IF(NODE.EQ.9) THEN
      CALL GAUSS_INTEGRATION(INTR,INTS,I,J,X,Y,WXY)    
      CALL SHAPE_QUADRANGLE_9(NODE,X,Y,IEL,VN,VDN)   
      ENDIF
      IF(NODE.EQ.3.OR.NODE.EQ.6) THEN
      CALL HAMMER_INTEGRATION(INTR,I,X,Y,Z,WXY)    
      CALL SHAPE_TRIANGLE(NODE,X,Y,Z,IEL,VN,VDN)
      ENDIF
      CALL ELEMENT_JACOBI(MND,VXY,VDN,SJ,VD0)
      CALL ELEMENT_VB(MPROB,MND,VXY,VN,VD0,VB,SR)
      DO 42 II=1,NFSTR
      DO 42 JJ=1,MND2
      VSN(II,JJ)=0.0
      DO 42 KK=1,NFSTR
      VSN(II,JJ)=VSN(II,JJ)+VD(II,KK)*VB(KK,JJ)
42    CONTINUE
      K=(I-1)*INTS+J
      DO 43 II=1,NFSTR
      VSS(K,II)=0.0
      DO 43 JJ=1,MND2
      VSS(K,II)=VSS(K,II)+VSN(II,JJ)*VU(JJ)      
43    CONTINUE
45    CONTINUE
46    CONTINUE
      RETURN
      END
C==============================SUB:6============================================
      SUBROUTINE DYNAM(GMM,GKK,GP,DAMP,U0,V0,A,AW,B)
C*******************************************************************************
C                 SOLVE THE DYNAMIC PROBLEM BY CENTRAL DIFFERENCE METHOD
C                         OR NEWMARK METHOD
C                 CALL SUBROUTINES:DECOMPOS,BACKSUBS,CENTER OR NEWMARK
C********************************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI
      COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA
      COMMON/COM3/MND2,NUMPT2
      COMMON/DYN/HUV,OMEGA,CC1,CC2,TT,DT,ALPHA,DELTA
      DIMENSION GMM(NUMPT2,MBAND),GKK(NUMPT2,MBAND),GP(NUMPT2),
     $       A(NUMPT2),U0(NUMPT2),V0(NUMPT2),AW(NUMPT2,MBAND),
     $       B(NUMPT2),DAMP(NUMPT2,MBAND)
C***
C*********************COMPUTE DAMPLING MATRIX(DAMP)
C***
      DO 50 I=1,NUMPT2
      DO 50 J=1,MBAND
      DAMP(I,J)=CC1*GMM(I,J)+CC2*GKK(I,J)
50    CONTINUE
C***
C**********************FORM INITIAL DISPLACEMENT(U0) AND VELOCITY(V0)
C***
      DO 60 I=1,NUMPT2
      U0(I)=0.0
      V0(I)=0.0
60    CONTINUE
      IF(HUV.EQ.1.OR.HUV.EQ.3) THEN
      READ(5,*)
      READ(5,*)(U0(I),I=1,NUMPT2)
      ENDIF
      IF(HUV.EQ.2.OR.HUV.EQ.3) THEN
      READ(5,*)
      READ(5,*)(V0(I),I=1,NUMPT2)
      ENDIF
      WRITE(6,*)'INITIAL DISPLACEMENTS-U0:'
      WRITE(6,62)(U0(I),I=1,NUMPT2)
      WRITE(6,*)'INITIAL VELOCITIES-V0:'
      WRITE(6,62)(V0(I),I=1,NUMPT2)
62    FORMAT(2X,4E15.6)                   
C**************************************************************************************
C                    COMPUTE INITIAL ACCELERATION-A(NUMPT2)
C***********************************************************************************
      DO 71 I=1,NUMPT2
71    A(I)=GP(I)
      DO 70 I=1,NUMPT2
      IK=I-MBAND+1
      IF(IK.LT.1) IK=1
      DO 70 J=IK,I
      A(I)=A(I)-GKK(J,I-J+1)*U0(J)-DAMP(I,J-I+1)*V0(J)
70    CONTINUE
      DO 72 I=1,NUMPT2
    IK=I+MBAND-1
    IF(IK.GT.NUMPT2) IK=NUMPT2
    DO 72 J=I+1,IK
    A(I)=A(I)-GKK(I,J-I+1)*U0(J)-DAMP(I,J-1)*V0(J)
72  CONTINUE
      DO 73 I=1,NUMPT2
      DO 73 J=1,MBAND
      AW(I,J)=GMM(I,J)
73    CONTINUE
      CALL DECOMPOS(NUMPT2,MBAND,AW)
      CALL BACKSUBS(NUMPT2,MBAND,AW,A)   !^^^^^^^^^
C***
C*******************BY USING GENTRAL DIFFERENCE METHOD WHEN MSOLV=2
C***
      IF(MSOLV.EQ.2) CALL CENTER(GMM,GKK,DAMP,GP,U0,V0,A,AW,B)
      RETURN
      END
C===========================SUB:6-1=====================================================
      SUBROUTINE CENTER(GMM,GKK,DAMP,GP,U0,V0,A,AW,B)
C********************************************************************************
C                      SOLVE DYNAMIC RESPONSE BY CENTRAL DIFFERENCE METHOD
C                          CALL SUBROUTINES:DECOMPOS,BACKSUBS
C*********************************************************************************
      IMPLICIT  REAL*8(A-H,O-Z)
      COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI
      COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA
      COMMON/COM3/MND2,NUMPT2
      COMMON/DYN/HUV,OMEGA,CC1,CC2,TT,DT,ALPHA,DELTA
      DIMENSION GMM(NUMPT2,MBAND),GKK(NUMPT2,MBAND),
     $          DAMP(NUMPT2,MBAND),GP(NUMPT2),U0(NUMPT2),
     $          V0(NUMPT2),AW(NUMPT2,MBAND),A(NUMPT2),B(NUMPT2)
      CLOSE(30,STATUS='DELETE')
      OPEN(30,FILE='OUT_CEN',STATUS='UNKNOWN')
      WRITE(30,*)'DYNAMIC RESPONSE RESULT BY CENTRAL DIFFERENCE'
      WRITE(30,*)'TOTAL TIME=',TT
      WRITE(*,101)
101   FORMAT(/6X,'#SOLVE DYNAMIC RESPONSE BY CENTRAL DIFFERENCE')
      WRITE(*,102)
102   FORMAT(/6X,'%OUTPUT DYNAMIC DISPLACEMENT IN FILE<OUT_PUT>')
C***
C***********************INITIAL COMPUTATIONS
C***
      C0=1./(DT*DT)
      C1=0.5/DT
      C2=2.*C0
      C3=1./C2
      DO 30 I=1,NUMPT2
      A(I)=U0(I)-DT*V0(I)+C3*A(I)
30    CONTINUE
      DO 40 I=1,NUMPT2
      B(I)=GP(I)
      DO 40 J=1,MBAND
      AW(I,J)=GMM(I,J)
40    GMM(I,J)=C0*GMM(I,J)+C1*DAMP(I,J)
C***
C************************COMPUTATIONS FOR EACH TIME STEP
C***
      CALL DECOMPOS(NUMPT2,MBAND,GMM)
      DO 300 Y=0,TT,DT
      IF(OMEGA.GT.0.0) AP=SIN(OMEGA*Y)
      IF(OMEGA.LT.0.0) AP=COS(OMEGA*Y)
      DO 41 I=1,NUMPT2
      GP(I)=B(I)
      IF(OMEGA.NE.0.0) GP(I)=GP(I)*AP
41    CONTINUE
      DO 50 I=1,NUMPT2
      IK=I-MBAND+1
      IF(IK.LT.1) IK=1
      DO 50 J=IK,I
      GP(I)=GP(I)-(GKK(J,I-J+1)-C2*AW(J,I-J+1))*U0(J)-
     $          (C0*AW(J,I-J+1)-C1*DAMP(J,I-J+1))*A(J)
50    CONTINUE
      DO 60 I=1,NUMPT2
      IK=I+MBAND-1
      IF(IK.GT.NUMPT2) IK=NUMPT2
      DO 60 J=I+1,IK
      GP(I)=GP(I)-(GKK(I,J-I+1)-C2*AW(I,J-I+1))*U0(J)-
     $ (C0*AW(I,J-I+1)-C1*DAMP(I,J-I+1))*A(J)
60    CONTINUE
      CALL BACKSUBS(NUMPT2,MBAND,GMM,GP)
      NSTEP=Y/DT+1
      WRITE(30,61) NSTEP,DT,Y+DT
      IF(MPROB.EQ.4) THEN
      WRITE(30,73)
      WRITE(30,74)(GP(I),I=1,NUMPT2)
      ELSE
      WRITE(30,71)
      WRITE(30,72)(GP(I),I=1,NUMPT2)
      ENDIF
      DO 70 I=1,NUMPT2
      A(I)=U0(I)
      U0(I)=GP(I)
70    CONTINUE
300   CONTINUE
61    FORMAT(2X,'NO.OF STEP=',I5,3X,'STEP LENGTH=',E15.6,3X,
     $           'AT TIME=',E15.6)
71    FORMAT(2X,'DISPLACEMENT:'/2(13X,'X-',13X,'Y-'))
72    FORMAT(4X,4E16.8)
73    FORMAT(2X,'DISPLACEMENT:'/11X,'THETA-X',11X,'THETA-Y',
     $            12X,'W-Z')
74    FORMAT(6X,3E16.8)
      CLOSE(30)
      RETURN
      END
C=============================SUB:6-2===============================================
      SUBROUTINE NEWMARK(GMM,GKK,DAMP,GP,U0,V0,A,AW,B)
C***************************************************************************
C                   SOLVE DYNAMIC RESPONSE BY NEWMARK METHOD
C                        CALL SUBROUTINES:DECOMPOS,BACKSUBS
C***************************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI
      COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA
      COMMON/COM3/MND2,NUMPT2
      COMMON/DYN/HUV,OMEGA,CC1,CC2,TT,DT,ALPHA,DELTA
      DIMENSION   GMM(NUMPT2,MBAND),GKK(NUMPT2,MBAND),
     $            AW(NUMPT2,MBAND),DAMP(NUMPT2,MBAND),GP(NUMPT2),
     $            U0(NUMPT2),V0(NUMPT2),A(NUMPT2),B(NUMPT2)
       CLOSE(31,STATUS='DELETE')
      OPEN(31,FILE='OUT_NMK',STATUS='UNKNOWN')
      WRITE(31,*)'DYNAMIC RESPONSE RESULT BY NEWMARK METHOD'
      WRITE(31,*)'TOTAL TIME=',TT
      WRITE(*,101)
101   FORMAT(/6X,'# SOLVE DYNAMIC RESPONSE BY NEWMARK METHOD##')
      WRITE(*,102)
102   FORMAT(/6X,'% OUTPUT DYNAMIC DISPLACEMENT IN FILE<OUT_NMK>')
C***
C*****************************INITIAL COMPUTATIONS
C***
      C0=1.0/(ALPHA*DT*DT)
      C1=DELTA/(ALPHA*DT)
      C2=1.0/(ALPHA*DT)
      C3=0.5/ALPHA-1.0
      C4=DELTA/ALPHA-1.0
      C5=DT/2.0*(DELTA/ALPHA-2.0)
      C6=DT*(1.0-DELTA)
      C7=DELTA*DT
C************************************************************************************
C               COMPUTE  K'=K+C0*M+C1*C
C************************************************************************************
      DO 40 I=1,NUMPT2
      B(I)=GP(I)
      DO 40 J=1,MBAND
      AW(I,J)=GKK(I,J)
40    GKK(I,J)=GKK(I,J)+C0*GMM(I,J)+C1*DAMP(I,J)
C***********************************************************************************
C               TRIANGLE DECOMPOSITION OF THE MATRIX:[GKK]
C***********************************************************************************
      CALL DECOMPOS(NUMPT2,MBAND,GKK)
C***
C*************************COMPUTATIONS FOR EACH TIME STEP
C***
      DO 300 Y=0,TT,DT
      IF(OMEGA.GT.0.0) AP=SIN(OMEGA*Y)
      IF(OMEGA.LT.0.0) AP=COS(OMEGA*Y)
      DO 41 I=1,NUMPT2
      GP(I)=B(I)
      IF(OMEGA.NE.0.0) GP(I)=GP(I)*AP
41    CONTINUE
      DO 50 I=1,NUMPT2
      IK=I-MBAND+1
      IF(IK.LT.1) IK=1
      DO 50 J=IK,I
      GP(I)=GP(I)+GMM(J,I-J+1)*(C0*U0(J)+C2*V0(J)+C3*A(J))+
     $            DAMP(J,I-J+1)*(C1*U0(J)+C4*V0(J)+C5*A(J))
50    CONTINUE
      DO 60 I=1,NUMPT2
      IK=I+MBAND-1
      IF(IK.GT.NUMPT2) IK=NUMPT2
      DO 60 J=I+1,IK
      GP(I)=GP(I)+GMM(J,I-J+1)*(C0*U0(J)+C2*V0(J)+C3*A(J))+
     $             DAMP(J,I-J+1)*(C1*U0(J)+C4*V0(J)+C5*A(J))
60    CONTINUE
      CALL BACKSUBS(NUMPT2,MBAND,GKK,GP)
      NSTEP=Y/DT+1
      WRITE(31,61) NSTEP,DT,Y+DT
      IF(MPROB.EQ.4) THEN
      WRITE(31,73)
      WRITE(31,74) (GP(I),I=1,NUMPT2)
      ELSE
      WRITE(31,71)
      WRITE(31,72) (GP(I),I=1,NUMPT2)
      ENDIF
      DO 70 I=1,NUMPT2
      A(I)=C0*(GP(I)-U0(I))-C2*V0(I)-C3*A(I)
      V0(I)=V0(I)+C6*A(I)+C7*A(I)
      U0(I)=GP(I)
70    CONTINUE
300   CONTINUE
61    FORMAT(2X,'NO.OF STEP=',I5,3X,'STEP LENGTH=',E15.6,3X,
     $                'AT TIME=',E15.6)
71    FORMAT(2X,'DIAPLACEMENT:'/2(13X,'X-',13X,'Y-'))
72    FORMAT(4X,4E16.8)
73    FORMAT(2X,'DISPLACEMENT:'/11X'THETA-X',11X,'THETA-Y',
     $                12X,'W-Z')
74    FORMAT(6X,3E16.8)
      CLOSE(31)
      RETURN
      END
C===================================SUB:7==============================================
      SUBROUTINE EIGEN(GMM,GKK,AA,BB,GM,GK,V,W1,W2)
C*************************************************************************************
C                    SOLVE EIGENVALUE PROBLEM BY INVERSE OR SUBSPACE METHODS
C                       CALL SUBROUTINE INVERSE OR SUBROUTINE SUBSPACE
C***************************************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI
      COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA
      COMMON/COM3/MND2,NUMPT2
      DIMENSION GMM(NUMPT2,MBAND),GKK(NUMPT2,MBAND),AA(NUMPT2,NVA),
     $             BB(NUMPT2,NVA),GM(NVA,NVA),GK(NVA,NVA),V(NVA,NVA),
     $             W1(NVA),W2(NVA)
C***
C***********************BY USING INVERSE METHOD WHEN MSOLV=4
C***
      IF(MSOLV.EQ.4) CALL INVERSE(GMM,GKK,AA,BB)
C***
C***********************BY USING SUBSPACE METHOD WHEN MSOLV=5
C***
      IF(MSOLV.EQ.5) CALL SUBSPACE(GMM,GKK,AA,BB,GM,GK,V,W1,W2)
      RETURN
      END
C===============================SUB:7-1==================================================
      SUBROUTINE INVERSE(GMM,GKK,AA,BB)
C******************************************************************************************
C               SOLVE EIGENVALUE BY INVERSE ITERATION METHOD
C                  CALL SUBROUTINES:DECOMPOS,BACKSUBS
C****************************************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI
      COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA
      COMMON/COM3/MND2,NUMPT2
      DIMENSION GKK(NUMPT2,MBAND),GMM(NUMPT2,MBAND),AA(NUMPT2,NVA),
     $                 BB(NUMPT2)
      CLOSE(32,STATUS='DELETE')
      OPEN(32,FILE='OUT_VERS',STATUS='UNKNOWN')
      WRITE(32,*)'DYNAMIC CHARACTER RESUITS BY INVERSE METHOD'
      WRITE(*,101)
101   FORMAT(/6X,'# SOLVE EIGENVALUE BY INVERSE METHOD#')
      WRITE(*,102)
102   FORMAT(/6X,'% OUTPUT EIGENVALUE AND MODE IN FILE<OUT_VERS>')
C***
C**************************FORM INITIAL ITERATION VECTORS A(I,J)
C***
      DO 5 I=1,NUMPT2
      DO 5 J=1,NVA
      K=I+J
    AA(I,J)=RAN(K)
5     CONTINUE
C***
C**************************TRIANGLE DECOMPOSITION OF STIFFNESS MATRIX
C***
      CALL DECOMPOS(NUMPT2,MBAND,GKK)
      DO 100 II=1,NVA
      W1=-1.0
      GK=0.0
      GM=0.0
C***
C*************************COMPUTE  Y=MX
C***
333   CONTINUE
      DO 8 I=1,NUMPT2
8     BB(I)=0.0
      DO 10 I=1,NUMPT2
      IK=I-MBAND+1
      IF(IK.LT.1) IK=1
      DO 10 J=IK,I
      BB(I)=BB(I)+GMM(J,I-J+1)*AA(J,II)
10    CONTINUE
      DO 12 I=1,NUMPT2
      IK=I+MBAND-1
      IF(IK.GT.NUMPT2) IK=NUMPT2
      DO 12 J=I+1,IK
      BB(I)=BB(I)+GMM(I,J-I+1)*AA(J,II)
12    CONTINUE
C***
C***********************SOLVE KX=Y
C***
111   CONTINUE
      NSTEP=NSTEP+1
      IF(NSTEP.GT.500) THEN
      WRITE(*,*)'NO.=',II,'STEP=',NSTEP
      WRITE(32,*)'NO.=',II,'STEP=',NSTEP
      RETURN
      ENDIF
      DO 20 I=1,NUMPT2
20    AA(I,II)=BB(I)
      CALL BACKSUBS(NUMPT2,MBAND,GKK,AA(I,II))
C***
C**********************COMPUTE W**2=K'/M'
C***
      DO 24 I=1,NUMPT2
      GK=GK+AA(I,II)*BB(I)
      BB(I)=0.0
24    CONTINUE
      DO 25 I=1,NUMPT2
      IK=I-MBAND+1
      IF(IK.LT.1) IK=1
      DO 25 J=IK,I
      BB(I)=BB(I)+GMM(I,J-I+1)*AA(J,II)
25    CONTINUE
      DO 26 I=NUMPT2,1,-1
      IK=I+MBAND-1
      IF(IK.GT.NUMPT2) IK=NUMPT2
      DO 26 J=I+1,IK
      BB(I)=BB(I)+GMM(I,J-I+1)*AA(J,II)
26    CONTINUE
      DO 27 I=1,NUMPT2
27    GM=GM+AA(I,II)*BB(I)
      W2=GK/GM
      W1=ABS((W1-W2)/W2)
      IF(W1.GT.1.0E-6) THEN
      W1=W2
C***
C*********************GRAM-SCHMIDT ORTHOGONIGATION
C***
      IF(II.EQ.1) GOTO 222
      IF(II.NE.1) THEN
      DO 30 JJ=1,II-1
      S=0.0
      DO 40 I=1,NUMPT2
      IK=I-MBAND+1
      IF(IK.LT.1) IK=1
      DO 40 J=IK,I
      S=S+AA(I,JJ)*GMM(J,I-J+1)*AA(J,II)
40    CONTINUE
      DO 50 I=1,NUMPT2
      IK=I+MBAND-1          
      IF(IK.GT.NUMPT2) IK=NUMPT2
      DO 50 J=I+1,IK
      S=S+AA(I,JJ)*GMM(I,J-I+1)*AA(J,II)
50    CONTINUE
      DO 51 K=1,NUMPT2
      AA(K,II)=AA(K,II)-AA(K,JJ)*S
51    CONTINUE
30    CONTINUE
      S=0.0
      DO 52 I=1,NUMPT2
      IK=I-MBAND+1
      IF(IK.LT.1) IK=1
      DO 52 J=IK,I
      S=S+AA(I,II)*GMM(J,I-J+1)*AA(J,II)
52    CONTINUE
      DO 54 I=1,NUMPT2
      IK=I+MBAND-1
      IF(IK.GT.NUMPT2) IK=NUMPT2
      DO 54 J=I+1,IK
      S=S+AA(I,II)*GMM(I,J-I+1)*AA(J,II)
54  CONTINUE
      DO 55 K=1,NUMPT2
      AA(K,II)=AA(K,II)/S
55    CONTINUE
      GK=0.0
      GM=0.0
      GOTO 333
      ENDIF
222   DO 80 I=1,NUMPT2
      W2=SQRT(GM)
      BB(I)=BB(I)/W2
80    CONTINUE
      GK=0.0
      GM=0.0
      GOTO 111
      ELSE
C***
C***************************COMPUTE FREQUENCY,PERIOD,EIGENVECTOR
C***
      WW=SQRT(W2)
      PD=2*3.1415926/WW
      WRITE(32,93) II,NSTEP,WW,PD
      DO 90 I=1,NUMPT2
      W2=SQRT(GM)
      AA(I,II)=AA(I,II)/W2
90    CONTINUE
      IF(MPROB.EQ.4) THEN
      WRITE(32,92)
      WRITE(32,95)(AA(I,II),I=1,NUMPT2)
      ELSE
      WRITE(32,91)
      WRITE(32,94)(AA(I,II),I=1,NUMPT2)
      ENDIF
      ENDIF
      NSTEP=0
100   CONTINUE
91    FORMAT(2X,'VIBRATION MODE:'/2(13X,'X-',13X,'Y-'))
92    FORMAT(2X,'VIBRATION MODE:'/11X,'THETA-X',11X,'THETA-Y',
     $            12X,'W-Z')
93    FORMAT('***',2X,'NO.OF EIGENVALUE=',I5,4X,'ITERATION TIMES=',
     $            I5/5X,'FREQUENCE=',F16.4,4X,'PERION='E16.6)
94    FORMAT(2X,4E16.8)
95    FORMAT(4X,3E16.8)
      CLOSE(32)
      RETURN
      END
C=======================================SUB:7-2=====================================
      SUBROUTINE SUBSPACE(GMM,GKK,AA,BB,GM,GK,V,W1,W2)
C********************************************************************************
C                SOLVE EIGENVALUE BY SUBSPACE ITERATION METHOD
C                CALL SUBROUTINES:DECOMPOS,BACKSUBS,JOCOBI,ARRANGE
C**************************************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI
      COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA
      COMMON/COM3/MND2,NUMPT2
      DIMENSION GMM(NUMPT2,MBAND),GKK(NUMPT2,MBAND),AA(NUMPT2,NVA),
     $            BB(NUMPT2,NVA),GM(NVA,NVA),GK(NVA,NVA),V(NVA,NVA),
     $            W2(NVA),GMK(NUMPT2),W1(NVA)
      CLOSE(33,STATUS='DELETE')
      OPEN(33,FILE='OUT_SUBS',STATUS='UNKNOWN')
      WRITE(*,101)
      WRITE(33,*)'RESULTS OF EIGENVALUES BY SUBSPACE METHOD:'
101   FORMAT(/6X,'#SOLVE EIGENVALUE BY SUBSPACE METHOD#')
      WRITE(*,102)
102   FORMAT(/6X,'% OUTPUT EIGENVALUE AND MODE INTO FILE<OUT_SUBS>')
C***
C***************************FORM INITIAL ITERATION VECTORS
C***
      DO 5 I=1,NUMPT2
      GMK(I)=GMM(I,1)/GKK(I,1)
      AA(I,1)=1.0
5     CONTINUE
      DO 7 I=1,NUMPT2
      DO 7 J=2,NVA
7     AA(I,J)=0.0
      DO 8 J=2,NVA
      QQ=0.0
      DO 9 I=1,NUMPT2
      IF(GMK(I).GT.QQ) THEN
      QQ=GMK(I)
      K=I
      ENDIF
9     CONTINUE
      GMK(K)=0.0
      AA(K,J)=1.0
8     CONTINUE
C***
C************************TRIANGLE DECOMPOSTITION OF STIFFNESS MATRIX
C***
      CALL DECOMPOS(NUMPT2,MBAND,GKK)
C***
C************************COMPUTE Y=MX
C***
      DO 10 K=1,NVA
      DO 11 I=1,NUMPT2
      IK=I-MBAND+1
      IF(IK.LT.1) IK=1
      DO 11 J=IK,I
      BB(I,K)=BB(I,K)+GMM(J,I-J+1)*AA(J,K)
11    CONTINUE
      DO 12 I=1,NUMPT2
      IK=I+MBAND-1
      IF(IK.GT.NUMPT2) IK=NUMPT2
      DO 12 J=I+1,IK
      BB(I,K)=BB(I,K)+GMM(I,J-I+1)*AA(J,K)
12    CONTINUE
C***
C**********************SOLVING EQUATION  KX=Y
C***
      DO 20 I=1,NUMPT2
20    AA(I,K)=BB(I,K)
10    CONTINUE
      NSTEP=0
111   CONTINUE
      NSTEP=NSTEP+1
      DO 30 J=1,NUMPT2
      DO 30 K=1,NVA
30    BB(J,K)=AA(J,K)
      DO 32 K=1,NVA
      CALL BACKSUBS(NUMPT2,MBAND,GKK,AA(1,K))
32    CONTINUE
C***
C***********************COMPUTE   K1=XT*Y
C***
      DO 40 I=1,NVA
      DO 40 J=1,NVA
      GK(I,J)=0.0
      GK(I,J)=GK(I,J)+AA(K,I)*BB(K,J)
40    CONTINUE
C***
C***********************COMPUTE Y1=M*X1
C***
      DO 45 K=1,NVA
      DO 50 I=1,NUMPT2
      IK=I-MBAND+1
      IF(IK.LT.1) IK=1
      DO 50 J=IK,I
      BB(I,K)=BB(I,K)+GMM(J,I-J+1)*AA(J,K)
50    CONTINUE
      DO 60 I=1,NUMPT2
      IK=I+MBAND-1
      IF(IK.GT.NUMPT2) IK=NUMPT2
      DO 60 J=I+1,IK
      BB(I,K)=BB(I,K)+GMM(J,I-J+1)*AA(J,K)
60    CONTINUE
45    CONTINUE
C***
C*******************COMPUTE M1=XT*Y1
C***
      DO 70 I=1,NVA
      DO 70 J=1,NVA
      GM(I,J)=0.0
      DO 70 K=1,NUMPT2
      GM(I,J)=GM(I,J)+AA(K,I)*BB(K,J)
70    CONTINUE
C***********************************************************************************
C          SOLVE GENERAL EIGENVALUE PROBLEM
C***********************************************************************************
      CALL JACOBI(GK,GM,V,W2,NVA)
      CALL ARRANGE(W2,V,NVA)
C***
C**********************CHECK IF THE CONVERGENT CONDITION IS SATISFIED
C***
      IF(NVA.EQ.NUMPT2) GOTO 222
      DO 80 J=1,NVA
      IF(ABS((W2(J)-W1(J))/W2(J)).GT.1.E-4) THEN
      DO 82 I=1,NUMPT2
      DO 82 K=1,NVA
      AA(I,K)=0.0
      DO 82 M=1,NVA
      AA(I,K)=AA(I,K)+BB(I,M)*V(M,K)
82    CONTINUE
      DO 85 K=1,NVA
85    W1(K)=W2(K)
      GOTO 111
      ENDIF
80    CONTINUE
222   CONTINUE
      DO 88 I=1,NUMPT2
      DO 88 J=1,NVA
      BB(I,J)=0.0
      DO 88 K=1,NVA
      BB(I,K)=BB(I,K)+AA(I,K)*V(K,J)
88    CONTINUE
      DO 90 J=1,NVA
      WW=SQRT(W2(J))
      PD=2*3.1415926/WW
      WRITE(33,91) J,NSTEP,WW,PD
      IF(MPROB.EQ.4) THEN
      WRITE(33,93)
      WRITE(33,95)(BB(I,J),I=1,NUMPT2)
      ELSE
      WRITE(33,92)
      WRITE(33,94)(BB(I,J),I=1,NUMPT2)
      ENDIF
90    CONTINUE
91    FORMAT('***',2X,'NO.OF EIGENVALUE=',I5,4X,'ITERATIN TIMES=',
     $       I5/5X,'FREQUENCY=',F16.4,4X,'PERIOD=',E16.6)
92    FORMAT(2X,'VIBRATION MODE:'/2(13X,'X-',13X,'Y-'))
93    FORMAT(2X,'VIBRATION MODE:'/11X,'THETA-X',11X,'THETA-Y',
     $       12X,'W-Z')
94    FORMAT(2X,4E16.8)
95    FORMAT(5X,3E16.8)
      CLOSE(33)
      RETURN
      END
C===================================SUB:7-2-1======================================== 
      SUBROUTINE JACOBI(GK,GM,V,EIGV,N)
C*******************************************************************************
C                  SOLVE EIGENVALUE BY JACOBI METHOD
C*********************************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION GK(N,N),GM(N,N),V(N,N),EIGV(N),D(N)
      RTOL=1.E-12
      NSMAX=15
      DO 10 I=1,N
      D(I)=0.0
      D(I)=GK(I,I)/GM(I,I)
10    EIGV(I)=D(I)
      DO 30 I=1,N
      DO 20 J=1,N
20    V(I,J)=0.0
30    V(I,I)=1.0
      NSWEEP=0
    NN=N-1
40    NSWEEP=NSWEEP+1
      EPS=(0.01**NSWEEP)**2
      DO 110 J=1,NN
      JJ=J+1
      DO 110 K=JJ,N
      EPTOLA=(GK(J,K)*GK(J,K))/(GK(J,J)*GK(K,K))
      EPTOLB=(GM(J,K)*GM(J,K))/(GM(J,J)*GM(K,K))
      IF((EPTOLA.LT.EPS).AND.(EPTOLB.LT.EPS)) GOTO 110
      AKK=GK(K,K)*GM(J,K)-GM(K,K)*GK(J,K)
      AJJ=GK(J,J)*GM(J,K)-GM(J,J)*GK(J,K)
      AB=GK(J,J)*GM(K,K)-GK(K,K)*GM(J,J)
      CHECK=(AB*AB+4.0*AKK*AJJ)/4.0
      IF(CHECK) 50,51,51
50    STOP 222
51    SQCH=SQRT(CHECK)
      D1=AB/2.+SQCH
      D2=AB/2.-SQCH
      DEN=D1
      IF(ABS(D2).GT.ABS(D1)) DEN=D2
      IF(DEN) 55,57,55
57    CA=0.0
      CG=-GK(J,K)/GK(K,K)
    GOTO 60
55  CA=AKK/DEN
    CG=-AJJ/DEN
60    IF(N-2) 61,90,61
61    JP1=J+1
      JM1=J-1
      KP1=K+1
      KM1=K-1
      IF(JM1-1) 63,62,62
62    DO 68 I=1,JM1
      AJ=GK(I,J)
      BJ=GM(I,J)
      AK=GK(I,K)
      BK=GM(I,K)
      GK(I,J)=AJ+CG*AK
      GM(I,J)=BJ+CG*BK
      GK(I,K)=AK+CA*AJ
68    GM(I,K)=BK+CA*BJ
63    IF(KP1-N) 64,64,66
64    DO 65 I=KP1,N
      AJ=GK(J,I)
      BJ=GM(J,I)
      AK=GK(K,I)
      BK=GM(K,I)
      GK(J,I)=AJ+CG*AK
      GM(J,I)=BJ+CG*BK
      GK(K,I)=AK+CA*AJ
65    GM(K,I)=BK+CA*BJ
66    IF(JP1-KM1) 70,70,90
70    DO 80 I=JP1,KM1
      AJ=GK(J,I)
      BJ=GM(J,I)
      AK=GK(I,K)
      BK=GM(I,K)
      GK(J,I)=AJ+CG*AK
      GM(J,I)=BJ+CG*BK
      GK(I,K)=AK+CA*AJ
80    GM(I,K)=BK+CA*BJ
90    AK=GK(K,K)
      BK=GM(K,K)
      GK(K,K)=AK+2.*CA*GK(J,K)+CA*CA*GK(J,J)
      GM(K,K)=BK+2.*CA*GM(J,K)+CA*CA*GM(J,J)
      GK(J,J)=GK(J,J)+2.*CG*GK(J,K)+CG*CG*AK
      GM(J,J)=GM(J,J)+2.*CG*GM(J,K)+CG*CG*BK
      GK(J,K)=0.0
      GM(J,K)=0.0
      DO 91 I=1,N
      XJ=V(I,J)
      XK=V(I,K)
      V(I,J)=XJ+CG*XK
91    V(I,K)=XK+CA*XJ
110   CONTINUE
      DO 92 I=1,N
      IF(GK(I,I).GT.0.0.AND.GM(I,I).GT.0.0) GOTO 92
      STOP 333
92    EIGV(I)=GK(I,I)/GM(I,I)
      DO 93 I=1,N
      TOL=RTOL*D(I)
      DIF=ABS(EIGV(I)-D(I))
      IF(DIF.GT.TOL) GOTO 97
93    CONTINUE
      EPS=TOL**2
      DO 94 J=1,NN
    JJ=J+1
    DO 94 K=JJ,N
      EPSA=(GK(J,K)*GK(J,K))/(GK(J,J)*GK(K,K))
      EPSB=(GM(J,K)*GK(J,K))/(GM(J,J)*GM(K,K))
      IF((EPSA.LT.EPS).AND.(EPSB.LT.EPS)) GOTO 94
      GOTO 97
94    CONTINUE
      DO 95 I=1,N
      DO 95 J=1,N
      GK(J,I)=GK(I,J)
95    GM(J,I)=GM(I,J)
      DO 96 J=1,N
      BB=SQRT(GM(J,J))
      DO 96 K=1,N
96    V(K,J)=V(K,J)/BB
      RETURN
97    DO 98 I=1,N
98    D(I)=EIGV(I)
      IF(NSWEEP.LT.NSMAX) GOTO 40
      GOTO 94
    RETURN
      END
C===========================SUB:7-2-2====================================================
      SUBROUTINE ARRANGE(W,V,N)
C*************************************************************************************
C                    ARRANGE EIGENVALUE INTO ORDER
C************************************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  W(N),V(N,N)                                                    
      DO 13 I=1,N-1
      K=I
      P=W(I)
      DO 11 J=I+1,N
      IF(W(J).LT.P) THEN
      K=J
      P=W(J)
    ENDIF
11    CONTINUE
      IF(K.NE.I) THEN
      W(K)=W(I)
      W(I)=P
      DO 12 J=1,N
      P=V(J,I)
      V(J,I)=V(J,K)
      V(J,K)=P
12    CONTINUE
      ENDIF
13    CONTINUE
      RETURN
      END