C                 2-D MICROPOLAR FINITE ELEMENT METHOD
C *                          SKYLINE=MICRO:2D                          *
C     NC        ....  NUMBER OF TOTAL CONSTRAINED NODE
C     ND        ....  NUMBER OF TOTAL DEGREE OF FREEDOM
C     NE        ....  NUMBER OF TOTAL ELEMENT
C     NN        ....  NUMBER OF TOTAL NODAL POINTS
C     TH        ....  THICKNESS OF THE ELEMENT
C
C                                 VECTOR
C     SDISP (ND)....  STRUCTURAL NODAL DISPLACEMENTS & ROTATIONS
C     EDISP ( 9)....  ELEMENT    NODAL DISPLACEMENTS & ROTATIONS
C     KSTRN (NC)....  CONSTRAINED NODE
C     KSTRT (NC)....  TYPE OF CONSTRAIN
C                        100  X-DISPLACEMENT  IS FIXED
C                        010  Y-DISPLACEMENT  IS FIXED
C                        001  Z-ROTATION      IS FIXED
C                        101  X-DIS. & Z-ROT. IS FIXED
C     INXY  (NE)....  LOCATION OF INTEREST IN ELEMENT FOR OUTPUT
C                        100  AT MIDPOINT BETWEEN I AND J NODES
C                        010  AT MIDPOINT BETWEEN J AND M NODES
C                        001  AT MIDPOINT BETWEEN M AND I NODES
C                        000  AT CENTROID OF TRIANGULAR ELEMNT
C
C                                 MATRIX
C     E   (NE,4)....  DISPLACEMENT-STRAIN
C     PHIJ(NE,4)....  ROTATION    -STRAIN
C     T   (NE,4)....  FORCE       -STRESS
C     CM  (NE,4)....  COUPLE      -STRESS
C     EK  ( 9,9)....  ELEMENT    STIFFNESS MATRIX
C     SK (ND,ND)....  STRUCTURAL STIFFNESS MATRIX
C     NNP (NE,3)....  NODE NUMBER FORMING THE ELEMENT
C
C                                 TENSOR
C     B0(NE,9,9)....  DISP-STRAIN DISPLACEMENT MATRIX B
C     B1(NE,9,9)....  ROTN-STRAIN ROTATION     MATRIX B
C
      IMPLICIT REAL*8(A-H,O-Z)
      INTEGER*4 ITIM
      INTEGER*2 TYPE,CODE
      COMMON/B01/B0(1400,9,9),B1(1400,9,9),INXY(1400),ID(1400)
      COMMON/XYZ/NNP(1400,3),X(800),Y(800)
      COMMON/CSTRN/KSTRN(80),KSTRT(80)
      COMMON/MAT/D0(4,4),D1(4,4),A11,A12,A22,A33,A34,A44,B22,B24,B44
      COMMON/AIN/NB,NC,ND,NE,NN,TH
      COMMON/FSK/A(200000),V(2400),EK(9,9),MAXA(2401),NEITR
C DIMENSION OF MAXA() MUST BE ONE GREATER THAN V().
C     CALL ASSIGN(5,'INPUT.IN')
C     CALL ASSIGN(6,'OUTPUT.OUT')
      OPEN (5,FILE='INPUT.IN')
      OPEN (6,FILE='OUTPUT.OUT')
C     SAVE
C
C I-------------------- INPUT & OUTPUT FILE OPEN ----------------------I
C$INSERT SYSCOM>A$KEYS.INS.FTN
C$INSERT SYSCOM>KEYS.INS.FTN
C$INSERT SYSCOM>ERRD.INS.FTN
C      CALL SRCH$$(K$DELE,'FEMO',INTS(4),INTS(2),TYPE,CODE)
C      CALL OPEN$A(A$WRIT+A$SAMF,'FEMO',INTS(4),INTS(2))
C      CALL SRCH$$(K$READ,'FEMI',INTS(4),INTS(1),TYPE,CODE)
C      CALL SRCH$$(K$WRIT,'FEMO',INTS(4),INTS(2),TYPE,CODE)
C      CT1=CTIM$A(ITIM)
C
C I---------------------------- DATA INPUT ----------------------------I
C NE? READ TOTAL ELEMENT NUMBER (NE=6)
      READ(5,100) NE
  100 FORMAT(I4)
C NP? READ ELEMENT CONNECTION INFORMATION
      DO 10 N=1,NE
      READ(5,200) NNP(N,1),NNP(N,2),NNP(N,3),INXY(N)
  200 FORMAT(4X,3I4,I6)
   10 CONTINUE
C NN? READ TOTAL NUMBER OF NODAL POINTS (NN=8)
      READ(5,100) NN
C XY? READ COORDINATE (X,Y) OF NODAL POINTS
      DO 20 N=1,NN
      READ(5,300) X(N),Y(N)
   20 CONTINUE
  300 FORMAT(4X,2F20.10)
C NC? READ NUMBER OF CONSTRAINED NODES (NC=3)
      READ(5,100) NC
C KS? READ CONSTRAINTS
      DO 30 N=1,NC
      READ(5,400) KSTRN(N),KSTRT(N)
   30 CONTINUE
  400 FORMAT(4X,2I4)
C TH? READ THICKNESS OF THE MATERIAL
      READ(5,500) TH
  500 FORMAT(E11.4)
