
C     ******************************************************************
C
      PROGRAM ULARC
C      
C     ******************************************************************
C
C     SAMPLE ELASTO-PLASTIC ANALYSIS OF PLANE FRAMES
C
C     PROGRAMMED BY A. SUDHAKAR, G. H. POWELL, G. ORR, R. WHEATON
C		    UNIVERSITY OF CALIFORNIA, BERKELEY, CALIFORNIA
C		    1972
C
C     ******************************************************************
C
      Implicit None
      character*60 ZFILE1,ZFILE2

      INTEGER IBUST           !   Failure Code
      INTEGER ICLOSE          !   Limi
      INTEGER NCYC
      INTEGER INEL
      INTEGER JNEL
      INTEGER NLC
      INTEGER NELS
      REAL XFAC
      INTEGER NEQ
      INTEGER NJTS
      INTEGER NSJTS,NTERM
      REAL DMAX
      REAL FAX
      INTEGER NLCS,MBAND
      Integer, parameter :: ii = 50

      REAL X(140),Y(140),RRTOT(140,3)

      REAL YMP(ii),YMN(ii),ROTCP(ii),ROTCN(ii)

      REAL PP(420),RR(420)
      INTEGER NA(420)
      REAL R(420)
      REAL AA(ii),E(ii),XI(ii),SKII(ii),SKJJ(ii),SKIJ(ii)
      INTEGER NSFSCI(ii),NSFSCJ(ii)
      REAL SPA(50),SPB(50),SPR(50),COSS(50),SINS(50)
      INTEGER NSPJ(50)
      REAL SPFRA(50),SPFRB(50),SPFRR(50)
      REAL EK(3,3),EKA(3,6),FK(3,3),FKA(3,6),AKA(6,6)
      INTEGER LM(6)
      REAL A(3,6),AB(3,6)
      INTEGER IHI(200),IHJ(200),NODI(200),NODJ(200),NSFTP(200)
      INTEGER NSECI(200),NSECJ(200)
      REAL STOT(200,3),VPTOT(200,3),VPTP(200,3),VPTN(200,3)
      REAL S(200,3),VP(200,3)

      REAL SK(6000)

      INTEGER MAXL,I
      Logical  ::   ex

      COMMON /CONTR/ IBUST,ICLOSE,NCYC,INEL,JNEL,NLC,XFAC,NELS,NEQ,NJTS,
     1NSJTS,DMAX,NTERM,FAX,NLCS,MBAND
      COMMON /JOINT/ X,Y,RRTOT
      COMMON /SECTAB/ YMP,YMN,ROTCP,ROTCN
      COMMON /EQNS/ PP,RR,NA,R
      COMMON /STFTAB/ AA,E,XI,SKII,SKJJ,SKIJ,NSFSCI,NSFSCJ
      COMMON /SUPJT/ SPA,SPB,SPR,COSS,SINS,NSPJ,SPFRA,SPFRB,SPFRR
      COMMON /STIFF/ EK,EKA,FK,FKA,AKA,LM,A,AB
      COMMON /MEMBER/ IHI,IHJ,NODI,NODJ,NSFTP,NSECI,NSECJ,
     2STOT,VPTOT,VPTP,VPTN,S,VP
C
      COMMON SK

      MAXL=6000                   ! Maximum Elements

      write (*,331)
