C      PROGRAM PLFNL2 BY MIKE FOLSE   LAST MODIFIED 12/19/01
C      ADDED DOUBLE PRECISION DUE TO NUMERICAL STABILITY PROBLEMS
C      STATIC ANALYSIS OF PLANE FRAMES AND TRUSSES,
C           INCLUDING P-DELTA EFFECTS  
C              (NOT VERY LARGE DISPLACEMENTS)
C       USE CONSISTENT UNITS
C
C      REVISED 1989 TO INCLUDE P-DELTA EFFECTS ON FIXED END
C         ACTIONS
C
C      INPUT FROM FILE PLFNL2.DAT OR FROM KEYBOARD
C       YOU MAY BUILD PLFNL2.DAT IN AN EDITOR OR WITHIN PROGRAM
C       PROGRAM AUTOMATICALLY STORES KEYBOARD INPUT IN PLFNL2.DAT
C      OUTPUT IS TO FILE PLFNL2.OUT
C
C      INPUT LIST FOR FILE PLFNL2.DAT(FREE FORMAT-SEPERATE INPUT
C           ITEMS WITH ONE COMMA OR ONE OR MORE BLANKS,INPUT
C           ALL ITEMS EVEN IF ZERO)
C
C      LINE 1: TITLE(15A4)
C      LINE 2: NUMBER OF JOINTS(NJT),NUMBER OF MEMBERS(NMEM),
C              NUMBER OF JOINTS WITH APPLIED LOADS(NLDJT),
C              MODULUS OF ELASTICITY(E),
C              NUMBER OF PIN END MEMBERS(NPE)
C               NPE=-1 INDICATES ALL MEMBERS PINNED IE. A TRUSS
C              NPD = P-DELTA CODE  =0 NO P-DELTA 
C                                  =1 P-DELTA
C
C      LINE 3 (REPEAT NJT TIMES): JOINT NUMBER, X COORDINATE,
C             Y COORDINATE, SUPPORT CODE(NS)
C          (=1 IF AT LEAST 1 DOF AT JOINT IS SUPPORT,0 OTHERWISE)
C         NOTE:INPUT JOINTS IN ASCENDING ORDER,MISSING JOINTS WILL
C               BE GENERATED BASED ON EQUAL SPACING
C      LINE 3A  IF NS=1,INPUT SPRING SUPPORT CONSTANTS IN X,Y, AND
C              Z ROTATION. 
C            IF NEGATIVE NUMBER INPUT,DEFAULT = 1.E30 ASSIGNED
C      LINE 4 (REPEAT NMEM TIMES): MEMBER NUMBER, JOINT NUMBER OF
C             JOINT TO BE CONSIDERED J NODE THIS MEMBER,
C             JOINT NUMBER OF K NODE THIS MEMBER,
C             CROSS-SECTION AREA THIS MEMBER, 
C             MOMENT OF INERTIA ABOUT BENDING AXIS,
C             MEMBER LOAD CODE(MLC)
C               (=1 IF MEMBER HAS LOAD, =0 OTHERWISE)
C      LINE 4A   IF MLC=1 INPUT VALUE OF DISTRIBUTED LOAD AT J NODE
C               THEN VALUE AT K NODE, 
C               THEN NUMBER OF CONCENTRATED LOADS (3 MAX)
C               INPUT POSITIVE LOAD IN DOWN DIRECTION
C      LINE 4B IF THERE ARE CONC LOADS,INPUT
C                ( LOAD VALUE , DISTANCE FROM J) FOR EACH LOAD
C      LINE 5 (REPEAT NLDJT TIMES): JOINT NUMBER OF LOADED JOINT,
C             APPLIED LOAD IN X DIRECTION,
C             APPLIED LOAD IN Y DIRECTION,
C             APPLIED COUNTERCLOCKWISE MOMENT ABOUT Z AXIS
C      LINE 6  (IF NPE .NE. 0) INPUT MEMBER NUMBERS OF
C                    PIN ENDED MEMBERS
       IMPLICIT REAL*8 (A-H, O-Z)
       DIMENSION X(98),Y(98),JEND(130),KEND(130),RI(130),D(6),
     1  SM(6,6),RT(6,6),RTRANS(6,6),SX(98),SY(98),SRZ(98),
     1  B(294),TITLE(15),AX(130),ALEN(130),SMR(6,6),
     1IN(6),SMRT(6,6),CXX(130),CYY(130),EACT(6),NS(98),MLC(130),
     1 WL(130),WR(130),FEM(6),P(6),NP(130),NPP(130),NCLD(130),
     1 PV(130,3),PD(130,3),AF(130),AFP(130),FEMS(130,6)
      COMMON/AA/A(294,63)
      DIMENSION AJ(294)
      OPEN(6,FILE='PLFNL2.OUT',STATUS='UNKNOWN')
C      PROGRAM DIMENSIONED FOR MAXIMUM 98 JOINTS AND 130 MEMBERS
      IDIM=294
      JDIM=63
   30 CONTINUE
      NCD=0
      IB=0