C AB? READ MATERIAL PROPERTIES
      READ(5,600) A11,A12,A22,A33,A34
      READ(5,600) A44,B22,B24,B44
  600 FORMAT(5E14.7)
C
C I--------------------- TOTAL DEGREE OF FREEDOM ----------------------I
C EACH NODE HAS 3-DOF OF (UX,UY,PHIZ)
      ND=NN*3
C CALCULATION OF BAND-WIDTH OF STRUCTURAL STIFFNESS MATRIX SK(ND,ND)
      NB=0
      DO 40 IB=1,NE
      IMAX=MAX0(NNP(IB,1),NNP(IB,2),NNP(IB,3))
      IMIN=MIN0(NNP(IB,1),NNP(IB,2),NNP(IB,3))
      NBCHEK=(IMAX-IMIN+1)*3
      IF(NBCHEK.GT.NB) NB=NBCHEK
      IF(NB .GT. 120) GO TO 99
   40 CONTINUE
C
      SCALE=1.0D0
      DO 41 I=1,NN
      X(I)=X(I)*SCALE
   41 Y(I)=Y(I)*SCALE
C
C I--------------------------- DATA OUTPUT ----------------------------I
C PRINT THE HEAD OF OUTPUT
      WRITE(6,700)
  700 FORMAT(1H1)
      WRITE(6,800)
  800 FORMAT(3X,'*******************************************************
     ****************')
      WRITE(6,900)
  900 FORMAT(3X,'*',68X,'*')
      WRITE(6,1000)
 1000 FORMAT(3X,'*',13X,'2-D ORTHOTROPIC MICROPOLAR STRESS ANALYSIS',
     *       13X,'*')
      WRITE(6,1050)
 1050 FORMAT(3X,'*',28X,' SKYLINEMICRO',27X,'*')
      WRITE(6,900)
      WRITE(6,800)
C
C WRITE THE INFORMATIONS
      WRITE(6,1100)
 1100 FORMAT(/20X,'**** DISCRETIZATION NUMBER ****')
      WRITE(6,1200)
 1200 FORMAT(/13X,'ELEMT.#',3X,'NODES.#',3X,'CONSTR.#',3X,'THIKNES',
     *       3X,'BAND-WIDTH')
      WRITE(6,1300) NE,NN,NC,TH,NB
 1300 FORMAT(13X,I4,7X,I3,7X,I3,5X,E10.3,5X,I5)
C
C WRITE ORTHOTROPIC MATERIAL PROPERTIES
      WRITE(6,1400)
 1400 FORMAT(//22X,'**** MATERIAL PROPERTIES ****',
     *    //6X,'A11',9X,'A12',9X,'A22',9X,'A33',9X,'A34',9X,'A44')
      WRITE(6,1500) A11,A12,A22,A33,A34,A44
 1500 FORMAT(6E12.3)
      WRITE(6,1600)
 1600 FORMAT(/6X,'B22',9X,'B24',9X,'B44')
      WRITE(6,1500) B22,B24,B44
C
C WRITE ELEMENT CONNECTION INFORMATIONS
      WRITE(6,1700)
 1700 FORMAT(//21X,'**** ELEMENT-NODE CONNECTION ****')
      WRITE(6,1800)
 1800 FORMAT(/3X,'ELM NP1 NP2 NP3 IXY   ELM NP1 NP2 NP3 IXY   ELM NP1 NP
     *2 NP3 IXY')
      DO 45 I=1,NE
   45 ID(I)=I
      LINE=NE/3
      IRESID=NE-3*LINE
      DO 50 N=1,LINE
   50 WRITE(6,1900) (ID(3*(N-1)+I),NNP(3*(N-1)+I,1),NNP(3*(N-1)+I,2),
     * NNP(3*(N-1)+I,3),INXY(3*(N-1)+I),I=1,3)
 1900 FORMAT(3(2X,I4,1X,I3,1X,I3,1X,I3,1X,I3))
      IF(IRESID .EQ. 0) GO TO 56
      WRITE(6,1900) (ID(3*LINE+I),NNP(3*LINE+I,1),NNP(3*LINE+I,2),
     *NNP(3*LINE+I,3),INXY(3*LINE+I),I=1,IRESID)
C
C WRITE NODAL COORDINATES
   56 WRITE(6,2000)
 2000 FORMAT(//24X,'**** NODAL COORDINATE ****')
      WRITE(6,2100)
 2100 FORMAT(/1X,'NODE',5X,'X',8X,'Y',4X,'NODE',5X,'X',8X,'Y',
     * 4X,'NODE',5X,'X',8X,'Y')
      DO 55 I=1,NN
   55 ID(I)=I
      LINE=NN/3
      IRESID=NN-3*LINE
      DO 60 N=1,LINE
   60 WRITE(6,2200) (ID(3*(N-1)+I),X(3*(N-1)+I),Y(3*(N-1)+I),I=1,3)
 2200 FORMAT(3(2X,I3,1X,F8.3,1X,F8.3))
      IF(IRESID .EQ. 0) GO TO 57
      WRITE(6,2200) (ID(3*LINE+I),X(3*LINE+I),Y(3*LINE+I),I=1,IRESID)
C
C WRITE CONSTRAINTS
   57 WRITE(6,2300)
 2300 FORMAT(//27X,'**** CONSTRAINT ****')
      WRITE(6,2400)
 2400 FORMAT(/24X,'CNSTRND-NODE',2X,'CNSTRND-CODE')
      WRITE(6,2500) (KSTRN(N),KSTRT(N),N=1,NC)
 2500 FORMAT(28X,I3,12X,I3)
C
C I--------------------------- MAIN PROGRAM ---------------------------I
      WRITE(1,1)
    1 FORMAT(' ENTERING INTO STSTIF')
                         CALL STSTIF
      WRITE(1,2)
    2 FORMAT(' FINISHED STSTIF')
      WRITE(1,3)
    3 FORMAT(' ENTERING INTO LOADER')
                         CALL LOADER
      WRITE(1,4)
    4 FORMAT('FINISHED LOADER')
      WRITE(1,5)
    5 FORMAT(' ENTERING INTO COLSOL')
                         CALL COLSOL
      WRITE(1,6)
    6 FORMAT(' FINISHED COLSOL')
                         CALL STRESS
C
C     CT3=CTIM$A(ITIM)
C     C=CT3-CT1
C     WRITE(6,9999) C
C
   99 WRITE(6,999) NB,IB
  999 FORMAT('*****STOP NB=',I5,'  AT ELEMENT=',I5)
 9999 FORMAT(2X,'COMP. TIME T3=',F10.3)
C
C      CALL SRCH$$(K$CLOS,'FEMI',INTS(4),INTS(1),TYPE,CODE)
C      CALL SRCH$$(K$CLOS,'FEMO',INTS(4),INTS(2),TYPE,CODE)
C
      STOP
                                    END
C
C **********************************************************************
C *                                                                    *
                         SUBROUTINE STSTIF
C *                                                                    *
C **********************************************************************
C TO CALCULATE STRUCTURAL STIFFNESS MATRIX A(NB*ND)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/CSTRN/KSTRN(80),KSTRT(80)
      COMMON/XYZ/NNP(1400,3),X(800),Y(800)
      COMMON/AIN/NB,NC,ND,NE,NN,TH
      COMMON/FSK/A(200000),V(2400),EK(9,9),MAXA(2401),NEITR
      COMMON/MHH/MHT(2400)
C
      WRITE(1,1)
    1 FORMAT(' BEGIN STSTIF')
C FROM THE ARRAY A(NB*ND) FOR STRUCTURAL MATRIX
      NBND=NB*ND
      DO 10 I=1,NBND
10    A(I)=0.0D0
C
C SUPERPOSE THE ELEMENT STIFFNESS MATRIX EK(9,9) OF THE
C ELEMENT "NEITR" TO THE STRUCTURAL STIFFNESS MATRIX A(NB*ND)
      DO 20 NEITR=1,NE
C ?????????????????????????????
C      WRITE(1,111) NEITR
C  111 FORMAT(5X,'STSTIF-MATRIX, NEITR=',I4)
C ?????????????????????????????
C
C   COLUMN DETERMINATION
                               CALL ELSTIF
      DO 20 INC=1,3
      INOC=NNP(NEITR,INC)
      IBC=(NNP(NEITR,INC)-1)*3
      DO 20 IDC=1,3
      ICEL=(INC-1)*3+IDC
      ICST=IBC+IDC
      IDI=0
      IF(ICST.GT.NB)IDI=ICST-NB
      IVC=(ICST-1)*NB
C
C   ROW DETERMINATION
      DO 18 INR=1,3
      INOR=NNP(NEITR,INR)
      IF(INOC.LT.INOR)GO TO 18
      IBR=(NNP(NEITR,INR)-1)*3
      IDVC=IDC
      IF(INOC.GT.INOR)IDVC=3
      DO 15 IDR=1,IDVC
      IREL=(INR-1)*3+IDR
      IVV=IVC+IBR+IDR-IDI
      SS=A(IVV)+EK(IREL,ICEL)
C     IF(DABS(SS).LT.1.0D-14)SS=0.0D0
15    A(IVV)=SS
18    CONTINUE
20    CONTINUE
C
C ELIMINATE THE CONSTRAINT POINTS
C      KSTRN(N)....  N-TH CONSTRAINED NODE
C      KSTRT(N)....  TYPE OF N-TH CONSTRAIN
C                      100  X-DISPLACEMENT  IS FIXED
C                      010  Y-DISPLACENENT  IS FIXED
C                      001  Z-ROTATION      IS FIXED
C                      101  X-DIS. & Z-ROT. IS FIXED
      DO 140 N=1,NC
      IRCX=KSTRN(N)*3-2
      IRCY=KSTRN(N)*3-1
      IRCZ=KSTRN(N)*3
      KCHK=KSTRT(N)
C
C  ELIMINATE THE X-DISPLACEMENT
      IF(KCHK.LT.100) GO TO 60
      ICB=ND-IRCX+1
      IF(ICB.GT.NB)ICB=NB
      DO 50 I=1,ICB
      IDI=0
      IRCXX=IRCX+I-1
      IF(IRCXX.GT.NB)IDI=IRCXX-NB
      IXV=(IRCX-2+I)*NB+IRCX-IDI
      IF(I.EQ.1)GO TO 30
      A(IXV)=0.0D0
      GO TO 50
30    A(IXV)=1.0D0
      IF(IRCX.EQ.1)GO TO 50
      LIMF=IRCX-IDI-1
      DO 40 J=1,LIMF
      IXXV=IXV-J
40    A(IXXV)=0.0D0
50    CONTINUE
      KCHK=KCHK-100
C
C  ELIMINATE THE Y-DISPLACEMENT
60    IF(KCHK.LT.010) GO TO 100
      ICB=ND-IRCY+1
      IF(ICB.GT.NB)ICB=NB
      DO 90 I=1,ICB
      IDI=0
      IRCYY=IRCY+I-1
      IF(IRCYY.GT.NB)IDI=IRCYY-NB
      IYV=(IRCY-2+I)*NB+IRCY-IDI
      IF(I.EQ.1)GO TO 70
      A(IYV)=0.0D0
      GO TO 90
70    A(IYV)=1.0D0
      IF(IRCY.EQ.1)GO TO 90
      LIMF2=IRCY-IDI-1
      DO 80 J=1,LIMF2
      IYYV=IYV-J
80    A(IYYV)=0.0D0
90    CONTINUE
      KCHK=KCHK-10
C
C ELIMINATE THE Z-ROTATION
100   IF(KCHK.LT.001) GO TO 140
      ICB=ND-IRCZ+1
      IF(ICB.GT.NB)ICB=NB
      DO 130 I=1,ICB
      IDI=0
      IRCZZ=IRCZ+I-1
      IF(IRCZZ.GT.NB)IDI=IRCZZ-NB
      IZV=(IRCZ-2+I)*NB+IRCZ-IDI
      IF(I.EQ.1)GO TO 110
      A(IZV)=0.0D0
      GO TO 130
110   A(IZV)=1.0D0
      IF(IRCZ.EQ.1)GO TO 130
      LIMF3=IRCZ-IDI-1
      DO 120 J=1,LIMF3
      IZZV=IZV-J
120   A(IZZV)=0.0D0
130   CONTINUE
140   CONTINUE
C
C   TO CALCULATE COLUMN HEIGHTS
C
      DO 160 I=1,ND
      IDI=0
      IF(I.GT.NB)IDI=I-NB
      IIV=(I-1)*NB
      DO 150 J=1,I
      IF(A(IIV+J).EQ.0.0D0)GO TO 150
      MHT(I)=I-J-IDI
      GO TO 160
150   CONTINUE
160   CONTINUE
C
C     PROGRAM
C        TO CALCULATE ADDRESSES OF DIAGONAL ELEMENTS IN
C        BANDED MATRIX WHOSE COLUMN HEIGHTS ARE KNOWN
C
C        MHT  = ACTIVE COLUMN HEIGHTS
C        MAXA = ADDRESSES OF DIAGONAL ELEMENTS
C
C
C
C     CLEAR ARRAY MAXA
C
      NM=ND+1
      DO 170 I=1,NM
170   MAXA(I)=0
C
      MAXA(1)=1
      MAXA(2)=2
      IF(ND.EQ.1)GO TO 190
      DO 180 I=2,ND
180   MAXA(I+1)=MAXA(I)+MHT(I)+1
190   NWK=MAXA(ND+1)-MAXA(1)
C
C   TO STORE STIFFNESS MATRIX A(NB*ND) IN COMPACTED FORM A(NWK)
C
      IAN=0
      DO 200 I=1,ND
      ICK=MAXA(I+1)-MAXA(I)
      IDI=0
      IF(I.GT.NB)IDI=I-NB
      INBB=(I-1)*NB+I-IDI
      DO 200 II=1,ICK
      IAN=IAN+1
      IAV=INBB-II+1
      A(IAN)=A(IAV)
200   CONTINUE
C
C   CLEAR THE REST OF A(NB*ND) ARRAY
C
      LIMS=NWK+1
      DO 210 I=LIMS,NBND
210   A(I)=0.0D0
      RETURN
                                     END
C
C **********************************************************************
C *                                                                    *
                       SUBROUTINE ELSTIF
C *                                                                    *
C **********************************************************************
C TO CALCULATE ELEMENT STIFFNESS MATRIX
C        EK(9,9) ... ELEMENT STIFFNES MATRIX FOR
C        NEITR   ... ELEMENT NUMBER
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/FSK/A(200000),V(2400),EK(9,9),MAXA(2401),NEITR
      COMMON/AIN/NB,NC,ND,NE,NN,TH
      COMMON/B01/B0(1400,9,9),B1(1400,9,9),INXY(1400),ID(1400)
      COMMON/XYZ/NNP(1400,3),X(800),Y(800)
      COMMON/MAT/D0(4,4),D1(4,4),A11,A12,A22,A33,A34,A44,B22,B24,B44
C
C DETERMINE THE NODAL POINTS OF THE ELEMENT N
      I=NNP(NEITR,1)
      J=NNP(NEITR,2)
      M=NNP(NEITR,3)
C
C DETERMINE PARAMETERS
      AI = X(J)*Y(M) - X(M)*Y(J)
      AJ = X(M)*Y(I) - X(I)*Y(M)
      AM = X(I)*Y(J) - X(J)*Y(I)
      BI = Y(J) - Y(M)
      BJ = Y(M) - Y(I)
      BM = Y(I) - Y(J)
C
      CI = X(M) - X(J)
      CJ = X(I) - X(M)
      CM = X(J) - X(I)
C
      DEL= AI + AJ + AM
      TETRA=TH/6.0D0
      RNI=TETRA
      RNJ=TETRA
      RNM=TETRA
      RNN=TH/(DEL+DEL)
C
      A11BI=A11*BI
      A11BJ=A11*BJ
      A11BM=A11*BM
      A12BI=A12*BI
      A12BJ=A12*BJ
      A12BM=A12*BM
      A12CI=A12*CI
      A12CJ=A12*CJ
      A12CM=A12*CM
      A22CI=A22*CI
      A22CJ=A22*CJ
      A22CM=A22*CM
      A33BI=A33*BI
      A33BJ=A33*BJ
      A33BM=A33*BM
      A34BI=A34*BI
      A34BJ=A34*BJ
      A34BM=A34*BM
      A34CI=A34*CI
      A34CJ=A34*CJ
      A34CM=A34*CM
      A44CI=A44*CI
      A44CJ=A44*CJ
      A44CM=A44*CM
C
      A343=A34-A33
      A343BI=A343*BI
      A343BJ=A343*BJ
      A343BM=A343*BM
C
      A434=A44-A34
      A434CI=A434*CI
      A434CJ=A434*CJ
      A434CM=A434*CM
C
      A4343=A44-A34-A34+A33
C
      B22BI=B22*BI
      B22BJ=B22*BJ
      B22BM=B22*BM
      B44CI=B44*CI
      B44CJ=B44*CJ
      B44CM=B44*CM
C
      BCII=BI*CI+BI*CI
      BCIJ=BI*CJ+BJ*CI
      BCIM=BI*CM+BM*CI
      BCJJ=BJ*CJ+BJ*CJ
      BCJM=BJ*CM+BM*CJ
      BCMM=BM*CM+BM*CM
C
      AUTO=TH*DEL/12.0D0
      CROS=TH*DEL/24.0D0
      RNINI=AUTO
      RNJNJ=AUTO
      RNMNM=AUTO
      RNINJ=CROS
      RNJNM=CROS
      RNMNI=CROS
C
C CONSTRUCT ELEMENT STIFFNESS MATRIX EK(9,9)
      EK(1,1) = (A11BI*BI + A44CI*CI)*RNN
      EK(1,4) = (A11BI*BJ + A44CI*CJ)*RNN
      EK(1,7) = (A11BI*BM + A44CI*CM)*RNN
      EK(4,4) = (A11BJ*BJ + A44CJ*CJ)*RNN
      EK(4,7) = (A11BJ*BM + A44CJ*CM)*RNN
      EK(7,7) = (A11BM*BM + A44CM*CM)*RNN
C
      EK(2,4) = (A12CI*BJ + A34BI*CJ)*RNN
      EK(2,7) = (A12CI*BM + A34BI*CM)*RNN
      EK(5,7) = (A12CJ*BM + A34BJ*CM)*RNN
C
      EK(1,2) = (A12BI*CI + A34CI*BI)*RNN
      EK(1,5) = (A12BI*CJ + A34CI*BJ)*RNN
      EK(1,8) = (A12BI*CM + A34CI*BM)*RNN
      EK(4,5) = (A12BJ*CJ + A34CJ*BJ)*RNN
      EK(4,8) = (A12BJ*CM + A34CJ*BM)*RNN
      EK(7,8) = (A12BM*CM + A34CM*BM)*RNN
C
      EK(2,2) = (A22CI*CI + A33BI*BI)*RNN
      EK(2,5) = (A22CI*CJ + A33BI*BJ)*RNN
      EK(2,8) = (A22CI*CM + A33BI*BM)*RNN
      EK(5,5) = (A22CJ*CJ + A33BJ*BJ)*RNN
      EK(5,8) = (A22CJ*CM + A33BJ*BM)*RNN
      EK(8,8) = (A22CM*CM + A33BM*BM)*RNN
C
      EK(1,3) =                A434CI*RNI
      EK(1,6) =                A434CI*RNJ
      EK(1,9) =                A434CI*RNM
      EK(4,6) =                A434CJ*RNJ
      EK(4,9) =                A434CJ*RNM
      EK(7,9) =                A434CM*RNM
C
      EK(2,3) =                A343BI*RNI
      EK(2,6) =                A343BI*RNJ
      EK(2,9) =                A343BI*RNM
      EK(5,6) =                A343BJ*RNJ
      EK(5,9) =                A343BJ*RNM
      EK(8,9) =                A343BM*RNM
C
      EK(3,4) =                A434*CJ*RNI
      EK(3,7) =                A434*CM*RNI
      EK(6,7) =                A434*CM*RNJ
C
      EK(3,5) =                A343*BJ*RNI
      EK(3,8) =                A343*BM*RNI
      EK(6,8) =                A343*BM*RNJ
C
      EK(3,3) = A4343*RNINI + (B22BI*BI+BCII*B24+B44CI*CI)*RNN
      EK(3,6) = A4343*RNINJ + (B22BJ*BI+BCIJ*B24+B44CI*CJ)*RNN
      EK(3,9) = A4343*RNMNI + (B22BM*BI+BCIM*B24+B44CI*CM)*RNN
      EK(6,6) = A4343*RNJNJ + (B22BJ*BJ+BCJJ*B24+B44CJ*CJ)*RNN
      EK(6,9) = A4343*RNJNM + (B22BJ*BM+BCJM*B24+B44CJ*CM)*RNN
      EK(9,9) = A4343*RNMNM + (B22BM*BM+BCMM*B24+B44CM*CM)*RNN
C
      DO 5 IR=1,9
      DO 5 IC=1,9
      AA=EK(IR,IC)
    5 EK(IC,IR)=AA
C
C CONSTRUCT B0 AND B1-MATRIX FOR LATER USE IN SUBROUTINE STRESS
C  CLEAR B0 AND B1-MATRIX
      DO 10 IR=1,9
      DO 10 IC=1,9
      B0(NEITR,IR,IC)=0.0D0
   10 B1(NEITR,IR,IC)=0.0D0
C
C CHECK THE LOCATION OF INTEREST FOR OUTPUT
      IF(INXY(NEITR).EQ.100) GO TO 20
      IF(INXY(NEITR).EQ.010) GO TO 30
      IF(INXY(NEITR).EQ.001) GO TO 40
      XX=(X(I)+X(J)+X(M))/3.0D0
      YY=(Y(I)+Y(J)+Y(M))/3.0D0
      GO TO 50
   20 XX=(X(I)+X(J))/2.0D0
      YY=(Y(I)+Y(J))/2.0D0
      GO TO 50
   30 XX=(X(J)+X(M))/2.0D0
      YY=(Y(J)+Y(M))/2.0D0
      GO TO 50
   40 XX=(X(M)+X(I))/2.0D0
      YY=(Y(M)+Y(I))/2.0D0
C
   50 VNI=AI+BI*XX+CI*YY
      VNJ=AJ+BJ*XX+CJ*YY
      VNM=AM+BM*XX+CM*YY
C B0-MATRIX
      B0(NEITR,1,1)=  BI/DEL
      B0(NEITR,1,4)=  BJ/DEL
      B0(NEITR,1,7)=  BM/DEL
      B0(NEITR,2,2)=  CI/DEL
      B0(NEITR,2,5)=  CJ/DEL
      B0(NEITR,2,8)=  CM/DEL
      B0(NEITR,3,2)=  BI/DEL
      B0(NEITR,3,3)=-VNI/DEL
      B0(NEITR,3,5)=  BJ/DEL
      B0(NEITR,3,6)=-VNJ/DEL
      B0(NEITR,3,8)=  BM/DEL
      B0(NEITR,3,9)=-VNM/DEL
      B0(NEITR,4,1)=  CI/DEL
      B0(NEITR,4,3)= VNI/DEL
      B0(NEITR,4,4)=  CJ/DEL
      B0(NEITR,4,6)= VNJ/DEL
      B0(NEITR,4,7)=  CM/DEL
      B0(NEITR,4,9)= VNM/DEL
C
C B1-MATRIX
      B1(NEITR,1,3)=BI/DEL
      B1(NEITR,1,6)=BJ/DEL
      B1(NEITR,1,9)=BM/DEL
      B1(NEITR,3,3)=CI/DEL
      B1(NEITR,3,6)=CJ/DEL
      B1(NEITR,3,9)=CM/DEL
      RETURN
      END
C
C **********************************************************************
C *                                                                    *
                        SUBROUTINE LOADER
C *                                                                    *
C **********************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/CSTRN/KSTRN(80),KSTRT(80)
      COMMON/AIN/NB,NC,ND,NE,NN,TH
      COMMON/FSK/A(200000),V(2400),EK(9,9),MAXA(2401),NEITR
C
C SET THE EXTERNAL LOAD TO BE ZERO FOR EACH DEGREE OF FREEDOM
      DO 10 LO=1,ND
   10 V(LO)=0.0D0
C
C LOAD THE STRUCTURE
      SCALE=1.0D0
      V( 10)= 22.50D0
      V(  4)= 60.0D0
      V(  1)= 67.5D0
      V( 31)= 60.0D0
      V( 28)= 60.0D0
      V( 25)= 60.0D0
      V( 22)= 30.0D0
C
C PRINT THE EXTERNAL LOAD
      WRITE(6,1000)
 1000 FORMAT(//26X,'**** EXTERNAL LOAD ****')
      WRITE(6,2000)
 2000 FORMAT(/1X,'NODE',4X,'X-FORCE',4X,'Y-FORCE',4X,'Z-MOMNT',
     *        1X,'NODE',4X,'X-FORCE',4X,'Y-FORCE',4X,'Z-MOMNT')
      IND=-1
      INE=0
      DO 20 N=1,ND,6
      IND=IND+2
      INE=INE+2
      IF(N.EQ.(ND-2))GO TO 15
      CHECK= DABS(V(N))+DABS(V(N+1))+DABS(V(N+2))
     1      +DABS(V(N+3))+DABS(V(N+4))+DABS(V(N+5))
      IF(CHECK .EQ. 0.0D0) GO TO 20
      WRITE(6,3000) IND,V(N),V(N+1),V(N+2),
     1               INE,V(N+3),V(N+4),V(N+5)
      GO TO 20
15    CHECK= DABS(V(N))+DABS(V(N+1))+DABS(V(N+2))
      IF(CHECK .EQ. 0.0D0) GO TO 20
      WRITE(6,3000) IND,V(N),V(N+1),V(N+2)
   20 CONTINUE
 3000 FORMAT(2(2X,I3,1X,E10.3,1X,E10.3,1X,E10.3))
C
      RETURN
                                    END
C
C **********************************************************************
C *                                                                    *
                       SUBROUTINE COLSOL
C *                                                                    *
C **********************************************************************
C************************************************************************
C*                                                                      *
C*    PROGRAM                                                           *
C*       TO SOLVE FINITE ELEMENT STATIC EQUILIBRIUM EQUATIONS IN        *
C*       CORE, USING COMPACTED STORAGE AND COLUMN REDUCTION SCHEME      *
C*                                                                      *
C*    --INPUT VARIABLES--                                               *
C*       A(NWK)     = STIFFNESS MATRIX STORED IN COMPACTED FORM         *
C*       V(ND)      = RIGHT-HAND-SIDE LOAD VECTOR                       *
C*       MAXA(ND+1) = VECTOR CONTAINING ADDRESSES OF DIAGONAL           *
C*                    ELEMENTS OF STIFFNESS MATRIX IN A                 *
C*       ND         = NUMBER OF EQUATIONS                               *
C*       NWK        = NUMBER OF ELEMENTS BELOW SKYLINE OF MATRIX        *
C*                                                                      *
C*    --OUTPUT--                                                        *
C*       A(NWK)     = D AND L - FACTORS OF STIFFNESS MATRIX             *
C*       V(ND)      = DISPLACEMENT VECTOR                               *
C*                                                                      *
C************************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/AIN/NB,NC,ND,NE,NN,TH
      COMMON/FSK/A(200000),V(2400),EK(9,9),MAXA(2401),NEITR
C
C     PERFORM L*D*L(T) FACTORIZATION OF STIFFNESS MATRIX
C
      DO 140 N=1,ND
      KN=MAXA(N)
      KL=KN+1
      KU=MAXA(N+1)-1
      KH=KU-KL
      IF(KH)110,90,50
50    K=N-KH
      IC=0
      KLT=KU
      DO 80 J=1,KH
      IC=IC+1
      KLT=KLT-1
      KI=MAXA(K)
      NND=MAXA(K+1)-KI-1
      IF(NND)80,80,60
60    KK=MIN0(IC,NND)
      C=0.0D0
      DO 70 L=1,KK
70    C=C+A(KI+L)*A(KLT+L)
      A(KLT)=A(KLT)-C
80    K=K+1
90    K=N
      B=0.0D0
      DO 100 KK=KL,KU
      K=K-1
      KI=MAXA(K)
      C=A(KK)/A(KI)
      B=B+C*A(KK)
100   A(KK)=C
      A(KN)=A(KN)-B
110   IF(A(KN))120,120,140
120   WRITE(6,2000)N,A(KN)
      STOP
140   CONTINUE
C
C     REDUCE RIGHT-HAND-SIDE LOAD VECTOR
C
      DO 180 N=1,ND
      KL=MAXA(N)+1
      KU=MAXA(N+1)-1
      IF(KU-KL)180,160,160
160   K=N
      C=0.0D0
      DO 170 KK=KL,KU
      K=K-1
170   C=C+A(KK)*V(K)
      V(N)=V(N)-C
180   CONTINUE
C
C     BACK-SUBSTITUTE
C
      DO 200 N=1,ND
      K=MAXA(N)
200   V(N)=V(N)/A(K)
      IF(ND.EQ.1)RETURN
      N=ND
      DO 230 L=2,ND
      KL=MAXA(N)+1
      KU=MAXA(N+1)-1
      IF(KU-KL)230,210,210
210   K=N
      DO 220 KK=KL,KU
      K=K-1
220   V(K)=V(K)-A(KK)*V(N)
230   N=N-1
      WRITE(6,1000)
 1000 FORMAT(//25X,'*** NODAL DISPLACEMENT ***')
      WRITE(6,1500)
 1500 FORMAT(/1X,'NODE',4X,'X-DISP',5X,'Y-DISP',5X,'Z-ROTN',
     *        3X,'NODE',4X,'X-DISP',5X,'Y-DISP',5X,'Z-ROTN')
      IND=-1
      INE=0
      DO 250 K=1,ND,6
      IND=IND+2
      INE=INE+2
      IF(K.EQ.(ND-2))GO TO 240
      WRITE(6,3000) IND,V(K  ),V(K+1),V(K+2),
     *              INE,V(K+3),V(K+4),V(K+5)
      GO TO 250
240   WRITE(6,3000) IND,V(K),V(K+1),V(K+2)
250   CONTINUE
      RETURN
2000  FORMAT(//48H STOP - STIFFNESS MATRIX NOT POSITIVE DEFINITE  ,//
     *           32H NONPOSITIVE PIVOT FOR EQUATION ,I4,//
     *           10H PIVOT =  ,E20.12)
3000  FORMAT(2(2X,I3,1X,E10.3,1X,E10.3,1X,E10.3))
      END
C
C **********************************************************************
C *                                                                    *
                       SUBROUTINE STRESS
C *                                                                    *
C **********************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/B01/B0(1400,9,9),B1(1400,9,9),INXY(1400),ID(1400)
      COMMON/XYZ/NNP(1400,3),X(800),Y(800)
      COMMON/MAT/D0(4,4),D1(4,4),A11,A12,A22,A33,A34,A44,B22,B24,B44
      COMMON/AIN/NB,NC,ND,NE,NN,TH
      COMMON/FSK/A(200000),V(2400),EK(9,9),MAXA(2401),NEITR
      COMMON/DSP/EDISP(9),E(1400,4),PHIJ(1400,4),T(1400,4),CM(1400,4),
     1           U(2400)
C
C
C CLEAR D0- AND D1-MATRIX
      DO 10 I=1,4
      DO 10 J=1,4
      D0(I,J)=0.0D0
   10 D1(I,J)=0.0D0
C
C D0-MATRIX
      D0(1,1)=A11
      D0(1,2)=A12
      D0(2,1)=A12
      D0(2,2)=A22
      D0(3,3)=A33
      D0(3,4)=A34
      D0(4,3)=A34
      D0(4,4)=A44
C
C D1-MATRIX
      D1(2,2)=B22
      D1(2,4)=B24
      D1(4,2)=B24
      D1(4,4)=B44
C
C PICK UP ELEMENT NODAL DISPLACEMENTS EDISP( 9)
C FROM STRUCTURAL NODAL DISPLACEMENTS SDISP(ND)
C
C FOR EACH ELEMENT "IE"
      DO 300 IE=1,NE
       DO 20 IJM=1,3
       IEB=(IJM-1)*3
       ISB=(NNP(IE,IJM)-1)*3
        DO 20 IDOF=1,3
   20   EDISP(IEB+IDOF)=V(ISB+IDOF)
C
C DISPLACEMENT STRAIN : EIJ=B0*UE
       DO 40 IC=1,4
       SUM=0.0D0
        DO 30 K=1,9
   30   SUM=SUM+B0(IE,IC,K)*EDISP(K)
   40  E(IE,IC)=SUM
C
C ROTATION STRAIN : PHI,J=B1*UE
       DO 60 IC=1,4
       SUM=0.0D0
        DO 50 K=1,9
   50   SUM=SUM+B1(IE,IC,K)*EDISP(K)
   60  PHIJ(IE,IC)=SUM
C
C FORCE STRESS : TIJ=D0*EIJ
       DO 80 IC=1,4
       SUM=0.0D0
        DO 70 K=1,4
   70   SUM=SUM+D0(IC,K)*E(IE,K)
   80  T(IE,IC)=SUM
C
C COUPLE STRESS : MIJ=D1*PHIJ
       DO 100 IC=1,4
       SUM=0.0D0
        DO 90 K=1,4
   90   SUM=SUM+D1(IC,K)*PHIJ(IE,K)
  100  CM(IE,IC)=SUM
C STRAIN-ENERGY : U
        SUM=0.0D00
        DO 200 IC=1,4
  200   SUM=SUM+E(IE,IC)*T(IE,IC)+PHIJ(IE,IC)*CM(IE,IC)
       U(IE)=SUM*0.50D00
  300 CONTINUE
C
C WRITE THE CALCULATED STRAINS AND STRESSES
      WRITE(6,1000)
 1000 FORMAT(//19X,'**** STRESSES & STRAINS CALCULATED ****')
      WRITE(6,2000)
 2000 FORMAT(/,'ELMT',1X,'COMP',1X,'DISP-STRAN',3X,'FORCE-STRS'
     *      ,2X,'COMP',1X,'ROTATN-GRAD',2X,'COUPLE-STRS'
     *      ,2X,'STRN-ENEGY')
      WRITE(6,3000)
 3000 FORMAT(14X,'E',12X,'T',15X,'PHI,J',9X,'M',12X,'U')
      DO 400 IE=1,NE
      IF(IE.GT.34.AND.IE.LT.900) GO TO 400
      WRITE(6,4000) IE
 4000 FORMAT(I6)
      WRITE(6,5000) E(IE,1),T(IE,1),PHIJ(IE,1),CM(IE,1),U(IE)
 5000 FORMAT(6X,'XX',1X,1PE11.4,2X,1PE11.4,2X,'Z,X',1X,1PE11.4,
     *       2X,1PE11.4,2X,1PE11.4)
      WRITE(6,6000) E(IE,2),T(IE,2),PHIJ(IE,2),CM(IE,2)
 6000 FORMAT(6X,'YY',1X,1PE11.4,2X,1PE11.4,2X,'X,Z',1X,1PE11.4,
     *       2X,1PE11.4)
      WRITE(6,7000) E(IE,3),T(IE,3),PHIJ(IE,3),CM(IE,3)
 7000 FORMAT(6X,'XY',1X,1PE11.4,2X,1PE11.4,2X,'Z,Y',1X,1PE11.4,
     *       2X,1PE11.4)
      WRITE(6,8000) E(IE,4),T(IE,4),PHIJ(IE,4),CM(IE,4)
 8000 FORMAT(6X,'YX',1X,1PE11.4,2X,1PE11.4,2X,'Y,Z',1X,1PE11.4,
     *       2X,1PE11.4)
  400 CONTINUE
      RETURN
      END
OK,