331   format(1x,' This is a Intel Microsoft Fortran Programme ',/
     11x,' amended by J. Nichols from ULARC - UC Berkeley',/
     31x,' Ver. 1.2  4 April 2015'//)
C
C     read input / output files and open files
C
      write(*,5)
    5 format(1x,'File name for input data :'\)
      read (*,6)zFILE1
      inquire(file = zFile1, exist=ex)
      IF (.NOT. ex) THEN
          WRITE (*,'(2A/)') ' >> Cannot find file ', zfile1
          GOTO 99
      END IF


      open (2,file=zfile1,err=99,status='unknown')
    6 format(a12)
      write (*,7)
    7 format(1x,'File name for output data :'\)
      read(*,8)zfile2
    8 format(a12)
      open (1,file=zfile2,status='unknown')
      open(3,file='R.csv',status='unknown')
c
c     set up structure
c
      write(3,331)
10    NLC=0

C---------------------------------------------------------------------------
C
C     Call Input Routine
C
C---------------------------------------------------------------------------

      write (*,332)
332   format(1x,'Input Subroutine')
      CALL INPT

C---------------------------------------------------------------------------
C
C     Call BANDC Routine
C
C---------------------------------------------------------------------------


      NEQ=NJTS*3
      write (*,333)
333   format(1x,'BandC Subroutine ')
      CALL BANDC

C---------------------------------------------------------------------------
C
C     Structure Stiffness Model Data
C
C---------------------------------------------------------------------------


      write(1,20) NA(NEQ),MAXL
   20 FORMAT(/////20H STRUCTURE STIFFNESS//
     110X,19H STORAGE REQUIRED =I6/
     210X,19H STORAGE PROVIDED =I6)

      IF (NA(NEQ).LE.MAXL) GO TO 30
      STOP

   30 NLC=NLC+1
C
C     SET UP LOAD CASE
C
      write (*,334)
334   format(1x,'Load Subroutine')

      IF (NLC.GT.NLCS) GO TO 100
      write (*,329)nlc,nlcs
329   format(1x,'Load Case :',i4,'      Number of Load Cases :',i4)
      NCYC=0
      XFAC=0.
      CALL LOAD
C
C     IGNORE LOAD CASE IF NOT APPROPRIATE
      IF(IBUST.GT.0) GO TO 30
C
C     CYCLE WITHIN LOAD CASE
C
   40 NCYC=NCYC+1
      write (*,335)
335   format(1x,'Stiffness Subroutine')
      IF(NCYC.LE.NTERM) GO TO 60
      write (1,50)
   50 FORMAT(/////31H NUMBER OF CYCLES EXCEEDS LIMIT)
      IBUST=2
      GO TO 30
C
C     FORM STRUCTURE STIFFNESS
C
   60 CALL STIF (SK)
C
C     SOLVE EQUILBRIUM EQUATIONS
C
      DO 70 I=1,NEQ
   70 RR(I)=PP(I)
      write (*,336)
336   format(1x,'Optimum Solution Subroutine')
      CALL OPTSOL (SK,RR,NA,NEQ,1,1,1)
      CALL OPTSOL (SK,RR,NA,NEQ,1,1,2)
C
C     DETERMINE LOAD FACTOR FOR NEXT EVENT,UPDATE STATUS
C
      IBUST=0
      ICLOSE=0
      write (*,337)
337   format(1x,'Force Subroutine')
      CALL FORCE
C
C     OUTPUT RESULTS
C
      IF(ICLOSE.EQ.1) GO TO 40
      IF(IBUST.EQ.1) GO TO 30
      write (*,338)
338   format(1x,'Output Subroutine')
      CALL OUTPT
      IF(IBUST.NE.0) GO TO 30
      GO TO 40
100   write (*,339)
339   format(1x,'We are all done!')
      write(1,341)
341   format(1x,//,72('-'))

      goto 101
99    Write(*,*)'Unknown file error 100!'
C
101   continue
      END

C     ******************************************************************
C     
      SUBROUTINE INPT
C
C     ******************************************************************

      character*60 ZHED
C

      INTEGER IBUST
      INTEGER ICLOSE
      INTEGER NCYC
      INTEGER INEL
      INTEGER JNEL
      INTEGER NLC
      INTEGER N_DUM
      INTEGER NELS
      REAL XFAC
      INTEGER NEQ
      INTEGER NJTS
      Integer, parameter :: ii = 50

      REAL X(140),Y(140),RRTOT(140,3)

      REAL YMP(ii),YMN(ii),ROTCP(ii),ROTCN(ii)

      REAL PP(420),RR(420)
      INTEGER NA(420)
      REAL R(420)
      REAL AA(ii),E(ii),XI(ii),SKII(ii),SKJJ(ii),SKIJ(ii),YBAR(ii)
      INTEGER NSFSCI(ii),NSFSCJ(ii)
      REAL SPA(50),SPB(50),SPR(50),COSS(50),SINS(50)
      INTEGER NSPJ(50)
      REAL SPFRA(50),SPFRB(50),SPFRR(50)
      REAL EK(3,3),EKA(3,6),FK(3,3),FKA(3,6),AKA(6,6)
      INTEGER LM(6)
      REAL A(3,6),AB(3,6)
      INTEGER IHI(200),IHJ(200),NODI(200),NODJ(200),NSFTP(200)
      INTEGER NSECI(200),NSECJ(200)
      REAL STOT(200,3),VPTOT(200,3),VPTP(200,3),VPTN(200,3)
      REAL S(200,3),VP(200,3)

      INTEGER NSJTS,NTERM
      REAL FAX,DMAX
      INTEGER NLCS,MBAND
      INTEGER NUMCP,NUMNG,NUMSFT,NUMCST,NUMEG,NUMSFG,NUMSTG,N,I
      INTEGER NDI,NDJ,NN,ND
      REAL PROP
      INTEGER M,NG,NFE,NFM,NTP,NLM,MD,J,NTPI,NTPJ

      COMMON /CONTR/ IBUST,ICLOSE,NCYC,INEL,JNEL,NLC,XFAC,NELS,NEQ,NJTS,
     1NSJTS,DMAX,NTERM,FAX,NLCS,MBAND
      COMMON /JOINT/ X,Y,RRTOT
      COMMON /SECTAB/ YMP,YMN,ROTCP,ROTCN
      COMMON /EQNS/ PP,RR,NA,R
      COMMON /STFTAB/ AA,E,XI,SKII,SKJJ,SKIJ,NSFSCI,NSFSCJ
      COMMON /SUPJT/ SPA,SPB,SPR,COSS,SINS,NSPJ,SPFRA,SPFRB,SPFRR
      COMMON /STIFF/ EK,EKA,FK,FKA,AKA,LM,A,AB
      COMMON /MEMBER/ IHI,IHJ,NODI,NODJ,NSFTP,NSECI,NSECJ,
     2STOT,VPTOT,VPTP,VPTN,S,VP
      Common /MODUL/ YBAR

C
C     READ AND PRINT CONTROL INFORMATION
C
      READ(2,5) zHED
    5 FORMAT(A60)

c     write(1,15)zHED
c  15 FORMAT(A60)

      READ(2,*) NJTS,NUMCP,NUMNG,NSJTS,NUMSFT,NUMCST,NELS,NUMEG,NUMSF
     1G,NUMSTG,NLCS,NTERM

   10 FORMAT(12I5)
      IF(NJTS.EQ.0)STOP
      write(1,20) zhed,NJTS,NUMCP,NUMNG,NSJTS,NUMSFT,NUMCST,NELS,NUMEG,
     1numsFG,NUMSTG,NLCS,NTERM
   20 FORMAT(1x,72('-'),/,4X,'ULARC STRUCTURAL PROGRAM',/
     11x,'   ELASTO-PLASTIC ANALYSIS OF PLANE FRAMES'2X,
     21x,'BERKELEY'/1x ,72('-')//,10x,A60,//1x,72('-'),/
     $1x,'   AMENDED J. NICHOLS, 2015. Complied Intel Fortran',/
     11x,72('-'),//
     $1x,' Control Data for the Model'//
     31x,' NUMBER OF NODES                 ='I5/
     41x,' NUMBER OF CONTROL POINTS        ='I5/
     51x,' NUMBER OF NODE GENERATIONS      ='I5//
     61x,' NUMBER OF SUPPORTED NODES       ='I5//
     71x,' NUMBER OF STIFFNESS TYPES       ='I5/
     81x,' NUMBER OF SECTION TYPES         ='I5//
     91x,' NUMBER OF MEMBERS               ='I5/
     $1x,' NUMBER OF MEMBER GENERATIONS    ='I5/
     $1x,' NUMBER OF STIFFNESS GENERATIONS ='I5/
     $1x,' NUMBER OF STRENGTH GENERATIONS  ='I5///
     $1x,' NUMBER OF LOAD CASES            ='I5//
     $1x,' MAXIMUM NUMBER OF CYCLES        ='I5)
C
C     INITIALIZE NODE COORDINATE ARRAY
C
      DO 30 N=1,NJTS
      X(N)=999999.99
   30 Y(N)=999999.99
C
C     CONTROL NODE COORDINATES
C
      write(1,40)
   40 FORMAT(/////,' CONTROL NODE COORDINATES',//
     1'  NODE        X-ORD        Y-ORD',/)
      DO 50 N=1,NUMCP
      READ(2,* ) I,X(I),Y(I)
   50 write(1,70) I,X(I),Y(I)
   60 FORMAT(I5,2F10.0)
   70 FORMAT(I4,1X,2F14.3)
C
C     EXECUTE GENERATION COMMANDS
C
      IF(NUMNG.EQ.0) GO TO 120
      write(1,80)
   80 FORMAT( ///4x,' COORDINATE GENERATION COMMANDS'//
     12x,' FIRST',3X,' LAST',3X,'NO. OF','    NODE','   DISTANCE'/
     2'    NODE',3X,' NODE',4X,'NODES','    DIFF'/)
      DO 110 N=1,NUMNG
          READ(2,*) NDI,NDJ,NN,ND,PROP
          write(1,100) NDI,ndj,nn,nd,prop
   90     FORMAT(4I5,F10.0)
  100     FORMAT(4i8,f10.0 )
          CALL LINGEN (X,Y,NDI,NDJ,NN,ND,PROP)
  110 CONTINUE
C
C     GENERATE UNSPECIFIED NODES
C
  120 N=1
  130 N=N+1
      IF (N.GE.NJTS) GO TO 150
      IF (X(N).NE.999999.99) GO TO 130
      NDJ=N
      NDI=N-1
  140 NDJ=NDJ+1
      IF (NDJ.GT.NJTS) GO TO 150
      IF (X(NDJ).EQ.999999.99) GO TO 140
      NN=NDJ-NDI-1
      PROP=0.
      CALL LINGEN (X,Y,NDI,NDJ,NN,1,FLOAT(NJTS))
      N=NDJ
      GO TO 130
  150 CONTINUE

C
C     COMPLETE NODE COORDINATES
C

      write (1,160) (N,X(N),Y(N),N=1,NJTS)
  160 FORMAT(/////' COMPLETE NODE COORDINATES',//
     1' NODE         X-ORD         Y-ORD',//
     2(I4,1X,2F14.3))

C
C     SUPPORT SPRING DATA
C

      write (1,170)
  170 FORMAT(/////' SUPPORTED NODES',//
     1' NODE    A-STIFFNESS    B-STIFFNESS',4X,
     2'R-STIFFNESS        ANGLE',/)

      READ (2,*) (NSPJ(M),SPA(M),SPB(M),SPR(M),COSS(M),M=1,NSJTS)
      write (1,190) (NSPJ(M),SPA(M),SPB(M),SPR(M),COSS(M),M=1,NSJTS)

  180 FORMAT(I5,4E10.3)
  190 FORMAT(I4,1X,3E15.5,F13.2)

      DO 200 N=1,NSJTS
          SINS(N)=SIN(COSS(N)*.01745329)
  200 COSS(N)=SQRT(1.-SINS(N)**2)

C
C     MEMBER STIFFNESS TYPE TABLE
C

      WRITE(1,215)

      DO 211 N=1,NUMSFT
          Read (2,*)N_DUM, AA(N), XI(N), YBAR(N), E(N),SKII(N),SKJJ(N),
     1    SKIJ(N),
     1    nsfsci(n),NSFSCJ(N)
      write (1,220)N,AA(N),XI(N),E(N),YBAR(N),SKII(N),SKJJ(N),SKIJ(N),
     1    nsfsci(n),nSFSCJ(N)
  211 CONTINUE

  210 FORMAT(I5,6F10.0,2I5)
  215 FORMAT(/////' STIFFNESS TYPE TABLE'//
     1' STIFF.',3X,'SECTION',2X,'MOMENT OF',3X,' YOUNGS    Y',
     26X,' STIFFNESS COEFFICIENTS' ,3X,'STRENGTH TYPES'/
     3' TYPE' ,3X,'    AREA' ,2X,'  INERTIA' ,2X,'   MODULUS',
     42X,'BAR',6X,'KII',5X,'KJJ',5X,'KIJ',5X,' END I  END J'//)
  220 FORMAT(I4,2X,2E11.3,E11.3,1X,F5.3,2X,3F8.3,3X,I4,2x,I6)

C
C     SECTION STRENGTH TYPE TABLE
C
      READ (2,*) (N,YMP(N),YMN(N),ROTCP(N),ROTCN(N),M=1,NUMCST)
      write (1,240) (N,YMP(N),YMN(N),ROTCP(N),ROTCN(N),N=1,NUMCST)
  230 FORMAT(I5,4F10.0)
  240 FORMAT(/////' SECTION STRENGTH TYPE TABLE',//
     1' SECTION           YIELD MOMENTS    ',
     2'         HINGE ROTATION LIMITS'/
     3'   TYPE        POSITIVE     NEGATIVE',8X,
     4'POSITIVE     NEGATIVE'//
     5(I6,5X,2F13.2,3X,2F13.5))
      DO 250 N=1,NUMCST
          IF (ROTCP(N).EQ.0.) ROTCP(N)=1.E10
          IF (ROTCN(N).EQ.0.) ROTCN(N)=1.E10
  250 ROTCN(N)=-ABS(ROTCN(N))

C
C     MEMBER LOCATION GENERATION
C

      write (1,260)
  260 FORMAT(/////' MEMBER LOCATION GENERATION',//
     1' FIRST     NODE    NODE    LAST     NODE',/
     2' MEMBER     I       J     MEMBER    DIFF',/)

      DO 300 NG=1,NUMEG
          READ (2,*) NFE,NDI,NDJ,NN,ND
          write (1,280) NFE,NDI,NDJ,NN,ND

  270     FORMAT(5I5)
  280     FORMAT(I5,I9,I8,I8,I9)
C
          IF (NN.EQ.0) NN=NFE
          NN=NN-NFE+1
          DO 290 N=1,NN
              NODI(NFE)=NDI
              NODJ(NFE)=NDJ
              NFE=NFE+1
              NDI=NDI+ND
              NDJ=NDJ+ND
  290     CONTINUE
  300 CONTINUE
C
C     MEMBER STIFFNESS TYPE GENERATION
C
      write (1,310)
  310 FORMAT(/////26H STIFFNESS TYPE GENERATION//
     17H  FIRST,5X,9HSTIFFNESS,5X,7H  LAST ,5X,6HMEMBER/
     27H MEMBER,5X,9H   TYPE  ,5X,7H MEMBER,5X,6H DIFF /)
      DO 360 I=1,NUMSFG
          READ (2,*) NFM,NTP,NLM,MD
          write (1,330) NFM,NTP,NLM,MD
  320     FORMAT(4I5)
  330     FORMAT(I5,7X,I6,8X,I5,7X,I4)
C
          IF (NLM.EQ.0) GO TO 340
          NLM=(NLM-NFM)/MD+1
          GO TO 350
  340     NLM=1
  350     DO 360 N=1,NLM
              NSFTP(NFM)=NTP
              NSECI(NFM)=NSFSCI(NTP)
              NSECJ(NFM)=NSFSCJ(NTP)
  360 NFM=NFM+MD
C
C     MEMBER STRENGTH GENERATION
C
      IF (NUMSTG.EQ.0) GO TO 430
      write (1,370)
  370 FORMAT(/////25H STRENGTH TYPE GENERATION//
     17H  FIRST,5X,8HSTRENGTH,5X,8HSTRENGTH,5X,7H  LAST ,5X,6HMEMBER/
     27H MEMBER,5X,8H TYPE I ,5X,8H TYPE J ,5X,7H MEMBER,5X,5H DIFF/)
      DO 420 J=1,NUMSTG
          REad (2,*) NFM,NTPI,NTPJ,NLM,MD
          write (1,390) NFM,NTPI,NTPJ,NLM,MD
  380     FORMAT(5I5)
  390     FORMAT(I5,7X,I6,7X,I6,7X,I5,7X,I4)
C
          IF (NLM.EQ.0) GO TO 400
          NLM=(NLM-NFM)/MD+1
          GO TO 410
  400     NLM=1
  410     DO 420 N=1,NLM
              NSECI(NFM)=NTPI
              NSECJ(NFM)=NTPJ
  420 NFM=NFM+MD
C
C     COMPLETE MEMBER TABLE
C
  430 write (1,440) (N,NODI(N),NODJ(N),NSFTP(N),NSECI(N),NSECJ(N),N=1,
     1NELS)
  440 FORMAT(/////21H COMPLETE MEMBER DATA//
     17H MEMBER,4X,4HNODE,4X,4HNODE,3X,9HSTIFFNESS,3X,7HSECTION,
     23X,7HSECTION/
     37H   NO. ,4X,4H  I ,4X,'  J ',3X,9H   TYPE  ,3X,7H   I   ,
     43X,7H   J	  //
     5(I5,6X,I3,5X,I3,4X,I6,6X,I4,6X,I4))
C
      RETURN
      END
      SUBROUTINE LINGEN (X,Y,NDI,NDJ,NN,ND,PROP)
C
C     LINEAR GENERATION OF NODE COORDINATES
C
      REAL X(1),Y(1)
      REAL XI,YI,DX,DY
      REAL PROP
      INTEGER N
C
      XI=X(NDI)
      YI=Y(NDI)
      IF (XI.EQ.999999.99) write (1,10) NDI
      IF (X(NDJ).EQ.999999.99) write (1,10) NDJ
   10 FORMAT(5H NODE,I4,16H NOT YET DEFINED)
      DX=X(NDJ)-XI
      DY=Y(NDJ)-YI
      IF (PROP.LT.1) GO TO 20
      PROP=PROP/SQRT(DX**2+DY**2)
   20 IF (PROP.EQ.0.) PROP=1./FLOAT(NN+1)
      DX=DX*PROP
      DY=DY*PROP
C
      DO 30 N=1,NN
          NDI=NDI+ND
          XI=XI+DX
          YI=YI+DY
          X(NDI)=XI
   30 Y(NDI)=YI
C
      RETURN
      END

C     ******************************************************************
C      
      SUBROUTINE BANDC
C
C     ******************************************************************      
C
      INTEGER IBUST
      INTEGER ICLOSE
      INTEGER NCYC
      INTEGER INEL
      INTEGER JNEL
      INTEGER NLC
      INTEGER NELS
      REAL XFAC
      INTEGER NEQ
      INTEGER NJTS

      INTEGER NSJTS,NTERM
      REAL FAX,DMAX
      INTEGER NLCS,MBAND
      INTEGER J,NE,NI,NJ
      INTEGER L,KK,JJ,K,II
      INTEGER NSUM

C      REAL YMP(20),YMN(20)
C     REAL       ROTCP(20)

      REAL PP(420),RR(420)
      INTEGER NA(420)
      REAL R(420)
C      REAL XI(15),SKII(15),SKJJ(15),SKIJ(15)
C      INTEGER NSFSCJ(15)
C      REAL SPA(50),SPB(50),SPR(50)
C       REAL COSS(50),SINS(50)
C      INTEGER NSPJ(50)
C      REAL SPFRA(50),SPFRB(50)
C      REAL SPFRR(50)
      REAL EK(3,3),EKA(3,6),FK(3,3),FKA(3,6),AKA(6,6)
      INTEGER LM(6)
      REAL A(3,6),AB(3,6)
      INTEGER IHI(200),IHJ(200),NODI(200),NODJ(200),NSFTP(200)
      INTEGER NSECI(200),NSECJ(200)
      REAL STOT(200,3),VPTOT(200,3),VPTP(200,3),VPTN(200,3)
      REAL S(200,3),VP(200,3)

C      REAL SK(6000)

      COMMON /CONTR/ IBUST,ICLOSE,NCYC,INEL,JNEL,NLC,XFAC,NELS,NEQ,NJTS,
     1NSJTS,DMAX,NTERM,FAX,NLCS,MBAND
      COMMON /EQNS/ PP,RR,NA,R
      COMMON /STIFF/ EK,EKA,FK,FKA,AKA,LM,A,AB
      COMMON /MEMBER/ IHI,IHJ,NODI,NODJ,NSFTP,NSECI,NSECJ,
     2STOT,VPTOT,VPTP,VPTN,S,VP
C
C     COMPUTATION OF COLUMN BAND WIDTH
C
      DO 10 J=2,NEQ
   10 NA(J)=J-1
      NA(1)=1
      DO 40 NE=1,NELS
          NI=NODI(NE)
          NJ=NODJ(NE)
          LM(6)=3*NJ
          LM(5)=LM(6)-1
          LM(4)=LM(5)-1
          LM(3)=3*NI
          LM(2)=LM(3)-1
          LM(1)=LM(2)-1
          DO 30 L=1,6
              KK=NEQ
              JJ=LM(L)
              DO 20 K=1,6
                  II=LM(K)
   20         KK=MIN0(KK,II)
              IF (KK.GT.NA(JJ)) GO TO 30
              NA(JJ)=KK
   30     CONTINUE
   40 CONTINUE
      NSUM=1
      DO 50 J=2,NEQ
          NSUM=NSUM+J-NA(J)+1
   50 NA(J)=NSUM
C
      RETURN
      END

C-------------------------------------------------------------------------
C
C     Loading Routine
C
C-------------------------------------------------------------------------

      SUBROUTINE LOAD
      implicit None
      character*60 ZHED

      INTEGER IBUST
      INTEGER ICLOSE
      INTEGER NCYC
      INTEGER INEL
      INTEGER JNEL
      INTEGER NLC
      INTEGER NELS
      REAL XFAC
      INTEGER NEQ
      INTEGER NJTS
      INTEGER NSJTS
      REAL DMAX
      INTEGER NTERM
      REAL FAX
      INTEGER NLCS,MBAND
      Integer, parameter :: ii = 50

      REAL X(140),Y(140),RRTOT(140,3)

      REAL YMP(ii),YMN(ii),ROTCP(ii),ROTCN(ii)

      REAL PP(420),RR(420)
      INTEGER NA(420)
      REAL R(420)
      REAL AA(ii),E(ii),XI(ii),SKII(ii),SKJJ(ii),SKIJ(ii)
      INTEGER NSFSCI(ii),NSFSCJ(ii)
      REAL SPA(50),SPB(50),SPR(50),COSS(50),SINS(50)
      INTEGER NSPJ(50)
      REAL SPFRA(50),SPFRB(50),SPFRR(50)
      REAL EK(3,3),EKA(3,6),FK(3,3),FKA(3,6),AKA(6,6)
      INTEGER LM(6)
      REAL A(3,6),AB(3,6)
      INTEGER IHI(200),IHJ(200),NODI(200),NODJ(200),NSFTP(200)
      INTEGER NSECI(200),NSECJ(200)
      REAL STOT(200,3),VPTOT(200,3),VPTP(200,3),VPTN(200,3)
      REAL S(200,3),VP(200,3)

      INTEGER I,NUMLG,NUMDG,LCODE,NE,M,NJ,NL,NFJ
      REAL PX,PY,PR
      INTEGER NLJ,INCR,N,NR,ND,NS,J
      REAL CC,SS,SFA,SFB,SFR,SFX,SFY


      COMMON /CONTR/ IBUST,ICLOSE,NCYC,INEL,JNEL,NLC,XFAC,NELS,NEQ,NJTS,
     1NSJTS,DMAX,NTERM,FAX,NLCS,MBAND
      COMMON /JOINT/ X,Y,RRTOT
      COMMON /SECTAB/ YMP,YMN,ROTCP,ROTCN
      COMMON /EQNS/ PP,RR,NA,R
      COMMON /STFTAB/ AA,E,XI,SKII,SKJJ,SKIJ,NSFSCI,NSFSCJ
      COMMON /SUPJT/ SPA,SPB,SPR,COSS,SINS,NSPJ,SPFRA,SPFRB,SPFRR
      COMMON /STIFF/ EK,EKA,FK,FKA,AKA,LM,A,AB
      COMMON /MEMBER/ IHI,IHJ,NODI,NODJ,NSFTP,NSECI,NSECJ,
     2STOT,VPTOT,VPTP,VPTN,S,VP
C     
C
C
C     INITIALIZE ARRAYS
C
      read(2,8)zhed
      write(*,8)zhed
8     format(1x,a60)
      DO 10 I=1,NEQ
          R(I)=0.
   10 PP(I)=0.
C
C     READ CONTROL INFORMATION
C
      READ (2,*) NUMLG,NUMDG,LCODE,FAX,DMAX
   20 FORMAT(3I5,2F10.0,11A4,A1)
      write (1,30) NLC,zhed,NUMLG,NUMDG,LCODE,FAX,DMAX
   30 FORMAT(1x///' LOAD CASE',I3,' ',A60/
     1' NO. OF NODAL LOAD GENERATIONS      =  'I4/
     2' NO. OF SUPPORT DISPL GENERATIONS   =  'I4/
     3' LOAD APPLICATION CODE              =  'I4/
     4' MAXIMUM LOAD FACTOR                ='F10.3/
     5' LIMITING DISPLACEMENT              ='F10.3 )
C
C     INITIALIZE IF NEW LOAD SEQUENCE
C
      IF (NLC.EQ.1) GO TO 40
      IF (LCODE.EQ.0) GO TO 90
   40 DO 50 NE=1,NELS
      IHI(NE)=0
   50 IHJ(NE)=0
      IBUST=0
      DO 60 M=1,NSJTS
          SPFRA(M)=0.
          SPFRB(M)=0.
   60 SPFRR(M)=0.
      DO 70 NJ=1,NJTS
          DO 70 I=1,3
   70 RRTOT(NJ,I)=0.
      DO 80 NE=1,NELS
          DO 80 I=1,3
              VPTOT(NE,I)=0.
              STOT(NE,I)=0.
              VPTP(NE,I)=0.
   80 VPTN(NE,I)=0.

C
C     READ AND PRINT LOAD GENERATION COMMANDS.
C

   90 IF(NUMLG.EQ.0) GO TO 200
      write (1,100)
  100 FORMAT(////' NODAL LOAD GENERATION',//
     1' FIRST       X-LOAD       Y-LOAD       MOMENT',
     2'     LAST      NODE',/
     3' NODE',45X,'NODE      DIFF',/)
      DO 160 NL=1,NUMLG
          READ (2,*) NFJ,PX,PY,PR,NLJ,INCR

          write (1,120) NFJ,PX,PY,PR,NLJ,INCR
  110     FORMAT(I5,3F10.0,2I5)
  120     FORMAT(I5,1X,3F13.3,2I8)
          IF (NLJ.EQ.0) GO TO 130
          NLJ=(NLJ-NFJ)/INCR+1
          GO TO 140
  130     NLJ=1
  140     DO 150 N=1,NLJ
              NR=3*NFJ-2
              PP(NR)=PP(NR)+PX
              PP(NR+1)=PP(NR+1)+PY
              PP(NR+2)=PP(NR+2)+PR
  150     NFJ=NFJ+INCR
  160 CONTINUE
C
C     PRINT MATRIX OF APPLIED LOADS
C
      write (1,180)
      DO 170 I=1,NJTS
          NR=3*I-2
  170 write (1,190) I,PP(NR),PP(NR+1),PP(NR+2)
  180 FORMAT(////' NODAL LOADS',//
     1' NODE       X-LOAD       Y-LOAD       MOMENT',/)
  190 FORMAT(I4,1X,3F13.3)
C
C     READ AND PRINT DISPLACEMENT GENERATION
C
  200 IF (NUMDG.EQ.0) GO TO 320
      write (1,210)
  210 FORMAT(////32H SUPPORT DISPLACEMENT GENERATION//
     16H FIRST,6X,7HX-DISPL,6X,7HY-DISPL,5X,8HROTATION,
     25X,4HLAST,5X,4HNODE/
     36H  NODE,39X,5X,4HNODE,5X,4HDIFF/)
      DO 270 ND=1,NUMDG
          READ (2,220) NFJ,PX,PY,PR,NLJ,INCR
          write (1,230) NFJ,PX,PY,PR,NLJ,INCR
  220     FORMAT(I5,3F10.0,2I5)
  230     FORMAT(I5,1X,3F13.4,2I8)
          IF (NLJ.EQ.0) GO TO 240
          NLJ=(NLJ-NFJ)/INCR+1
          GO TO 250
  240     NLJ=1
  250     DO 260 N=1,NLJ
              NR=3*NFJ-2
              R(NR)=R(NR)+PX
              R(NR+1)=R(NR+1)+PY
              R(NR+2)=R(NR+2)+PR
  260     NFJ=NFJ+INCR
  270 CONTINUE
C
C     PRINT APPLIED DISPLACEMENTS
C
      write (1,280)
  280 FORMAT(////22H SUPPORT DISPLACEMENTS//
     15H NODE,6X,7HX-DISPL,6X,7HY-DISPL,5X,8HROTATION/)
      DO 290 I=1,NSJTS
          NS=NSPJ(I)
          NR=3*NS-2
  290 write (1,300) NS,R(NR),R(NR+1),R(NR+2)
  300 FORMAT(I4,1X,3F13.4)
C
C     STORE SPRING PRELOADS, SET UP NODAL LOADS
C
      DO 310 J=1,NSJTS
          NS=NSPJ(J)
          NR=3*NS-2
          CC=COSS(J)
          SS=SINS(J)
          SFA=(R(NR)*CC+R(NR+1)*SS)*SPA(J)
          SFB=(-R(NR)*SS+R(NR+1)*CC)*SPB(J)
          SFR=R(NR+2)*SPR(J)
          SFX=SFA*CC-SFB*SS
          SFY=SFA*SS+SFB*CC
          PP(NR)=PP(NR)+SFX
          PP(NR+1)=PP(NR+1)+SFY
          PP(NR+2)=PP(NR+2)+SFR
          R(NR)=SFA
          R(NR+1)=SFB
          R(NR+2)=SFR
  310 CONTINUE
C
  320 IF (IBUST.EQ.1) write (1,330)
  330 FORMAT(/////45H STRUCTURE COLLAPSED.  THIS LOAD CASE IGNORED)
      IF (IBUST.EQ.2) write (1,340)
  340 FORMAT(/////45H CYCLE LIMIT EXCEEDED. THIS LOAD CASE IGNORED)
C
      RETURN
      END

C------------------------------------------------------------------------
C
C
C
C------------------------------------------------------------------------

      SUBROUTINE STIF (SK)


      INTEGER IBUST
      INTEGER ICLOSE
      INTEGER NCYC
      INTEGER INEL
      INTEGER JNEL
      INTEGER NLC
      INTEGER NELS
      REAL XFAC
      INTEGER NEQ
      INTEGER NJTS
      INTEGER NSJTS,NTERM
      REAL FAX,DMAX
      INTEGER NLCS,MBAND
      Integer, parameter :: iii = 50

      REAL X(140),Y(140),RRTOT(140,3)

      REAL YMP(iii),YMN(iii),ROTCP(iii),ROTCN(iii)

      REAL PP(420),RR(420)
      INTEGER NA(420)
      REAL R(420)
      REAL AA(iii),E(iii),XI(iii),SKII(iii),SKJJ(iii),SKIJ(iii)
      INTEGER NSFSCI(iii),NSFSCJ(iii)
      REAL SPA(50),SPB(50),SPR(50),COSS(50),SINS(50)
      INTEGER NSPJ(50)
      REAL SPFRA(50),SPFRB(50),SPFRR(50)
      REAL EK(3,3),EKA(3,6),FK(3,3),FKA(3,6),AKA(6,6)
      INTEGER LM(6)
      REAL A(3,6),AB(3,6)
      INTEGER IHI(200),IHJ(200),NODI(200),NODJ(200),NSFTP(200)
      INTEGER NSECI(200),NSECJ(200)
      REAL STOT(200,3),VPTOT(200,3),VPTP(200,3),VPTN(200,3)
      REAL S(200,3),VP(200,3)

      INTEGER JJ,I,NE,NI,NJ
      REAL DX,DY,ELL,COSA,SINA,XM,EKAA
      INTEGER J,K
      REAL FKAA
      INTEGER II,JJJ,N,NS
      REAL CS,SS
      INTEGER NR
      REAL SK(1)

      COMMON /CONTR/ IBUST,ICLOSE,NCYC,INEL,JNEL,NLC,XFAC,NELS,NEQ,NJTS,
     1NSJTS,DMAX,NTERM,FAX,NLCS,MBAND
      COMMON /JOINT/ X,Y,RRTOT
      COMMON /SECTAB/ YMP,YMN,ROTCP,ROTCN
      COMMON /EQNS/ PP,RR,NA,R
      COMMON /STFTAB/ AA,E,XI,SKII,SKJJ,SKIJ,NSFSCI,NSFSCJ
      COMMON /SUPJT/ SPA,SPB,SPR,COSS,SINS,NSPJ,SPFRA,SPFRB,SPFRR
      COMMON /STIFF/ EK,EKA,FK,FKA,AKA,LM,A,AB
      COMMON /MEMBER/ IHI,IHJ,NODI,NODJ,NSFTP,NSECI,NSECJ,
     2STOT,VPTOT,VPTP,VPTN,S,VP
C
C
C     INITIALIZE STRUCTURE STIFFNESS MATRIX
C
      JJ=NA(NEQ)
      DO 10 I=1,JJ
   10 SK(I)=0.0
C
C     SET UP ELEMENT STIFFNESS
C
      DO 80 NE=1,NELS
C
C     COMPUTE TERMS TO SET UP 3*3 ELEMENT STIFFNESS MATRIX
C
                      I=NSFTP(NE)
                      NI=NODI(NE)
                      NJ=NODJ(NE)
                      DX=X(NJ)-X(NI)
                      DY=Y(NJ)-Y(NI)
                      ELL=SQRT(DX**2+DY**2)
                      COSA=DX/ELL
                      SINA=DY/ELL
C
C     FORM ELEMENT STIFFNESS MATRIX IN LOCAL COORDINATES
C
                      XM=E(I)*XI(I)/ELL
                      EK(1,1)=E(I)*AA(I)/ELL
                      EK(1,2)=0.
                      EK(1,3)=0.
                      EK(2,2)=XM*SKII(I)
                      EK(2,3)=XM*SKIJ(I)
                      EK(3,3)=XM*SKJJ(I)
                      FK(1,1)=EK(1,1)*1.E-8
                      FK(1,2)=0.
                      FK(1,3)=0.
                      FK(2,2)=XM*4.*1.E-8
                      FK(2,3)=FK(2,2)*0.5
                      FK(3,3)=FK(2,2)
                      DO 20 I=1,3
                          DO 20 J=1,3
                              EK(J,I)=EK(I,J)
   20                 FK(J,I)=FK(I,J)
C
C     FORM ELEMENT STIFFNESS MATRIX IN GLOBAL COORDINATES
C
          CALL FORMA (A,IHI(NE),IHJ(NE),COSA,SINA,ELL)
          CALL FORMA (AB,0,0,COSA,SINA,ELL)
          DO 40 I=1,3
              DO 40 J=1,6
                  EKAA=0.
                  FKAA=0.
                  DO 30 K=1,3
                      FKAA=FKAA+FK(I,K)*AB(K,J)
   30             EKAA=EKAA+EK(I,K)*A(K,J)
                  EKA(I,J)=EKAA
   40     FKA(I,J)=FKAA
C
C     ADD SMALL ELASTIC STIFFNESS TO ELEMENT STIFFNESS
C
          DO 60 I=1,6
              DO 60 J=1,6
                  EKAA=0.
                  DO 50 K=1,3
   50             EKAA=EKAA+A(K,I)*EKA(K,J)+AB(K,I)*FKA(K,J)
   60     AKA(I,J)=EKAA
C
C     ADD ELEMENT STIFFNESS TO STRUCTURE STIFFNESS MATRIX
C
          LM(6)=3*NJ
          LM(5)=LM(6)-1
          LM(4)=LM(5)-1
          LM(3)=3*NI
          LM(2)=LM(3)-1
          LM(1)=LM(2)-1
          DO 70 I=1,6
              II=LM(I)
              DO 70 J=1,6
                  JJ=LM(J)
                  IF (JJ.LT.II) GO TO 70
                  JJJ=NA(JJ)-(JJ-II)
                  SK(JJJ)=SK(JJJ)+AKA(I,J)
   70     CONTINUE
C
   80 CONTINUE
C
C     ADD SUPPORT STIFFNESS
C
      DO 90 N=1,NSJTS
          NS=NSPJ(N)
          CS=COSS(N)
          SS=SINS(N)
          NR=3*NS-2
          JJ=NA(NR)
          SK(JJ)=SK(JJ)+SPA(N)*(CS**2)+SPB(N)*(SS**2)
          JJ=NA(NR+1)
          SK(JJ)=SK(JJ)+SPA(N)*(SS**2)+SPB(N)*(CS**2)
          JJ=JJ-1
          SK(JJ)=SK(JJ)+CS*SS*(SPA(N)-SPB(N))
          JJ=NA(NR+2)
          SK(JJ)=SK(JJ)+SPR(N)
   90 CONTINUE
C
      RETURN
      END

C--------------------------------------------------------------------------
C
C
C
C--------------------------------------------------------------------------

      SUBROUTINE FORCE
C

      INTEGER IBUST
      INTEGER ICLOSE
      INTEGER NCYC
      INTEGER INEL
      INTEGER JNEL
      INTEGER NLC
      INTEGER NELS
      REAL XFAC
      INTEGER NEQ
      INTEGER NJTS
      INTEGER NSJTS,NTERM
      REAL FAX,DMAX
      INTEGER NLCS,MBAND
      Integer, parameter :: ii = 50

      REAL X(140),Y(140),RRTOT(140,3)

      REAL YMP(ii),YMN(ii),ROTCP(ii),ROTCN(ii)

      REAL PP(420),RR(420)
      INTEGER NA(420)
      REAL R(420)
      REAL AA(ii),E(ii),XI(ii),SKII(ii),SKJJ(ii),SKIJ(ii)
      INTEGER NSFSCI(ii),NSFSCJ(ii)
      REAL SPA(50),SPB(50),SPR(50),COSS(50),SINS(50)
      INTEGER NSPJ(50)
      REAL SPFRA(50),SPFRB(50),SPFRR(50)
      REAL EK(3,3),EKA(3,6),FK(3,3),FKA(3,6),AKA(6,6)
      INTEGER LM(6)
      REAL A(3,6),AB(3,6)
      INTEGER IHI(200),IHJ(200),NODI(200),NODJ(200),NSFTP(200)
      INTEGER NSECI(200),NSECJ(200)
      REAL STOT(200,3),VPTOT(200,3),VPTP(200,3),VPTN(200,3)
      REAL S(200,3),VP(200,3)

      INTEGER NE,N,NI,NJ
      REAL DX,DY,ELL,COSA,SINA
      INTEGER I,NR,J
      REAL XM,RATIOI,RATIOJ
      INTEGER NS
      REAL DELS,RAT,RATIO,DIFFR
      INTEGER LIMJT
      REAL DELR,XXFAC
      INTEGER M
      REAL CC,SS,DSA,DSB,RMAX,RM

      REAL V(3),VE(3)
C
      COMMON /CONTR/ IBUST,ICLOSE,NCYC,INEL,JNEL,NLC,XFAC,NELS,NEQ,NJTS,
     1NSJTS,DMAX,NTERM,FAX,NLCS,MBAND
      COMMON /JOINT/ X,Y,RRTOT
      COMMON /SECTAB/ YMP,YMN,ROTCP,ROTCN
      COMMON /EQNS/ PP,RR,NA,R
      COMMON /STFTAB/ AA,E,XI,SKII,SKJJ,SKIJ,NSFSCI,NSFSCJ
      COMMON /SUPJT/ SPA,SPB,SPR,COSS,SINS,NSPJ,SPFRA,SPFRB,SPFRR
      COMMON /STIFF/ EK,EKA,FK,FKA,AKA,LM,A,AB
      COMMON /MEMBER/ IHI,IHJ,NODI,NODJ,NSFTP,NSECI,NSECJ,
     2STOT,VPTOT,VPTP,VPTN,S,VP

      write (1,10) NLC,NCYC
   10 FORMAT(/////////'  RESULTS - LOAD CASE :',I3,'  CYCLE :',I4)
C
C     MEMBER DEFORMATION INCREMENTS IN LOCAL COORDINATES
C
      DO 80 NE=1,NELS
C
C     COMPUTE 3*6 DISPLACEMENT TRANSFORMATION MATRIX
C
                      N=NSFTP(NE)
                      NI=NODI(NE)
                      NJ=NODJ(NE)
                      DX=X(NJ)-X(NI)
                      DY=Y(NJ)-Y(NI)
                      ELL=SQRT(DX**2+DY**2)
                      COSA=DX/ELL
                      SINA=DY/ELL
                      CALL FORMA (A,IHI(NE),IHJ(NE),COSA,SINA,ELL)
                      CALL FORMA (AB,0,0,COSA,SINA,ELL)
C
C     COMPUTE PLASTIC DEFORMATIONS
C
                      DO 20 I=1,3
                          V(I)=0.
   20                 VE(I)=0.
                      DO 50 I=1,3
                          NR=3*NI-3
                          DO 30 J=1,3
                              NR=NR+1
                              VE(I)=VE(I)+AB(I,J)*RR(NR)
   30                     V(I)=V(I)+A(I,J)*RR(NR)
                          NR=3*NJ-3
                          DO 40 J=4,6
                              NR=NR+1
                              VE(I)=VE(I)+AB(I,J)*RR(NR)
   40                     V(I)=V(I)+A(I,J)*RR(NR)
   50                 VP(NE,I)=VE(I)-V(I)
C
C     COMPUTE ELEMENT FORCES IN LOCAL COORDINATES
C
          XM=E(N)*XI(N)/ELL
          S(NE,1)=E(N)*AA(N)*V(1)/ELL
          S(NE,2)=XM*(SKII(N)*V(2)+SKIJ(N)*V(3))
          S(NE,3)=XM*(SKIJ(N)*V(2)+SKJJ(N)*V(3))
C
C     CHECK PLASTIC WORK FOR CLOSED HINGES
C
          IF (STOT(NE,2)*VP(NE,2).GE.0.) GO TO 70
          ICLOSE=1
          IHI(NE)=0
          write (1,60) NE,NODI(NE)
   60     FORMAT(//21H HINGE CLOSES, MEMBER,I4,6H, NODE,I4)
   70     IF (STOT(NE,3)*VP(NE,3).GE.0.) GO TO 80
          ICLOSE=1
          IHJ(NE)=0
          write (1,60) NE,NODJ(NE)
C
   80 CONTINUE
C
C     REPEAT IF HINGES CLOSED
C
      IF (ICLOSE.EQ.1) RETURN
C
C     DETERMINE MINIMUM LOAD FACTOR INCREMENT FOR NEXT HINGE
C
      RATIOI=1.E10
      RATIOJ=1.E10
      INEL=0
      JNEL=0
      DO 160 NE=1,NELS
C
              IF (IHI(NE).GT.0) GO TO 120
              NS=NSECI(NE)
              DELS=S(NE,2)
              IF (DELS) 100,120,90
   90         RAT=(YMP(NS)-STOT(NE,2))/DELS
              GO TO 110
  100         RAT=(-YMN(NS)-STOT(NE,2))/DELS
  110         IF (RAT.GT.RATIOI) GO TO 120
              RATIOI=RAT
              INEL=NE
  120         IF (IHJ(NE).GT.0) GO TO 160
              NS=NSECJ(NE)
              DELS=S(NE,3)
              IF (DELS) 140,160,130
  130         RAT=(YMP(NS)-STOT(NE,3))/DELS
              GO TO 150
  140         RAT=(-YMN(NS)-STOT(NE,3))/DELS
  150         IF (RAT.GT.RATIOJ) GO TO 160
              RATIOJ=RAT
              JNEL=NE
  160 CONTINUE
      DIFFR=RATIOI-RATIOJ
      IF (ABS(DIFFR).GT.1.0E-8) GO TO 170
      RATIO=RATIOI
      GO TO 190
  170 IF (DIFFR.LT.0.) GO TO 180
      RATIO=RATIOJ
      INEL=0
      GO TO 190
  180 RATIO=RATIOI
      JNEL=0
C
C     CHECK LOAD FACTOR FOR LIMITING DISPLACEMENT
C
  190 IF (DMAX.EQ.0.) GO TO 240
      LIMJT=0
      DO 230 NJ=1,NJTS
          NR=3*NJ-3
          DO 230 I=1,3
              NR=NR+1
              DELR=RR(NR)
              IF (DELR) 210,230,200
  200         RAT=(DMAX-RRTOT(NJ,I))/DELR
              IF (RAT.LT.0.) RAT=0.
              IF (RAT.GT.RATIO) GO TO 230
              GO TO 220
  210         RAT=(-DMAX-RRTOT(NJ,I))/DELR
              IF (RAT.LT.0.) RAT=0.
              IF (RAT.GT.RATIO) GO TO 230
  220         RATIO=RAT
              LIMJT=NJ
              IBUST=-1
  230 CONTINUE
C
C     CHECK FOR MAXIMUM LOAD FACTOR
C
  240 XXFAC=XFAC
      XFAC=XFAC+RATIO
      IF (XFAC.LT.FAX) GO TO 250
      RATIO=FAX-XXFAC
      XFAC=FAX
      IBUST=-2
C
C     UPDATE DISPLACEMENTS, CHECK COLLAPSE
C
  250 DO 260 NJ=1,NJTS
          NR=3*NJ-3
          DO 260 I=1,3
              NR=NR+1
              RRTOT(NJ,I)=RRTOT(NJ,I)+RATIO*RR(NR)
              IF (ABS(RRTOT(NJ,I)).LT.1.E4) GO TO 260
              IBUST=1
              GO TO 370
  260 CONTINUE
C
C     UPDATE MEMBER QUANTITIES
C
      DO 310 NE=1,NELS
C
              DO 270 I=1,3
                  STOT(NE,I)=STOT(NE,I)+RATIO*S(NE,I)
                  VP(NE,I)=VP(NE,I)*RATIO
  270         VPTOT(NE,I)=VPTOT(NE,I)+VP(NE,I)
C
C     ACCUMULATE MAXIMUM ROTATIONS
C
          DO 280 I=2,3
              IF (VPTOT(NE,I).GT.VPTP(NE,I)) VPTP(NE,I)=VPTOT(NE,I)
              IF (VPTOT(NE,I).LT.VPTN(NE,I)) VPTN(NE,I)=VPTOT(NE,I)
  280     CONTINUE
C
C     SET NEW HINGE CODES
C
          IF (IBUST.NE.0) GO TO 310
          IF (NE.NE.INEL) GO TO 300
          IHI(NE)=1
          write (1,290) NE,NODI(NE)
  290     FORMAT(//20H HINGE FORMS, MEMBER,I4,6H, NODE,I4)
  300     IF (NE.NE.JNEL) GO TO 310
          IHJ(NE)=1
          write (1,290) NE,NODJ(NE)
  310 CONTINUE
  315 CONTINUE
C
C     LIMITING DISPLACEMENT OR MAXIMUM LOAD FACTOR
C
      IF (IBUST.EQ.0) GO TO 350
      IF (IBUST.EQ.-2) GO TO 330
      write (1,320) LIMJT
  320 FORMAT(//'  LIMITING DISPLACEMENT REACHED AT NODE : ',I4)
      GO TO 330
  330 write (1,340)
  340 FORMAT(//'  MAXIMUM LOAD FACTOR REACHED',/)
C
C     UPDATE SUPPORT REACTIONS
C
  350 DO 360 M=1,NSJTS
          NR=3*NSPJ(M)-2
          CC=COSS(M)
          SS=SINS(M)
          DSA=RR(NR)*CC+RR(NR+1)*SS
          DSB=-RR(NR)*SS+RR(NR+1)*CC
          SPFRA(M)=SPFRA(M)-RATIO*(SPA(M)*DSA-R(NR))
          SPFRB(M)=SPFRB(M)-RATIO*(SPB(M)*DSB-R(NR+1))
  360 SPFRR(M)=SPFRR(M)-RATIO*(SPR(M)*RR(NR+2)-R(NR+2))
      RETURN
C
C     PRINT NORMALIZED COLLAPSE MODE
C
  370 RMAX=0.
      DO 380 NJ=1,NJTS
          NR=3*NJ-3
          DO 380 I=1,2
              NR=NR+1
              RM=ABS(RR(NR))
              IF (RM.GT.RMAX) RMAX=RM
  380 CONTINUE
      DO 390 N=1,NEQ
  390 RR(N)=RR(N)/RMAX
      write (1,400) (N,RR(3*N-2),RR(3*N-1),RR(3*N),N=1,NJTS)
  400 FORMAT(////20H STRUCTURE COLLAPSES///
     125H NORMALIZED COLLAPSE MODE//
     25H NODE,7X,7HX-DISPL,7X,7HY-DISPL,6X,8HROTATION/
     3(I4,1X,3F14.4))
C
      RETURN
      END

C-----------------------------------------------------------------------
C
C     Subroutine FormA
C
C-----------------------------------------------------------------------

      SUBROUTINE FORMA (A,IHI,IHJ,COSA,SINA,ELL)
      implicit None
C
      REAL A(3,6)
      INTEGER I,J
      Integer IHI, IHJ
      Real COSA
      REAL SinA
      REal ELL
C
C     FORM TRANSFORMATION MATRIX--- NO HINGE AT I AND J END
C
      A(1,1)=-COSA
      A(1,2)=-SINA
      A(1,3)=0.
      A(1,4)=COSA
      A(1,5)=SINA
      A(1,6)=0.
      A(2,1)=SINA/ELL
      A(2,2)=-COSA/ELL
      A(2,3)=1.
      A(2,4)=-SINA/ELL
      A(2,5)=COSA/ELL
      A(2,6)=0.
      A(3,1)=SINA/ELL
      A(3,2)=-COSA/ELL
      A(3,3)=0.
      A(3,4)=-SINA/ELL
      A(3,5)=COSA/ELL
      A(3,6)=1.
      IF (IHI.EQ.0.AND.IHJ.EQ.0) GO TO 40
      IF (IHI.GT.0.AND.IHJ.GT.0) GO TO 20
      IF (IHJ.GT.0) GO TO 10
C
C      FORM TRANSFORMATION MATRIX---HINGE AT I END ONLY
C
      A(2,1)=-SINA/(2.*ELL)
      A(2,2)=COSA/(2.*ELL)
      A(2,3)=0.
      A(2,4)=SINA/(2.*ELL)
      A(2,5)=-COSA/(2.*ELL)
      A(2,6)=-.5
      GO TO 40
C
C      FORM TRANSFORMATION MATRIX---HINGE AT J END ONLY
C
   10 A(3,1)=-SINA/(2.*ELL)
      A(3,2)=COSA/(2.*ELL)
      A(3,3)=-.5
      A(3,4)=SINA/(2.*ELL)
      A(3,5)=-COSA/(2.*ELL)
      A(3,6)=0.
      GO TO 40
C
C      FORM TRANSFORMATION MATRIX---HINGE AT I AND J ENDS
C
   20 DO 30 I=1,3
          DO 30 J=1,6
   30 A(I,J)=0.
      A(1,1)=-COSA
      A(1,2)=-SINA
      A(1,4)=COSA
      A(1,5)=SINA
C
   40 RETURN
      END

C-------------------------------------------------------------------------
C
C
C
C-------------------------------------------------------------------------

      SUBROUTINE OUTPT
      implicit none
C

      INTEGER IBUST
      INTEGER ICLOSE
      INTEGER NCYC
      INTEGER INEL
      INTEGER JNEL
      INTEGER NLC
      INTEGER NELS
      REAL XFAC
      INTEGER NEQ
      INTEGER NJTS
      INTEGER NSJTS,NTERM
      REAL FAX,DMAX
      INTEGER NLCS,MBAND
      Integer, parameter :: ii = 50

      REAL X(140),Y(140),RRTOT(140,3)

      REAL YMP(ii),YMN(ii),ROTCP(ii),ROTCN(ii)

      REAL PP(420),RR(420)
      INTEGER NA(420)
      REAL R(420)
      REAL AA(ii),E(ii),XI(ii),SKII(ii),SKJJ(ii),SKIJ(ii), YBAR(ii)
      INTEGER NSFSCI(ii),NSFSCJ(ii)
      REAL SPA(50),SPB(50),SPR(50),COSS(50),SINS(50)
      INTEGER NSPJ(50)
      REAL SPFRA(50),SPFRB(50),SPFRR(50)
      REAL EK(3,3),EKA(3,6),FK(3,3),FKA(3,6),AKA(6,6)
      INTEGER LM(6)
      REAL A(3,6),AB(3,6)
      INTEGER IHI(200),IHJ(200),NODI(200),NODJ(200),NSFTP(200)
      INTEGER NSECI(200),NSECJ(200)
      REAL STOT(200,3),VPTOT(200,3),VPTP(200,3),VPTN(200,3)
      REAL S(200,3),VP(200,3)
      real million

      INTEGER NJ,I,NE,IN,JN,N
      REAL CC,SS,FFX,FFY

      COMMON /CONTR/ IBUST,ICLOSE,NCYC,INEL,JNEL,NLC,XFAC,NELS,NEQ,NJTS,
     1NSJTS,DMAX,NTERM,FAX,NLCS,MBAND
      COMMON /JOINT/ X,Y,RRTOT
      COMMON /SECTAB/ YMP,YMN,ROTCP,ROTCN
      COMMON /EQNS/ PP,RR,NA,R
      COMMON /STFTAB/ AA,E,XI,SKII,SKJJ,SKIJ,NSFSCI,NSFSCJ
      COMMON /SUPJT/ SPA,SPB,SPR,COSS,SINS,NSPJ,SPFRA,SPFRB,SPFRR
      COMMON /STIFF/ EK,EKA,FK,FKA,AKA,LM,A,AB
      COMMON /MEMBER/ IHI,IHJ,NODI,NODJ,NSFTP,NSECI,NSECJ,
     2STOT,VPTOT,VPTP,VPTN,S,VP

      Common /MODUL/ YBAR
C
C
      character*3 xno,yes,yno(4)
      DATA YES /'***'/
      DATA XNO /'   '/
      Million = 1000000.0
C
C
      write (1,10) XFAC
   10 FORMAT(//' LOAD FACTOR AT END OF CYCLE = ',F10.3)
C
C     PRINT NODE DISPLACEMENTS
C
      write (1,20) (NJ,(1000.0d0*RRTOT(NJ,I),I=1,3),NJ=1,NJTS)
   20 FORMAT(/////' NODE DISPLACEMENTS AT END OF CYCLE'//
     1' NODE',7X,'X-DISPL',7X,'Y-DISPL',6X,'ROTATION',//
     2(I4,1X,3F14.3))
     
      write (3,120) (NJ,(1000.0d0*RRTOT(NJ,I),I=1,3),NJ=1,NJTS)
120   FORMAT(/////' NODE DISPLACEMENTS AT END OF CYCLE'//
     1' NODE',7X,'X-DISPL',7X,'Y-DISPL',6X,'ROTATION',//
     2(I4,3(",",F14.3)))
C
C     PRINT MEMBER FORCES
C
      write (1,30) (NE,NODI(NE),NODJ(NE),(STOT(NE,I),I=1,3),IHI(NE),
     1ihj(ne),nseci(ne),nsecj(ne) ,NE=1,NELS)
   30 FORMAT(/////'MEMBER ACTIONS AND HINGE CODES AT END OF CYCLE',//,
     1'MEMBER   NODE   NODE       AXIAL,       MOMENT',
     2'      MOMENT           HINGE CODES'/,
     3'  NO.     I   	 J         FORCE           I     ',
     4'     J	             I     J ' //
     5(I5,I8,I7,1X,3F13.3,I10,I6,2i4))
      Write(1,32)
32    Format(///'   Stress Results'//
     1'   Member   Node   Node     Node        Node    ',
     1'    Axial       Section     Bending',/
     1'    No.      I      J       Area        Area    '
     1'    Stress      Depth       Stress',/)
      write(1,31)(NE,NODI(NE),NODJ(NE),aa(nseci(ne)),aa(nsecj(ne)),
     1((stot(ne,1)/aa(nseci(ne)))/million),YBAR(nseci(ne)),
     1((((stot(ne,2))*YBAR(nseci(ne)))/XI(nseci(ne))/million)),
     1((((stot(ne,3))*YBAR(nseci(ne)))/XI(nseci(ne))/million)),NE=1,NELS
     1)
31    Format(2x,I4,4x,I4,4x,I4,2x,F10.5,2x,F10.5,2x,F14.5,3x, F5.3, 2x,
     1F10.3,2x, f10.3)

      write(3,33)(NE,NODI(NE),NODJ(NE),aa(nseci(ne)),aa(nsecj(ne)),
     1((stot(ne,1)/aa(nseci(ne)))/million),YBAR(nseci(ne)),
     1((((stot(ne,2))*YBAR(nseci(ne)))/XI(nseci(ne))/million)),
     1((((stot(ne,3))*YBAR(nseci(ne)))/XI(nseci(ne))/million)),NE=1,NELS
     1)

33    Format(I4,",",I4,",",I4,",",F10.5,",",F10.5,",",F14.5,",", F5.3,
     1",",F10.3,",", f10.3)

C
C     PLASTIC DEFORMATIONS
C
      write (1,40)
   40 FORMAT(/////34H PLASTIC HINGE ROTATION INCREMENTS
     128H AND MAXIMUM HINGE ROTATIONS//
     27H MEMBER,3X,4HNODE,3X,4HNODE,8X,19HROTATION INCREMENTS,
     310X,22HMAX POSITIVE ROTATIONS,10X,
     422HMAX NEGATIVE ROTATIONS/
     57H   NO. ,3X,4H  I ,3X,4H	J ,8X,5HEND I,9X,5HEND J,
     610X,5HEND I,12X,5HEND J,10X,5HEND I,12X,5HEND J/)
      DO 70 NE=1,NELS
C
              DO 50 I=1,4
   50         YNO(I)=XNO
              IN=NSECI(NE)
              JN=NSECJ(NE)
              IF (VPTP(NE,2).GT.ROTCP(IN)) YNO(1)=YES
              IF (VPTP(NE,3).GT.ROTCP(JN)) YNO(2)=YES
              IF (VPTN(NE,2).LT.ROTCN(IN)) YNO(3)=YES
              IF (VPTN(NE,3).LT.ROTCN(JN)) YNO(4)=YES
C
      write (1,60) NE,NODI(NE),NODJ(NE),VP(NE,2),VP(NE,3),VPTP(NE,2),
     1      yno(1),VPTP(NE,3),YNO(2),VPTN(NE,2),YNO(3),VPTN(NE,3),YNO(4)
   60         FORMAT(I5,I8,I7,1X,2F13.5,3X,4(F13.5,A3))
   70 CONTINUE
C
C     SUPPORT REACTIONS
C
      write (1,80)
   80 FORMAT(/////34H SUPPORT REACTIONS AT END OF CYCLE//
     15H NODE,6X,7HA-FORCE,6X,7HB-FORCE,7X,6HMOMENT,
     212X,7HX-FORCE,6X,7HY-FORCE/)
C
      DO 90 N=1,NSJTS
      NJ=NSPJ(N)
      CC=COSS(N)
      SS=SINS(N)
      FFX=SPFRA(N)*CC-SPFRB(N)*SS
      FFY=SPFRA(N)*SS+SPFRB(N)*CC
   90 write (1,100) NJ,SPFRA(N),SPFRB(N),SPFRR(N),FFX,FFY
  100 FORMAT(I4,1X,3F13.3,6X,2F13.3)
C
      RETURN
      END

C------------------------------------------------------------------------------
C
C
C
C------------------------------------------------------------------------------

      SUBROUTINE OPTSOL (A,B,NA,NEQ,LIMS,LIMF,KKK)

      Implicit None
C
      REAL A(1),B(1)
      INTEGER NA(1)
      INTEGER NEQQ,J1,IL,NAJP,J,NAJ,IF,IF1,JK,JIA,NAIP,I1,KL,I
      INTEGER NAI,IK,II,KF,JKA,IKA,K
      Integer NEQ, KK,LIMS, LIMF, MAX0, KKK
      REAL AA
      REAL CC
      INTEGER N,N1
      REAL BB
C
      NEQQ=NEQ-1
      GO TO (10,70), KKK
C     TRIANGULARIZE MATRIX A
   10 J1=3
      IL=1
      NAJP=NA(1)
      DO 60 J=2,NEQ
          NAJ=NA(J)
          IF=J1-NAJ+NAJP
          IF1=IF+1
          JK=NAJ-J
          IF (IF1.GT.IL) GO TO 40
          JIA=2+NAJP
          NAIP=NA(IF)
          I1=IF1+1
          KL=IF
          DO 30 I=IF1,IL
              NAI=NA(I)
              IK=NAI-I
              II=I1-NAI+NAIP
              KF=MAX0(II,IF)
              JKA=JK+KF
              IKA=IK+KF
              AA=A(JIA)
              DO 20 K=KF,KL
                  AA=AA-A(JKA)*A(IKA)
                  JKA=JKA+1
   20         IKA=IKA+1
              A(JIA)=AA
              JIA=JIA+1
              I1=I1+1
              KL=KL+1
   30     NAIP=NAI
C
   40     JKA=JK+IF
          AA=A(NAJ)
          DO 50 K=IF,IL
              NAI=NA(K)
              CC=A(JKA)/A(NAI)
              AA=AA-A(JKA)*CC
              A(JKA)=CC
   50     JKA=JKA+1
          A(NAJ)=AA
          IL=IL+1
          J1=J1+1
   60 NAJP=NAJ
C
      GO TO 150
C     REDUCE VECTOR B AND BACK SUBSTITUTE
   70 DO 80 N=LIMS,NEQQ
          IF (B(N).NE.0.) GO TO 90
   80 CONTINUE
      N=NEQQ
   90 N1=N+1
      I1=N1+1
      KL=N
      NAIP=NA(N)
      DO 110 I=N1,NEQ
          NAI=NA(I)
          II=I1-NAI+NAIP
          KF=MAX0(II,N)
          IK=NAI-I
          IKA=IK+KF
          BB=B(I)
          DO 100 K=KF,KL
              BB=BB-A(IKA)*B(K)
  100     IKA=IKA+1
          B(I)=BB
          I1=I1+1
          KL=KL+1
  110 NAIP=NAI
      DO 120 I=N,NEQ
          NAI=NA(I)
  120 B(I)=B(I)/A(NAI)
C
      J=NEQ
      J1=J+1
      KL=NEQQ
      NAJ=NA(NEQ)
      DO 140 I=LIMF,NEQQ
          NAJP=NA(J-1)
          II=J1-NAJ+NAJP
          JK=NAJ-J
          KF=MAX0(II,LIMF)
          JKA=JK+KF
          BB=B(J)
          DO 130 K=KF,KL
              B(K)=B(K)-A(JKA)*BB
  130     JKA=JKA+1
          J1=J1-1
          KL=KL-1
          J=J-1
  140 NAJ=NAJP
C
  150 RETURN
      END