C     MAXIMUM DIFFERENCE IN END NODE NUMBERS=20
      INP=5
      IOUT=6
      IP=7
      NGMP=0
      NGM=0
      WRITE(*,5000)
 5000 FORMAT(5X,'PROGRAM PLFNL2 BY MIKE FOLSE, P-DELTA ANALYSIS', 
     1' PLANE FRAMES',/,5X,'SELECT INPUT MODE',
     1 //,15X,'(1)  KEYBOARD',
     1//,15X,'(0)  FILE PLFNL2.DAT',//,15X,'(2)  TERMINATE',/)
      READ(*,*)NCD
      IF(NCD.NE.0.AND.NCD.NE.1)STOP
      WRITE(*,5005)
 5005 FORMAT(/,5X,'OUTPUT IS TO FILE PLFNL2.OUT',//,5X,
     1 ' "EDIT PLFNL2.OUT" TO VIEW CONTENTS',//,5X,
     1 ' "PRINT PLFNL2.OUT" TO GET PRINTED OUTPUT')
      IF(NCD.EQ.1)GO TO 50
      OPEN(5,FILE='PLFNL2.DAT',STATUS='OLD',ACCESS='SEQUENTIAL')
      READ(INP,1000)TITLE
      GO TO 60
   50 CONTINUE
      OPEN(7,FILE='PLFNL2.DAT',STATUS='UNKNOWN')
      WRITE(*,5010)
 5010 FORMAT(/,5X,'ENTER TITLE UP TO 60 CHARACTERS',/)
      READ(*,1000)TITLE
      WRITE(IP,1001)TITLE
 1001 FORMAT(1X,15A4)
 1000 FORMAT(15A4)
   60 CONTINUE
      WRITE(IOUT,2000)TITLE
 2000 FORMAT(//,10X,15A4)
      IF(NCD.EQ.1)GO TO 80
      READ(INP,*)NJT,NMEM,NLDJT,E,NPE,NPD
      GO TO 85
   80 CONTINUE
      WRITE(*,5020)
 5020 FORMAT(/,5X,'ENTER NUMBER OF JOINTS(98 MAXIMUM)',/)
      READ(*,*)NJT
      WRITE(*,5021)
 5021 FORMAT(/,5X,'ENTER NUMBER OF MEMBERS(130 MAXIMUM)',/)
      READ(*,*)NMEM
      WRITE(*,5022)
 5022 FORMAT(/,5X,'ENTER NUMBER OF JOINTS WITH APPLIED LOADING',/)
      READ(*,*)NLDJT
      WRITE(*,5023)
 5023 FORMAT(/,5X,'ENTER MODULUS OF ELASTICITY,E (USE CONSISTENT',
     1 ' UNITS)',/)
      READ(*,*)E
      WRITE(*,5025)
 5025 FORMAT(/,5X,'ENTER NUMBER OF PIN ENDED MEMBERS',
     1 /,10X,'VALUE OF -1 INDICATES ALL PINNED(A TRUSS)',/)
      READ(*,*)NPE
      WRITE(*,5078)
 5078 FORMAT(/,5X,'DO YOU WANT TO CALCULATE P-DELTA EFFECTS?',
     1' (0=NO,1=YES)',/)
      READ(*,*)NPD
      WRITE(IP,1010)NJT,NMEM,NLDJT,E,NPE,NPD
 1010 FORMAT(3I5,E12.3,2I5)
      WRITE(*,5024)
 5024 FORMAT(///,5X,'INPUT JOINT COORDINATES AND SPRING ',
     1 'SUPPORT CONSTANTS',//)
   85 CONTINUE
      WRITE(IOUT,2010)NJT,NMEM,NLDJT,E,NPE,NPD
 2010 FORMAT(//,10X,'NUMBER OF JOINTS =',I5,
     1 /,10X,'NUMBER OF MEMBERS =',I5,
     1 /,10X,'NUMBER OF JOINTS WITH APPLIED LOAD(S) =',I5,
     1/,10X,'MODULUS OF ELASTICITY =',E12.3,/,10X,
     1'NUMBER OF PIN END MEMBERS =',I5,'  ( -1 INDICATES TRUSS)',/
     110X,'NPD=',I5,'(1 MEANS CALCULATE P-DELTA EFFECTS)'
     1 //,40X,'SPRING CONSTANTS',/,1X,'JOINT       X           Y',
     1 '           SX         SY         SRZ',/)
      NUMEQ=3*NJT
      DO 87 I=1,NJT
      SX(I)=0.
      SY(I)=0.
      SRZ(I)=0.
   87 CONTINUE
      MM=0
   88 CONTINUE
      IF(NCD.EQ.1)GO TO 90
      READ(INP,*)M,X(M),Y(M),NS(M)
      IF(NS(M).EQ.1)READ(INP,*)SX(M),SY(M),SRZ(M)
      GO TO 100
   90 CONTINUE
      WRITE(*,5030)  
 5030 FORMAT(/,2X,'ENTER JOINT NUMBER,X COORD,',
     1'Y COORD AND SUPPORT CODE (=1 IF ANY SUPPORTS)',/)
      READ(*,*)M,X(M),Y(M),NS(M)
      IF(NS(M).EQ.1)WRITE(*,5031)
 5031 FORMAT(/,5X,'ENTER SPRING SUPPORT CONSTANTS X,Y,Z',/,
     19X,'(ENTER -1 FOR RIGID SUPPORT,O IF DOF NOT SUPPORTED)',/)
      IF(NS(M).EQ.1)READ(*,*)SX(M),SY(M),SRZ(M)
      WRITE(IP,1020)M,X(M),Y(M),NS(M)
 1020 FORMAT(I5,2E13.4,I5)
      IF(NS(M).EQ.1)WRITE(IP,1021)SX(M),SY(M),SRZ(M)
 1021 FORMAT(6E12.3)
  100 CONTINUE
      IF(SX(M).LT.0.)SX(M)=1.E30
      IF(SY(M).LT.0.)SY(M)=1.E30
      IF(SRZ(M).LT.0.)SRZ(M)=1.E30
      IF(M.EQ.(MM+1))GO TO 116
C     GENERATE MISSING NODES
      M1=MM+1
      M2=M-1
      DIFF=FLOAT(M-MM)
      XINC=(X(M)-X(MM))/DIFF
      YINC=(Y(M)-Y(MM))/DIFF
      AINC=0.
      DO 115 K=M1,M2
      AINC=AINC+1.
      X(K)=X(MM)+AINC*XINC
      Y(K)=Y(MM)+AINC*YINC
C      NOTE:GENERATED NODES CANNOT BE SUPPORTS
      NS(K)=0
  115 CONTINUE
  116 IF(M.EQ.NJT)GO TO 119
      MM=M
      GO TO 88
  119 CONTINUE
      DO 120 I=1,NJT
  120 WRITE(IOUT,2020)I,X(I),Y(I),SX(I),SY(I),SRZ(I)
 2020 FORMAT(1X,I5,2E13.4,3E12.3)
      DO 140 I=1,NMEM
      NP(I)=0
      IF(NCD.EQ.1)GO TO 130
      READ(INP,*)M,JEND(M),KEND(M),AX(M),RI(M),MLC(M)
      IF(MLC(M).EQ.0)GO TO 135
      READ(INP,*)WL(M),WR(M),NCLD(M)
      NCL=NCLD(M)
      IF(NCLD(M).GT.0)READ(INP,*)(PV(M,J),PD(M,J),J=1,NCL)
      GO TO 135
  130 CONTINUE
      WRITE(*,5040)
 5040 FORMAT(/,5X,'ENTER MEM NUM,J NODE,K NODE,AREA,INERTIA',
     1 ' AND LOAD CODE',
     1 /,6X,'(LOAD CODE =1 IF ANY MEMBER LOAD, 0 OTHERWISE)',/)
      READ(*,*)M,JEND(M),KEND(M),AX(M),RI(M),MLC(M)
      WRITE(IP,1030)M,JEND(M),KEND(M),AX(M),RI(M),MLC(M)
 1030 FORMAT(3I5,2E12.3,I5)
      IF(MLC(M).EQ.0)GO TO 135
      WRITE(*,5042)
 5042 FORMAT(/,2X,'ENTER DISTRIBUTED LOAD VALUE J END AND K END',
     1 /,5X,'(-Y LOCAL AXIS IS POSITIVE)',/,5X,
     1 'AND NUMBER OF CONCENTRATED LOADS(3 MAX)',/)
      READ(*,*)WL(M),WR(M),NCLD(M)
      WRITE(IP,6021)WL(M),WR(M),NCLD(M)
 6021 FORMAT(2E12.3,I5)
      NCL=NCLD(M)
      IF(NCL.EQ.0)GO TO 135
      WRITE(*,5043)
 5043 FORMAT(/,2X,'ENTER (LOAD VALUE,DISTANCE FROM J NODE)FOR',
     1 ' UP TO 3 CONC LOADS')
      READ(*,*)(PV(M,J),PD(M,J),J=1,NCL)
      WRITE(IP,1021)(PV(M,J),PD(M,J),J=1,NCL)
  135 CONTINUE
C
C     CALCULATE BAND WIDTH OF STIFFNESS MATRIX
      NBW=3*(IABS(JEND(M)-KEND(M))+1)
      IF(NBW.GT.IB)IB=NBW
  140 CONTINUE
      IF(IB.LE.JDIM)GO TO 150
      WRITE(IOUT,2025)IB
 2025 FORMAT(/,10X,'BAND WIDTH EXCEEDS DIMENSION ALLOWED',
     1 'IB=',I5)
  150 CONTINUE
      WRITE(IOUT,2030)
 2030 FORMAT(//,5X,'MEMBER   J NODE    K NODE    AREA',
     17X, 'INERTIA')
      DO 160 I=1,NMEM
  160 WRITE(IOUT,2040)I,JEND(I),KEND(I),AX(I),RI(I)
 2040 FORMAT(3(5X,I5),4E12.3)
      WRITE(IOUT,2035)
 2035 FORMAT(//,5X,'MEMBER LOADS',/,2X,'MEMBER   W AT J',
     1'   W AT K   P1       D1       P2       D2       P3',
     1'      D3')
      DO 170 I=1,NMEM
      IF(MLC(I).EQ.0)GO TO 170
      NCL=NCLD(I)
      WRITE(IOUT,2041)I,WL(I),WR(I),(PV(I,J),PD(I,J),J=1,NCL)
  170 CONTINUE
 2041 FORMAT(1X,I5,2X,8E9.2)
      IF(IB.GT.JDIM)STOP
C     INITIALIZE LOAD MATRIX TO ZERO
      DO 180 I=1,NUMEQ
  180 B(I)=0.
      IF(NLDJT.EQ.0)GO TO 192
C     INPUT JOINT LOADS
      WRITE(IOUT,2045)
 2045 FORMAT(//,5X,'JOINT      XFORCE    YFORCE      ZMOM',/)
      DO 190 I=1,NLDJT
      IF(NCD.EQ.1)GO TO 185
      READ(INP,*)JT,XFORCE,YFORCE,ZMOM
      GO TO 188
  185 CONTINUE
      WRITE(*,5050)
 5050 FORMAT(/,5X,'ENTER JOINT,FORCE X DIRECTION,',
     1 'FORCE Y DIRECTION, AND CC MOMENT',/)
      READ(*,*)JT,XFORCE,YFORCE,ZMOM
      WRITE(IP,1040)JT,XFORCE,YFORCE,ZMOM
 1040 FORMAT(I5,3E12.3)
  188 CONTINUE
      WRITE(IOUT,2050)JT,XFORCE,YFORCE,ZMOM
 2050 FORMAT(5X,I5,3E12.3)
C     ADD FORCES TO APPLIED LOAD MATRIX
      NXDF=3*JT-2
      NYDF=NXDF+1
      NZDF=NXDF+2
      B(NXDF)=XFORCE
      B(NYDF)=YFORCE
      B(NZDF)=ZMOM
  190 CONTINUE
      DO 191 I =1,NUMEQ
  191 AJ(I)=DBLE(B(I))
  192 CONTINUE
      IF(NPE.LE.0)GO TO 209
      IF(NCD.EQ.1)GO TO 200
      READ(INP,*)(NPP(I),I=1,NPE)
      GO TO 205
  200 CONTINUE
      WRITE(*,5093)
 5093 FORMAT(/,5X,'ENTER MEMBER NUMBERS OF PIN ENDED MEMBERS',/)
      READ(*,*)(NPP(I),I=1,NPE)
      WRITE(IP,5055)(NPP(I),I=1,NPE)
 5055 FORMAT(10I5)
  205 CONTINUE
      WRITE(IOUT,5056)(NPP(I),I=1,NPE)
 5056 FORMAT(//,5X,'THE FOLLOWING MEMBERS ARE PIN ENDED',
     1 (/,3X,16I5))
      DO 207 I=1,NPE
      NP(NPP(I))=1
  207 CONTINUE
  209 CONTINUE
      DO 208 I=1,NMEM
      IF(NPE.EQ.-1)NP(I)=1
      JJ=JEND(I)
      KK=KEND(I)
      XC=X(KK)-X(JJ)
      YC=Y(KK)-Y(JJ)
      AL=SQRT(XC*XC+YC*YC)
      CX=XC/AL
      CY=YC/AL
      CXX(I)=CX
      CYY(I)=CY
C     CX,CY, ARE DIRECTION COSINES FOR MEMBER I, I.E. COSINES OF
C      ANGLES BETWEEN MEMBER AXIS AND GLOBAL X,Y AXES
      ALEN(I)=AL
C     ADD MEMBER END ACTIONS TO LOAD MATRIX (FOR 0 AXIAL LOAD)
      IF(MLC(I).EQ.0)GO TO 208
      CALL FORMRT(RT,RTRANS,CX,CY)
      CALL INDX(JJ,KK,IN)
      CALL FM(WL(I),WR(I),AL,FEM,NP(I))
      NCL=NCLD(I)
      IF(NCL.EQ.0)GO TO 265
      DO 263 J=1,NCL
  263 CALL FM2(PV(I,J),PD(I,J),AL,FEM,NP(I))
  265 CALL MULT(RTRANS,FEM,1,P)
C     STORE FIXED END ACTIONS 
      DO 267 J = 1,6
  267 FEMS(I,J)=FEM(J)
      DO 270 J=1,6
      AJ(IN(J))=AJ(IN(J))-P(J)
  270 CONTINUE
  208 CONTINUE
  210 CONTINUE
      DO 220 I=1,NUMEQ    
C     REINITIALIZE GLOBAL STIFFNESS MATRIX TO ZERO
      DO 220 J=1,IB
      A(I,J)=0.
  220 CONTINUE
C     FORM MEMBER STIFFNESS MATRICES, ROTATE TO GLOBAL COORDINATES,
C       INDEX TO GLOBAL DEGREE OF FREEDOM NUMBERS, AND ADD TO
C       GLOBAL STIFFNESS MATRIX
      DO 300 I=1,NMEM
      CX=CXX(I)
      CY=CYY(I)
      AL=ALEN(I)
      CALL FORMRT(RT,RTRANS,CX,CY)
      CALL FORMSM(SM,AL,AX(I),E,RI(I),NP(I))
      IF(NGM.GT.0.AND.AX(I).GT.0.)CALL GSM(SM,AL,NP(I),AF(I))
      CALL MULT(SM,RT,6,SMRT)
      CALL MULT(RTRANS,SMRT,6,SMR)
      CALL INDX(JEND(I),KEND(I),IN)
C     FORM UPPER HALF OF STIFFNESS MATRIX,STORE IN RECTANGULAR
C     ARRAY WITH DIAGONAL ELEMENT IN COLUMN 1
      DO 260 J=1,6
      JG=IN(J)
      DO 260 K=1,6
      KG=IN(K)
      IF(KG.LT.JG)GO TO 260
      LG=KG-JG+1
      A(JG,LG)=A(JG,LG)+SMR(J,K)
  260 CONTINUE
  300 CONTINUE
C     CORRECT GLOBAL STIFFNESS MATRIX FOR BOUNDARY CONDITIONS
C       I.E. SUPPORT SPRINGS WHERE INPUT
      DO 360 I=1,NJT
      NXDF=3*I-2
      NYDF=NXDF+1
      NZDF=NXDF+2
      IF(NS(I).EQ.0)GO TO 350
      A(NXDF,1)=A(NXDF,1)+SX(I)
      A(NYDF,1)=A(NYDF,1)+SY(I)
      A(NZDF,1)=A(NZDF,1)+SRZ(I)
  350 IF(NPE.EQ.-1)A(NZDF,1)=1.E30
  360 CONTINUE
C     SOLVE   AX=B FOR DISPLACEMENT MATRIX X
      CALL UDU(AJ,NUMEQ,IB,IDIM,JDIM,NER)
      IF(NER.EQ.0)GO TO 370
      WRITE(IOUT,4200)
 4200 FORMAT(/,10X,'STIFFNESS MATRIX IS SINGULAR')
      GO TO 550
  370 CONTINUE
C     NOTE: DEFLECTIONS ARE RETURNED IN B MATRIX WHICH WENT INTO
C         UDU AS APPLIED LOADS
      IF(NGM.GT.0)GO TO 420
      WRITE(6,4203)
 4203 FORMAT(///,5X,'OUTPUT NOT INCLUDING P-DELTA EFFECTS',//)
  375 WRITE(IOUT,4100)
 4100 FORMAT(//,9X,'JOINT DISPLACEMENTS AND REACTIONS',//,5X,'JT '
     1,' X DISPL   Y DISPL    Z ROT      X REACT    Y REACT',
     1'    ZM REACT')
      DO 400 I=1,NJT
      NXDF=3*I-2
      NYDF=NXDF+1
      NZDF=NXDF+2
      XREACT=-SX(I)*AJ(NXDF)
      YREACT=-SY(I)*AJ(NYDF)
      ZREACT=-SRZ(I)*AJ(NZDF)
      WRITE(IOUT,4110)I,AJ(NXDF),AJ(NYDF),AJ(NZDF),
     1 XREACT,YREACT,ZREACT
 4110 FORMAT(1X,I5,3D11.3,3E11.3)
  400 CONTINUE
      WRITE(IOUT,4115)
 4115 FORMAT(//,13X,'MEMBER END ACTIONS',//,1X,'MEMBER J   ',
     1'K     AX J     SHR J      MOM J        AX K     SHR K', 
     1 6X,'MOM K')
  420 CONTINUE
      DO 500 I=1,NMEM
      AL=ALEN(I)
      JJ=JEND(I)
      KK=KEND(I)
      CALL INDX(JJ,KK,IN)
      DO 450 J=1,6
  450 D(J)=AJ(IN(J))
      CALL FORMRT(RT,RTRANS,CXX(I),CYY(I))
      CALL FORMSM(SM,AL,AX(I),E,RI(I),NP(I))
      IF(AX(I).LE.0.)GO TO 455
      IF(NGM.GT.0.OR.NGMP.EQ.1)CALL GSM(SM,AL,NP(I),AF(I))
  455 CONTINUE
      CALL MULT(SM,RT,6,SMRT)
      CALL MULT(SMRT,D,1,EACT)
      AFP(I)=EACT(4)
      IF(NGM.GT.0)GO TO 500
C     ADD FIXED END ACTIONS
      IF(MLC(I).EQ.0)GO TO 480
      DO 460 K=1,6
  460 EACT(K)=EACT(K)+FEMS(I,K)
  480 CONTINUE
      WRITE(IOUT,4120)I,JJ,KK,(EACT(J),J=1,6)
 4120 FORMAT(1X,I3,2I4,6E11.3)
  500 CONTINUE
  550 CONTINUE
      IF(NPD.EQ.0)STOP
      NGM=NGM+1
      IF(NGM.LT.5)GO TO 555
      WRITE(IOUT,*)'  DOES NOT CONVERGE'
      STOP
  555 CONTINUE
      IF(NGM.EQ.1)GO TO 570
C     COMPARE NEW AND OLD VALUES OF MEMBER AXIAL FORCE
      DIFMAX=0.
      DO 560 I=1,NMEM
      PCR=9.87*E*RI(I)/ALEN(I)**2
      RAT=AFP(I)/PCR
      IF(RAT.GT.-.02)GO TO 560
      DIF=ABS((AFP(I)-AF(I))/AFP(I))
      IF(DIF.GT.DIFMAX)DIFMAX=DIF
  560 CONTINUE
      IF(DIFMAX.LT.0.01)GO TO 590
  570 CONTINUE
      DO 580 I=1,NMEM
  580 AF(I)=AFP(I)
C
C     RECREATE LOAD MATRIX INCLUDING EFFECTS OF AXIAL LOAD ON
C       FIXED END ACTIONS
      DO 582 I=1,NUMEQ
C     INITIALIZE WITH JOINT LOADS
      AJ(I)=B(I)
  582 CONTINUE
C     ADD MEMBER FIXED END ACTIONS
C
      DO 585 I = 1,NMEM
      IF(MLC(I).EQ.0)GO TO 585
C     IF AXIAL STRESS IS LESS THAN .1 PCR COMPRESSION, IGNORE
C        EFFECTS OF AXIAL FORCE ON FIXED END ACTIONS
      AL=ALEN(I)
      EI=E*RI(I)
      PCR=9.87*EI/AL**2
      RAT=AFP(I)/PCR
      CX=CXX(I)
      CY=CYY(I)
      JJ=JEND(I)
      KK=KEND(I)
      CALL FORMRT(RT,RTRANS,CX,CY)
      CALL INDX(JJ,KK,IN)
      NCL=NCLD(I)
      IF(RAT.LT.-0.1.AND.AX(I).GT.0.)GO TO 588
      CALL FM(WL(I),WR(I),AL,FEM,NP(I))
      IF(NCL.EQ.0)GO TO 665
      DO 663 J=1,NCL
  663 CALL FM2(PV(I,J),PD(I,J),AL,FEM,NP(I))
  665 CONTINUE
      GO TO 680
  588 CONTINUE
C     INCLUDE AXIAL EFFECTS OF FIXED END ACTIONS
      CALL FMN(WL(I),WR(I),AL,FEM,NP(I),AFP(I),EI)
      IF(NCL.EQ.0)GO TO 675
      DO 670 J=1,NCL
  670 CALL FM2N(PV(I,J),PD(I,J),AL,FEM,NP(I),AFP(I),EI)      
  675 CONTINUE
C      
  680 CALL MULT(RTRANS,FEM,1,P)
C     SAVE FIXED END ACTIONS
      DO 690 J=1,6
  690 FEMS(I,J)=FEM(J)
C
      DO 700 J=1,6
  700 AJ(IN(J))=AJ(IN(J))-P(J)
  585 CONTINUE
  586 CONTINUE
      GO TO 210
  590 CONTINUE
      NPD=0
      NGMP=1
      WRITE(6,5001)NGM
 5001 FORMAT(///,10X,'OUTPUT INCLUDING P-DELTA EFFECTS',
     1 /,5X,'NUMBER OF CYCLES=',I5,//)
      NGM=0       
      GO TO 375
      END
C
C
C     SUBROUTINE UDU
      SUBROUTINE UDU(X,NUMEQ,IB,IDIM,JDIM,NER)
C     UDU  FOR BANDED SYMMETRIC OR  MATRICES
C      BY ERIK THOMPSON OF CSU
      IMPLICIT REAL*8 (A-H, O-Z)
      COMMON/AA/A(294,63)
      DIMENSION X(IDIM)
      NER=0
      NEM1=NUMEQ-1
      DO 450 I=1,NEM1
      JEND=NUMEQ-I+1
      IF(JEND.GT.IB)JEND=IB
      DO 440 J=2,JEND
      J1=J+I-1
      IF(A(I,1).EQ.0.)GO TO 480
      FAC=A(I,J)/A(I,1)
      K1=0
      DO 430 K=J,JEND
      K1=K1+1
      A(J1,K1)=A(J1,K1)-FAC*A(I,K)
  430 CONTINUE
      A(I,J)=FAC
  435 CONTINUE
      X(J1)=X(J1)-A(I,J)*X(I)
  440 CONTINUE
  450 CONTINUE
      X(NUMEQ)=X(NUMEQ)/A(NUMEQ,1)
      DO 470 I=1,NEM1
      I1=NUMEQ-I
      X(I1)=X(I1)/A(I1,1)
      JEND=NUMEQ-I1+1
      IF(JEND.GT.IB)JEND=IB
      DO 460 J=2,JEND
      J1=I1+J-1
      X(I1)=X(I1)-A(I1,J)*X(J1)
  460 CONTINUE
  470 CONTINUE
      RETURN
  480 NER=1
      RETURN
      END
      SUBROUTINE FORMSM(SM,AL,AX,E,RI,NP)
C     STIFFNESS MATRIX PLANE FRAME OR TRUSS
      IMPLICIT REAL*8 (A-H, O-Z)
      DIMENSION SM(6,6)
      DO 10 I = 1,6
      DO 10 N = I,6
      SM(I,N)=0.
   10 SM(N,I)=0.0
      SM(1,1)=ABS(AX)*E/AL
      SM(4,4)=SM(1,1)
      SM(1,4)=-SM(1,1)
      SM(4,1)=SM(1,4)
      SM(3,3)=0.000001*E*RI/AL
      SM(6,6)=SM(3,3)
      IF(NP.EQ.1)RETURN
      SM(2,2)=12.*E*RI/AL**3
      SM(5,5)=SM(2,2)
      SM(2,5)=-SM(2,2)
      SM(2,3)=6.*E*RI/AL**2
      SM(2,6)=SM(2,3)
      SM(3,5)=-SM(2,3)
      SM(5,6)=SM(3,5)
      SM(3,3)=4.*E*RI/AL
      SM(6,6)=SM(3,3)
      SM(3,6)=.5*SM(3,3)
      DO 20 I=1,5
      J=I+1
      DO 20 M=J,6
   20 SM(M,I)=SM(I,M)
      RETURN
      END
      SUBROUTINE FORMRT(RT,RTRANS,CX,CY)
C     FORM ROTATION MATRICES
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION RT(6,6),RTRANS(6,6)
      DO 10 I = 1,6
      DO 10 N = I,6
      RT(I,N)=0.
   10 RT(N,I)=0.0
      RT(1,1)=CX
      RT(2,2)=CX
      RT(4,4)=CX
      RT(5,5)=CX
      RT(1,2)=CY
      RT(4,5)=CY
      RT(2,1)=-CY
      RT(5,4)=-CY
      RT(3,3)=1.
      RT(6,6)=1.
C      FORM RTRANS
      DO 30 I = 1,6
      DO 30 N = 1,6
   30 RTRANS(N,I)=RT(I,N)
      RETURN
      END
      SUBROUTINE MULT(A,B,JB,C)
      IMPLICIT REAL*8 (A-H,O-Z)
C      MATRIX MULTIPLY
      DIMENSION A(6,6),B(6,JB),C(6,JB)
      DO 10 I =1,6
      DO 10 N =1,JB
      SUM=0.0
      DO 20 L=1,6
   20 SUM=SUM+A(I,L)*B(L,N)
   10 C(I,N)=SUM
      RETURN
      END
      SUBROUTINE INDX(JJ,KK,IN)
      DIMENSION IN(6)
C     INDX GIVES GLOBAL DEGREE OF FREEDOM NUMBERS CORRESPONDING
C      TO THE SIX MEMBER LOCAL DEGREES OF FREEDOM 
      J3=3*JJ
      K3=3*KK
      IN(3)=J3
      IN(2)=J3-1
      IN(1)=J3-2
      IN(6)=K3
      IN(5)=K3-1
      IN(4)=K3-2
      RETURN
      END
      SUBROUTINE FM(WL,WR,AL,FEM,NP)
C     SUBROUTINE TO CALC FIXED END ACTIONS FOR LINEAR LOAD
      IMPLICIT REAL*8(A-H, O-Z)
      DIMENSION FEM(6)
      FEM(1)=0.
      FEM(4)=0.
      IF(NP.EQ.1)GO TO 100
      FEM(3)=1./12.*WL*AL**2+1./30.*(WR-WL)*AL**2
      FEM(6)=-1./12.*WL*AL**2-1./20.*(WR-WL)*AL**2
      GO TO 110
  100 CONTINUE
      FEM(3)=0.
      FEM(6)=0.
  110 CONTINUE
      FEM(2)=WL*AL*.5+(WR-WL)*AL*.5/3.+(FEM(3)+FEM(6))/AL
      FEM(5)=WL*AL*.5+(WR-WL)*AL*.333333-(FEM(3)+FEM(6))/AL
      RETURN
      END
      SUBROUTINE FM2(PV,PD,AL,FEM,NP)
C     SUBROUTINE TO CALC FEM DUE TO CONC LOAD
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION FEM(6)
      PDP=AL-PD
      IF(NP.EQ.1)GO TO 100
      FEM3=PV*PD*PDP**2/AL**2
      FEM6=-PV*PD**2*PDP/AL**2
      GO TO 110
  100 CONTINUE
      FEM3=0.
      FEM6=0.
  110 CONTINUE
      FEM(2)=FEM(2)+PV*PDP/AL+(FEM3+FEM6)/AL
      FEM(5)=FEM(5)+PV*PD/AL-(FEM3+FEM6)/AL
      FEM(3)=FEM(3)+FEM3
      FEM(6)=FEM(6)+FEM6
      RETURN
      END
      SUBROUTINE FMN(WL,WR,AL,FEM,NP,PS,EI)
      IMPLICIT REAL*8 (A-H,O-Z)
C     SUBROUTINE TO CALC FIXED END ACTIONS FOR LINEAR LOAD
C      INCLUDING P-DELTA EFFECTS
      DIMENSION FEM(6)
      FEM(1)=0.
      FEM(4)=0.
      IF(NP.EQ.1)GO TO 100
      P=-PS
      RK=SQRT(P/EI)
      RKL=RK*AL
      RL2=AL*AL
      CO=COS(RKL)
      SI=SIN(RKL)
      DENOM=(SI-RKL)*(RK*SI) + RK*(1.-CO)**2
      Q=WR-WL
      A=1./DENOM*(-Q*RL2/6./P*RK*SI +Q*AL/2./P*(1.-CO))
      B=1./DENOM*((SI-RKL)*Q*AL/2./P + RK*(1.-CO)*Q*RL2/6./P)
      U=RK*AL*.5
      F1=3.*(TAN(U)-U)/(U**3)
C     END MOMENTS DUE TO UNIFORM WL
      FEM(3)=1./12.*WL*(AL**2)*F1*U/TAN(U)
      FEM(6)=-FEM(3)
C     END MOMENTS DUE TO TRIANGULAR LOAD 0 AT LT, WR-WL AT RT
C      Q0=WL-WR
C      IF(ABS(Q0).LT.0.001)GO TO 80
C      CON=Q0/(P*AL*SI)*(1./RK**2*SI-AL/RK*CO)
C      TH0A=CON- Q0*AL/3./P 
C      TH0B=-Q0*CO/P/RK/SI + Q0/P/RK/SI - CON - Q0*AL/6./P
C      F2=1.5/U*(.5/U - 1./TAN(2.*U))
C      F3 = 3./U*(1./SIN(2.*U) - .5/U)
C      AA=AL*F2/3./EI
C      BB=AL*F3/6./EI
C      CC=BB
C      DD=AA
C      DENOM=AA*DD-BB*CC
C      RML=(-TH0A*DD+TH0B*BB)/DENOM
C      RMR=(-TH0B*AA+TH0A*CC)/DENOM
      RML=B*P
      RMR=-A*P*SI-B*P*CO+Q/P*EI
      FEM(3)=FEM(3)-RML
      FEM(6)=FEM(6)-RMR
   80 CONTINUE
      GO TO 110
  100 CONTINUE
      FEM(3)=0.
      FEM(6)=0.
  110 CONTINUE
      FEM(2)=WL*AL*.5+(WR-WL)*AL*.5/3.+(FEM(3)+FEM(6))/AL
      FEM(5)=WL*AL*.5+(WR-WL)*AL*.333333-(FEM(3)+FEM(6))/AL
      RETURN
      END
      SUBROUTINE FM2N(PV,PD,AL,FEM,NP,PS,EI)
      IMPLICIT REAL*8 (A-H,O-Z)
C     SUBROUTINE TO CALC FEM DUE TO CONC LOAD
      DIMENSION FEM(6),A(4,4),B(4)
      IF(NP.EQ.1)GO TO 100
      Q=PV
      P=-PS
      AA=PD
      RL=AL
      B(1)=0.
      B(2)=Q/P
      B(3)=Q/P*(AA-RL)
      B(4)=B(2)
      RK=SQRT(P/EI)
      AK=AA*RK
      RKL=RL*RK
      SI=SIN(RKL)
      CO=COS(RKL)
      SIA=SIN(AK)
      COA=COS(AK)
      A(1,1)=-COA/P
      A(1,2)=SIA/P/RK
      A(1,3)=-COA
      A(1,4)=-SIA
      A(2,1)=RK*SIA/P
      A(2,2)=COA/P
      A(2,3)=RK*SIA
      A(2,4)= - RK*COA
      A(3,1)= 1./P
      A(3,2)= -RL/P
      A(3,3)=CO
      A(3,4)=SI
      A(4,1)=0.
      A(4,2)=1./P
      A(4,3)=RK*SI
      A(4,4)=-RK*CO
      CALL LDU2(A,B,1,4,4)
      RML=B(1)
      VL=B(2)
C      SUM MOMENTS ABOUT RIGHT END
      RMR=VL*RL - RML - Q*(RL - AA)
C
C      PREVIOUSLY USED PROCEDURE COMMENTED OUT - BREAKS DOWN IF 
C      AXIAL LOAD GREATER THAN SIMPLE BEAM BUCKLING LOAD
      PDP=AL-PD
C      IF(NP.EQ.1)GO TO 100
C      P=-PS
C      RK=SQRT(P/EI)
C      RL=RK*AL
C      CO=COS(RL)
C      SI=SIN(RL)
C      U=RK*AL*.5
C     END MOMENTS DUE TO CONC LD
C      TH0A=PV*SIN(RK*PDP)/P/SI - PV*PDP/P/AL
C      TH0B= PV*SIN(RK*PD)/P/SI - PV*PD/P/AL
C      F2=1.5/U*(.5/U - 1./TAN(2.*U))
C      F3 = 3./U*(1./SIN(2.*U) - .5/U)
C      AA=AL*F2/3./EI
C      BB=AL*F3/6./EI
C      CC=BB
C      DD=AA
C      DENOM=AA*DD-BB*CC
C      RML=(-TH0A*DD+TH0B*BB)/DENOM
C      RMR=(-TH0B*AA+TH0A*CC)/DENOM
C      FEM3=-RML
      FEM3=RML
      FEM6=RMR
      GO TO 110
  100 CONTINUE
      FEM3=0.
      FEM6=0.
  110 CONTINUE
      FEM(2)=FEM(2)+PV*PDP/AL+(FEM3+FEM6)/AL
      FEM(5)=FEM(5)+PV*PD/AL-(FEM3+FEM6)/AL
      FEM(3)=FEM(3)+FEM3
      FEM(6)=FEM(6)+FEM6
      RETURN
      END
C     SUBROUTINE GSM IS TO CALCULATE P-DELTA EFFECTS
      SUBROUTINE GSM(SM,AL,NP,AF)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION SM(6,6)
      IF(NP.EQ.0)GO TO 50
C      NP = 1   PIN ENDED MEMBER
      SM(2,2) = AF/AL
      SM(2,5) = - SM(2,2)
      SM(5,5) = SM(2,2)
      SM(5,2) = SM(2,5)
      RETURN
   50 CONTINUE
      A=1.2*AF/AL
      B=0.1*AF
      C=2.*AL/15.*AF
      D=AL/30.*AF
      SM(2,2)=SM(2,2)+A
      IF(SM(2,2).LE.0.)WRITE(6,100)     
      IF(SM(2,2).LT.0.)SM(2,2)=0.
      SM(5,5)=SM(2,2)
      SM(3,3)=SM(3,3)+C
      IF(SM(3,3).LE.0.)WRITE(6,100)
  100 FORMAT(5X,' MEMBER STIFFNESS TERM.LT.0.')
      IF(SM(3,3).LT.0.)SM(3,3)=0.
      SM(6,6)=SM(3,3)
      SM(2,5)=-SM(2,2)
      SM(2,3)=SM(2,3)+B
      SM(2,6)=SM(2,3)
      SM(3,5)=-SM(2,3)
      SM(5,6)=SM(3,5)
      SM(3,6)=SM(3,6)-D
      DO 20 I=1,5
      J=I+1
      DO 20 M=J,6
   20 SM(M,I)=SM(I,M)
      RETURN
      END
C     SUBROUTINE LDU2
      SUBROUTINE LDU2(A,X,LU,NUMEQ,IDIM)
      IMPLICIT REAL*8 (A-H,O-Z)
C     LDU2 BEST FOR UNBANDED SYMMETRIC OR NON-SYMMETRIC MATRICES
      DIMENSION A(IDIM,IDIM),X(IDIM)
      IOUT=6
      NEM1=NUMEQ-1
      DO 450 I=1,NEM1
C      WRITE(IOUT,5000)I,A(I,I)
C 5000 FORMAT(2X,I5,F15.3)
      JEND=NUMEQ-I
      DO 440 J=1,JEND
      JP1=J+I
      IF(LU.EQ.0)GO TO 435
      IF(A(I,I).EQ.0.)GO TO 480
      FAC=A(JP1,I)/A(I,I)
      A(JP1,I)=FAC
      DO 430 K=1,JEND
      KP1=K+I
      A(JP1,KP1)=A(JP1,KP1)-FAC*A(I,KP1)
  430 CONTINUE
  435 CONTINUE
      X(JP1)=X(JP1)-A(JP1,I)*X(I)
  440 CONTINUE
  450 CONTINUE
      X(NUMEQ)=X(NUMEQ)/A(NUMEQ,NUMEQ)
      DO 470 I=1,NEM1
      IROW=NUMEQ-I
      DO 460 J=1,I
      JROW=IROW+J
      X(IROW)=X(IROW)-A(IROW,JROW)*X(JROW)
  460 CONTINUE
      X(IROW)=X(IROW)/A(IROW,IROW)  
  470 CONTINUE
      GO TO 490
  480 WRITE(IOUT,4000)
 4000 FORMAT(10X,35HDENOMINATOR EQUAL ZERO-END OF JOB     )
      STOP
  490 CONTINUE
      RETURN
      END

