Ask Question

Name:
Title:
Your Question:

Answer Question

Name:
Your Answer:
User Submitted Source Code!


Description:
  1
Language: FORTRAN
Code:
      PROGRAM PCSAP4
C***********************************************************************
C
C                              PC-SAP4
C                     A STRUCTURAL ANALYSIS PROGRAM
C           FOR STATIC AND DYNAMIC RESPONSE OF LINEAR SYSTEMS
C
C                DEDICATED TO THE MEMORY OF FRED PETERSON
C
C                                BY
C                          BRUCE F. MAISON
C                            AUGUST 1994
C
C
C                     BASED ON SAP4 DEVELOPED BY
C
C                K.J. BATHE , E.L. WILSON , F.E. PETERSON
C                  UNIVERSITY OF CALIFORNIA , BERKELEY
C
C       IBM MAIN FRAME VERSION BY UNIVERSITY OF SOUTHERN CALIFORNIA
C                            AUGUST, 1973
C                         REVISED JULY, 1974
C
C***********************************************************************
C
C        ***** NOTICE TO USERS OF THIS SOFTWARE *****
C
C        COPYRIGHT (C) 1973
C        THE REGENTS OF THE UNIVERSITY OF CALIFORNIA (REGENTS)
C        ALL RIGHTS RESERVED
C
C  IN NO EVENT SHALL REGENTS OF THE UNIVERSITY OF CALIFORNIA BE
C  LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL,
C  OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF
C  THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IT REGENTS
C  HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
C
C  REGENTS SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT
C  NOT LIMITED TO, THE IMPLIED WARRANTIES OR MERCHANTABILITY AND
C  FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE AND ACCOMPANYING
C  DOCUMENTATION, IF ANY, PROVIDED HEREUNDER IS PROVIDED "AS IS".
C  REGENTS HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT,
C  UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
C
C     NO RESPONSIBILITY IS ASSUMED BY THE AUTHORS OR BY THE UNIVERSITY
C     OF CALIFORNIA FOR ANY ERRORS, MISTAKES, OR MISREPRESENTATIONS
C     THAT MAY OCCUR FROM THE USE OF THIS COMPUTER PROGRAM.  ALL
C     SOFTWARE PROVIDED ARE IN 'AS IS' CONDITION.  NO WARRANTIES OF
C     ANY KIND, WHETHER STATUTORY, WRITTEN, ORAL, EXPRESSED, OR
C     IMPLIED (INCLUDING WARRANTIES OF FITNESS AND MERCHANTABILITY)
C     SHALL APPLY.
C
C     PUBLIC DISTRIBUTION OF THIS PROGRAM THROUGH NISEE WAS MADE
C     POSSIBLE BY THE AUTHORS AND THE NATIONAL SCIENCE FOUNDATION
C     WITH THE STIPULATION THAT THE PROGRAM NEITHER BE SOLD IN WHOLE
C     OR IN PART FOR DIRECT PROFIT NOR ROYALTIES OR DEVELOPMENT
C     CHARGES MADE FOR ITS USE.  BY ACCEPTANCE OF DELIVERY OF THIS
C     PROGRAM PACKAGE, THE PURCHASER UNDERSTANDS THE RESTRICTIONS
C     ON THE USE AND DISTRIBUTION OF THE PROGRAM.  THE FEE PAID TO
C     NISEE REPRESENTS A CHARGE FOR DUPLICATION, MAILING, AND
C     DOCUMENTATION.  THE LEGAL OWNERSHIP OF THE PROGRAM REMAINS WITH
C     THE DEVELOPERS.
C
C     IF THIS PROGRAM HAS GIVEN YOU NEW ANALYSIS CAPABILITY RESULTING
C     IN INDIRECT PROFITS, YOU MAY WISH TO SUPPORT FURTHER WORK IN
C     THIS AREA BY GIVING AN UNRESTRICTED GRANT TO THE AUTHORS
C     UNIVERSITY.
C
C***********************************************************************
C
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 T,TT
      COMMON /JUNK / HED(12),JUK(460)
      COMMON /ELPAR/ NPAR(14),NUMNP,MBAND,NELTYP,N1,N2,N3,N4,N5,MTOT,NEQ
      COMMON /EM/  QQQ(2847)
      COMMON /DYN/   IDU5(11),NDYN
      COMMON /TAPES/ NQQ(6)
      COMMON /EXTRA/ MODEX,NT8,N10SV,NT10,KEQB,NUMEL,T(10)
      COMMON /SOL/   NBLOCK,NEQB,LL,NF,IDUM,NEIG,NAD,NVV,ANORM,NFO
      CHARACTER*20 NAM1,NAM2
C
C     PROGRAM CAPACITY CONTROLLED BY THE FOLLOWING TWO STATEMENTS ...
C
      COMMON A(100000001)
      MTOT=    100000000
C
      WRITE(*,6001)
      WRITE(*,6002)
      WRITE(*,6003)
      WRITE(*,6004)
      WRITE(*,6005)
      WRITE(*,6006)
      WRITE(*,6007)
      WRITE(*,6008)
      WRITE(*,6009)
C
      OPEN(1,FORM='UNFORMATTED')
      OPEN(2,FILE='NT2',FORM='UNFORMATTED')
      OPEN(3,FILE='NT3',FORM='UNFORMATTED')
      OPEN(4,FILE='NT4',FORM='UNFORMATTED')
      OPEN(7,FORM='UNFORMATTED')
      NT8 = 8
      OPEN(8,FILE='NT8',FORM='UNFORMATTED')
      REWIND NT8
      OPEN(9,FILE='NT9',FORM='UNFORMATTED')
      NT10= 10
      OPEN(10,FILE='NT10',FORM='UNFORMATTED')
      REWIND NT10
      N1=1
C
      WRITE(*,7700)
7700  FORMAT( /,1X,'* ENTER INPUT FILE NAME *',/ )
      READ(*,9900) NAM1
9900  FORMAT( A20 )
      WRITE(*,8800)
8800  FORMAT( 1X,'* ENTER OUTPUT FILE NAME *',/ )
      READ(*,9900) NAM2
      OPEN( 5, FILE = NAM1, FORM = 'FORMATTED' )
      OPEN( 6, FILE = NAM2, FORM = 'FORMATTED', STATUS = 'REPLACE' )
C
      WRITE(6,7001) NAM1,NAM2
      WRITE(6,7002)
      WRITE(6,7003)
      WRITE(6,7004)
      WRITE(6,7005)
      WRITE(6,7006)
      WRITE(6,7007)
      WRITE(6,7008)
      WRITE(6,7009)
      WRITE(6,7010)
      WRITE(6,7011)
      WRITE(6,7012)
C
C     P R O G R A M   C O N T R O L   D A T A
C
    5 CALL TTIME(T(1))
      READ (5,100,END=990) HED,NUMNP,NELTYP,LL,NF,NDYN,MODEX,NAD,
     1                     KEQB,N10SV
      IF(MODEX.GT.0) MODEX = 1
C     IF (NUMNP.EQ.0) STOP
      IF (NUMNP.EQ.0) GO TO 990
      WRITE (6,200) HED,NUMNP,NELTYP,LL,NF,NDYN,MODEX,NAD,KEQB,N10SV,
     1              MTOT
      IF(KEQB.LT.2) KEQB = 99999
      IF (NDYN.NE.0) LL=1
      IF(LL.GE.1) GO TO 10
      WRITE (6,300)
      STOP
C***  DATA PORTHOLE SAVE
   10 IF(MODEX.EQ.1)
     *WRITE (NT8)     HED,NUMNP,NELTYP,LL,NF,NDYN
C
      KDYN = IABS(NDYN) +1
      IF(KDYN.LE.5) GO TO 14
      WRITE (6,310) NDYN
      STOP
C
C     RE-START MODE ACTIVATED IF  NDYN.EQ.-2  OR  NDYN.EQ.-3
C
   14 IF(NDYN.LT.0) GO TO 20
C
C     I N P U T   J O I N T   D A T A
C
      N2=N1+6*NUMNP
      N3=N2+NUMNP
      N4=N3+NUMNP
      N5=N4+NUMNP
      N6=N5+NUMNP
      IF(N6.GT.MTOT) CALL ERROR(N6-MTOT)
C
      CALL INPUTJ(A(N1),A(N2),A(N3),A(N4),A(N5),NUMNP,NEQ)
C
C     F O R M   E L E M E N T   S T I F F N E S S E S
C
      CALL TTIME(T(2))
C
      MBAND=0
      NUMEL=0
C     OPEN(1,FORM='UNFORMATTED')
C     OPEN(1,FORM='UNFORMATTED',STATUS='SCRATCH')
      REWIND 1
C     OPEN(2,FORM='UNFORMATTED')
C     OPEN(2,FORM='UNFORMATTED',STATUS='SCRATCH')
      REWIND 2
C
      DO 900 M=1,NELTYP
      READ  (5,1001) NPAR
C***  DATA PORTHOLE SAVE
      IF(MODEX.EQ.1) WRITE (NT8) NPAR
      WRITE (1) NPAR
      NUMEL=NUMEL+NPAR(2)
      MTYPE=NPAR(1)
C
      CALL ELTYPE(MTYPE)
C
  900 CONTINUE
C
C     D E T E R M I N E   B L O C K S I Z E
C
C     ADDSTF
C
      NEQB=(MTOT - 4*LL)/(MBAND + LL + 1)/2
C
C     OVER-RIDE THE SYSTEM MATRIX BLOCKSIZE WITH THE INPUT (NON-ZERO)
C     VALUE, KEQB.
C     THIS OVER-RIDE ENTRY IS TO ALLOW PROGRAM CHECKING OF MULTI-
C     BLOCK ALGORITHMS WITH WHAT WOULD NORMALLY BE ONE BLOCK DATA.
C
      IF(KEQB.LT.NEQB) NEQB = KEQB
C
      GO TO  (690,700,700,700,730), KDYN
C
C     STATIC SOLUTION
C
  690 CONTINUE
      NEQB1=(MTOT - MBAND)/(2*(MBAND+LL) + 1)
      NEQB2=(MTOT - MBAND - LL*(MBAND-2))/(3*LL + MBAND + 1)
      IF (NEQB1.LT.NEQB) NEQB=NEQB1
      IF (NEQB2.LT.NEQB) NEQB=NEQB2
      NBLOCK = (NEQ-1)/NEQB +1
      IF(NEQB.GT.NEQ) NEQB=NEQ
      GO TO 790
C
C     EIGENSOLUTION
C
C        1.  DETERMINANT SEARCH ALGORITHM
C
  700 IF (NEQB.LT.NEQ) GO TO 710
      NIM=3
      NC=NF + NIM
      NVM=6
      NCA=NEQ*MAX0(MBAND,NC)
      NTOT=NCA + 4*NEQ + 2*NVM*NEQ + 5*NC
      NEIG=0
      IF(NTOT.LE.MTOT) GO TO 720
C
C        2.  SUBSPACE ITERATION ALGORITHM
C
  710 NV=MIN0(2*NF,NF+8)
      IF (NAD.NE.0) NV=NAD
      NEQB1=(MTOT - MBAND)/(2*MBAND + 1)
      NEQB2=(MTOT - MBAND - 2*NV - NV*(MBAND-2))/(3*NV + MBAND + 1)
      NEQB3=(MTOT - 3*NV*NV - 3*NV)/(2*NV + 1)
      NEQB4=(MTOT - 6*NV)/(1 + MBAND)
      IF (NEQB1.LT.NEQB) NEQB=NEQB1
      IF (NEQB2.LT.NEQB) NEQB=NEQB2
      IF (NEQB3.LT.NEQB) NEQB=NEQB3
      IF (NEQB4.LT.NEQB) NEQB=NEQB4
      NEIG=1
C
  720 CONTINUE
      NBLOCK = (NEQ-1)/NEQB +1
      IF (NEQB.GE.NEQ) NEQB=NEQ
C
C     HISTORY OR SPECTRUM ANALYSIS
C
      KREM = 1000
      NTOT = NBLOCK*NEQB*NF + KREM
      IF(MTOT.LT.NTOT)
     *WRITE (6,320)
      GO TO 790
C
C     STEP-BY-STEP DIRECT INTEGRATION
C
  730 CONTINUE
C     DISPLACEMENT COMPONENTS FOR DIRECT OUTPUT (*NSD*)
      NN2 = NEQ
C     DISPLACEMENT COMPONENTS REQUIRED FOR RECOVERY OF ALL OF THE
C     REQUESTED ELEMENT STRESS COMPONENTS (*NSS*)
      NN3 = NEQ
C
C        1.  DECOMPOSITION
C
      NEQB1 = (MTOT-NN2-NN3-NEQ-MBAND)/(2*MBAND+1)
C
C        2.  TIME INTEGRATION PHASE
C
      NEQB2 = (MTOT-MBAND-2*(NN2+NN3)-5*NEQ)/(MBAND+1)
C
      IF(NEQB1.LT.NEQB) NEQB = NEQB1
      IF(NEQB2.LT.NEQB) NEQB = NEQB2
      IF(NEQB.GT.NEQ)   NEQB = NEQ
      NBLOCK = (NEQ-1)/NEQB +1
C
C        3.  INPUT PHASE
C
C     NUMBER OF TIME FUNCTIONS (*NFN*)
      NN2 = 10
C     MAXIMUM NUMBER OF FUNCTION DEFINITION POINTS (*MXLP*)
      NN3 = 40
C
      NN4 = 6*NUMNP + 2*NN2*NEQ
      IF(NN4.GT.MTOT)
     *WRITE (6,320)
      NN4 = NEQ*2*(NN2+1) + NN2*(1+2*NN3)
      IF(NN4.GT.MTOT)
     *WRITE (6,320)
C
  790 CONTINUE
C
C     I N P U T   N O D A L   L O A D S
C
      N3=N2+NEQB*LL
      N4=N3+6*LL
      WRITE (6,201) NEQ,MBAND,NEQB,NBLOCK
C
      CALL TTIME(T(3))
C
      CALL INL(A(N1),A(N2),A(N3),A(N4),NUMNP,NEQB,LL)
C
      CALL TTIME(T(4))
C
C     F O R M   T O T A L   S T I F F N E S S
C
      NE2B=2*NEQB
      N2=N1+NEQB*MBAND
      N3=N2+NEQB*LL
      N4=N3+4*LL
      NN2=N1+NE2B*MBAND
      NN3=NN2+NE2B*LL
      NN4=NN3+4*LL
C
      CALL ADDSTF (A(N1),A(NN2),A(NN3),A(NN4),NUMEL,NBLOCK,NE2B,LL,MBAND
     1,ANORM,NVV)
C
      CALL TTIME(T(5))
C
C     S O L U T I O N   P H A S E
C
   20 GO TO (30,40,50,60,70), KDYN
C
C     STATIC SOLUTION
C
   30 IF(MODEX.EQ.0) GO TO 32
      DO 31 I=6,10
   31 T(I) = T(5)
      GO TO 90
C
   32 CALL SOLEQ
      CALL TTIME(T(6))
      DO 33 I=7,10
   33 T(I) = T(6)
      GO TO 90
C
C     EIGENVALUE EXTRACTION
C
   40 T(6) = T(5)
      CALL SOLEIG
      CALL TTIME(T(7))
      T(8) = T(7)
      T(9) = T(7)
      T(10)= T(7)
      GO TO 90
C
C     FORCED DYNAMIC RESPONSE ANALYSIS
C
   50 T(6) = T(5)
      IF(NDYN.LT.0) GO TO 52
      CALL SOLEIG
      CALL TTIME (T(7))
      GO TO 54
   52 DO 53 I=1,6
   53 T(I+1)=T(I)
      REWIND 2
      READ (2) NEQ,NBLOCK,NEQB,MBAND,N1,NF,(QQQ(I),I=1,NF)
      REWIND 7
      IMAX=NEQB*NF
      READ (7) (A(I),I=1,NF)
      DO 56 L=1,NBLOCK
   56 READ (7) (A(I),I=1,IMAX)
   54 CALL HISTRY
      CALL TTIME (T(8))
      T(9) = T(8)
      T(10)= T(8)
      GO TO 90
C
C     RESPONSE SPECTRUM ANALYSIS
C
   60 T(6) = T(5)
      IF(NDYN.LT.0) GO TO 62
      CALL SOLEIG
      CALL TTIME (T(7))
      T(8) = T(7)
      GO TO 64
   62 DO 63 I=1,7
   63 T(I+1)=T(I)
      REWIND 2
      READ (2) NEQ,NBLOCK,NEQB,MBAND,N1,NF
      REWIND 7
      IMAX=NEQB*NF
      READ (7) (A(I),I=1,NF)
      DO 66 L=1,NBLOCK
   66 READ (7) (A(I),I=1,IMAX)
   64 CALL RESPEC
      CALL TTIME (T(9))
      T(10)= T(9)
      GO TO 90
C
C     STEP-BY-STEP (DIRECT INTEGRATION) ANALYSIS
C
   70 DO 71 I=6,9
   71 T(I) = T(5)
      CALL STEP
      CALL TTIME(T(10))
C
C     COMPUTE AND PRINT OVERALL TIME LOG
C
C##90 TT = 0.0
   90 TT = 0.0D0
      DO 95 I=1,9
      T(I) = T(I+1)-T(I)
      TT = TT + T(I)
   95 CONTINUE
C
      WRITE (6,203) (T(K),K=1,9),TT
C
      GO TO 5
C 990 STOP
  990 CONTINUE
      CLOSE(1,STATUS='DELETE')
      CLOSE(2,STATUS='DELETE')
      CLOSE(3,STATUS='DELETE')
      CLOSE(4,STATUS='DELETE')
      CLOSE(7,STATUS='DELETE')
      CLOSE(8,STATUS='DELETE')
      CLOSE(9,STATUS='DELETE')
      CLOSE(10,STATUS='DELETE')
C
  100 FORMAT (12A6/9I5)
  200 FORMAT(1H ,12A6///
     1  38H C O N T R O L   I N F O R M A T I O N, // 4X,
     2  27H NUMBER OF NODAL POINTS   =, I5 / 4X,
     3  27H NUMBER OF ELEMENT TYPES  =, I5 / 4X,
     4  27H NUMBER OF LOAD CASES     =, I5 / 4X,
     5  27H NUMBER OF FREQUENCIES    =, I5 / 4X,
     6  27H ANALYSIS CODE (NDYN)     =, I5 / 4X,
     7  16H   EQ.0,  STATIC,               / 4X,
     8  26H   EQ.1,  MODAL EXTRACTION,     / 4X,
     9  25H   EQ.2,  FORCED RESPONSE,      / 4X,
     A  27H   EQ.3,  RESPONSE SPECTRUM,    / 4X,
     *  28H   EQ.4,  DIRECT INTEGRATION,   / 4X,
     B  27H SOLUTION MODE (MODEX)    =, I5 / 4X,
     C  19H   EQ.0,  EXECUTION,            / 4X,
     D  20H   EQ.1,  DATA CHECK,           / 4X,
     E  19H NUMBER OF SUBSPACE,            / 4X,
     F  27H ITERATION VECTORS (NAD)  =, I5 / 4X,
     G  27H EQUATIONS PER BLOCK      =, I5 / 4X,
     H  27H TAPE10 SAVE FLAG (N10SV) =, I5 / 4X,
     I  27H STORAGE CAPACITY (MTOT)  =, I10 / )
  201 FORMAT (38H E Q U A T I O N   P A R A M E T E R S, //
     *       34H  TOTAL NUMBER OF EQUATIONS      =,I5,
     1      /34H  BANDWIDTH                      =,I5,
     2      /34H  NUMBER OF EQUATIONS IN A BLOCK =,I5,
     3      /34H  NUMBER OF BLOCKS               =,I5)
  203 FORMAT (// 1H ,31HO V E R A L L   T I M E   L O G, //
     1 5X,30HNODAL POINT INPUT            =, F8.3 /
     2 5X,30HELEMENT STIFFNESS FORMATION  =, F8.3 /
     3 5X,30HNODAL LOAD INPUT             =, F8.3 /
     4 5X,30HTOTAL STIFFNESS FORMATION    =, F8.3 /
     5 5X,30HSTATIC ANALYSIS              =, F8.3 /
     6 5X,30HEIGENVALUE EXTRACTION        =, F8.3 /
     7 5X,30HFORCED RESPONSE ANALYSIS     =, F8.3 /
     8 5X,30HRESPONSE SPECTRUM ANALYSIS   =, F8.3 /
     * 5X,30HSTEP-BY-STEP INTEGRATION     =, F8.3 //
     9 5X,30HTOTAL SOLUTION TIME          =, F8.3 /)
C
  300 FORMAT (// 48H ** ERROR.  (AT LEAST ONE LOAD CASE IS REQUIRED)  )
  310 FORMAT (// 33H ** ERROR.  ANALYSIS CODE (NDYN =,I3,9H) IS BAD.  )
  320 FORMAT (// 47H ** WARNING.  ESTIMATE OF STORAGE FOR A DYNAMIC,
     1           32H ANALYSIS EXCEEDS AVAILABLE CORE, // 1X)
C
 1001 FORMAT (14I5)
6001  FORMAT( //,35X,'Welcome To',/// )
6002  FORMAT( 9X,'                       SSSSS    AAAA    PPPPP       
     *   ' )
6003  FORMAT( 8X,'                       S       A    A   P    P      
     *   ' )
6004  FORMAT( 7X,'                       SSSSS   AAAAAA   PPPPP        
     *     ' )
6005  FORMAT( 6X,'                           S   A    A   P
     *   ' )
6006  FORMAT( 5X,'                       SSSSS   A    A   P
     *   ',// )
6007  FORMAT( 39X,'By',/ )
6008  FORMAT( 27X,'Dan Marinescu & Tim Benson',/ )
6009  FORMAT( 36X,'May 2010',/ )
7001  FORMAT( 1X,'INPUT FILE: ',A20,20X,'OUTPUT FILE: ',A20,5(/) )
7002  FORMAT( 9X,'                       SSSSS    AAAA    PPPPP       
     *   ' )
7003  FORMAT( 8X,'                       S       A    A   P    P      
     *   ' )
7004  FORMAT( 7X,'                       SSSSS   AAAAAA   PPPPP        
     *     ' )
7005  FORMAT( 6X,'                           S   A    A   P
     *   ' )
7006  FORMAT( 5X,'                       SSSSS   A    A   P
     *   ',/ )
7007  FORMAT( 38X,'By',/ )
7008  FORMAT( 27X,'Dan Marinescu & Tim Benson' )
7009  FORMAT( 36X,'May 2010',2(/) )
7010  FORMAT( 26X,'Based On SAP-4 Developed By' )
7011  FORMAT( 21X,'K.J. Bathe, E.L. Wilson, F.E. Peterson' )
7012  FORMAT( 23X,'University of California, Berkeley',/// )
      END
      SUBROUTINE ADDMAS (TMASS,BLKMAS,NEQ,NEQB,NBLOCK)
      IMPLICIT REAL*8(A-H,O-Z)
C
C     CALLED BY:  STEP
C
C     THIS ROUTINE READS THE SYSTEM MASS MATRIX IN BLOCKED FORM
C     FROM *TAPE9* AND ASSEMBLES THE BLOCKS INTO A SINGLE VECTOR
C     *NEQ* WORDS IN LENGTH -- I.E., SYSTEM MASS MATRIX (DIAGONAL)
C     IS STORED IN CORE.  SYSTEM MASS MATRIX *TMASS* IS SAVED ON
C     *TAPE3*.
C
      DIMENSION      TMASS(NEQ),BLKMAS(NEQB)
C
      NT3 = 3
      REWIND NT3
      NT9 = 9
      REWIND NT9
C
      KSHIFT = 0
C
C     LOOP ON THE TOTAL NUMBER OF SYSTEM EQUATION BLOCKS
C
      DO 200 K=1,NBLOCK
      READ (NT9) BLKMAS
      K1 = KSHIFT
      DO 100 L=1,NEQB
      K1 = K1+1
      IF(K1.GT.NEQ) GO TO 250
      TMASS(K1) = BLKMAS(L)
  100 CONTINUE
      KSHIFT = KSHIFT+NEQB
  200 CONTINUE
C
  250 WRITE (NT3) TMASS
C
      RETURN
      END
      SUBROUTINE ADDSTF (A,B,STR,TMASS,NUMEL,NBLOCK,NE2B,LL,MBAND,ANORM,
     1NVV)
C
      IMPLICIT REAL*8(A-H,O-Z)
C
C     CALLED BY:  MAIN
C
C     FORMS GLOBAL EQUILIBRIUM EQUATIONS IN BLOCKS
C
      DIMENSION A(NE2B,MBAND),B(NE2B,LL),STR(4,LL),TMASS(NE2B)
C
      COMMON /DYN/   NT,NOT,ALFA,DT,BETA,NFN,NGM,NAT,NDYN
      COMMON /EM/ LRD,ND,LM(63),IPAD,SS(2331),IDUM(966)
      COMMON /EXTRA/ MODEX,NT8,IFILL(24)
C
      NEQB=NE2B/2
      K=NEQB+1
      X=NBLOCK
      MB=DSQRT(X)
      MB=MB/2+1
      NEBB=MB*NE2B
      MM=1
      NDEG=0
      NVV=0
C##   ANORM=0.
      ANORM=0.D0
      NSHIFT=0
      REWIND 3
C     OPEN(4,FORM='UNFORMATTED')
C     OPEN(4,FORM='UNFORMATTED',STATUS='SCRATCH')
      REWIND 4
C     OPEN(9,FORM='UNFORMATTED')
C     OPEN(9,FORM='UNFORMATTED',STATUS='SCRATCH')
      REWIND 9
C
C     READ ELEMENT LOAD MULTIPLIERS
C
      WRITE (6,2000)
      DO 50 L=1,LL
      READ  (5,1002)   (STR(I,L),I=1,4)
   50 WRITE (6,2002) L,(STR(I,L),I=1,4)
      IF(MODEX.EQ.0) WRITE (8) STR
C
C     FOR A STEP-BY-STEP ANALYSIS  (NDYN.EQ.4)  READ THE SOLUTION
C     CONTROL CARD.  THE TIME STEP (DT) AND THE DAMPING COEFFICIENTS
C     (ALFA/BETA) ARE REQUIRED FOR THE ASSEMBLY OF THE EFFECTIVE
C     SYSTEM STIFFNESS MATRIX IN THIS ROUTINE.
C
      IF(NDYN.NE.4) GO TO 65
C
      READ  (5,1004) NFN,NGM,NAT,NT,NOT,DT,ALFA,BETA
      WRITE (6,2004) NFN,NGM,NAT,NT,NOT,DT,ALFA,BETA
      IF(NAT.EQ.0) NAT = 1
      IF(NOT.EQ.0) NOT = 1
C##   IF(DT.GT.1.0E-12) GO TO 55
      IF(DT.GT.1.0D-12) GO TO 55
      WRITE (6,3000)
      STOP
C
C     COMPUTE INTEGRATION COEFFICIENTS FOR ASSEMBLY OF EFFECTIVE
C     SYSTEM STIFFNESS (STEP-BY-STEP ANALYSIS ONLY)
C
C##55 TETA = 1.4
   55 TETA = 1.4D0
      DT1  = TETA*DT
      DT2  = DT1**2
C##   A0   = (6.+3.*ALFA*DT1)/(DT2+3.*BETA*DT1)
      A0   = (6.D0+3.D0*ALFA*DT1)/(DT2+3.D0*BETA*DT1)
C
   65 IF(MODEX.EQ.1) RETURN
C
C     FORM EQUATIONS IN BLOCKS     ( 2 BLOCKS AT A TIME)
C
      DO 1000 M=1,NBLOCK ,2
      DO 100 I=1,NE2B
      DO 100 J=1,MBAND
C##00 A(I,J)=0.
  100 A(I,J)=0.D0
      READ (3) ((B(I,L),I=1,NEQB),L=1,LL),(TMASS(I),I=1,NEQB)
      IF (M.EQ.NBLOCK) GO TO 200
      READ (3) ((B(I,L),I=K,NE2B),L=1,LL),(TMASS(I),I=K,NE2B)
  200 CONTINUE
C
C     OPEN(7,FORM='UNFORMATTED')
C     OPEN(7,FORM='UNFORMATTED',STATUS='SCRATCH')
      REWIND 7
      REWIND 2
      NA=7
      NUME=NUM7
      IF (MM.NE.1) GO TO 75
      NA=2
      NUME=NUMEL
      NUM7 =0
C
   75 DO 700 N=1,NUME
      READ (NA) LRD,ND,(LM(I),I=1,ND),(SS(I),I=1,LRD)
      MSHFT = ND * (ND+1)/2 +4 *ND
      DO 600 I=1,ND
      LMN=1-LM(I)
      II=LM(I)-NSHIFT
      IF (II.LE.0.OR.II.GT.NE2B) GO TO 600
      IMS=I+MSHFT
      TMASS(II)=TMASS(II)+ SS(IMS)
      DO 300 L=1,LL
      DO 300 J=1,4
      KK = ND *(ND+1)/2 + ND*(J-1)
  300 B(II,L)=B(II,L)+SS(I+KK)*STR(J,L)
      DO 500 J=1,ND
      JJ=LM(J)+LMN
      IF(JJ) 500,500,390
  390 IF(J-I) 396,394,394
  394 KK = ND*I-(I-1)*I/2 +J-ND
      GO TO 400
  396 KK =ND*J -(J-1)*J/2+I-ND
  400 A(II,JJ)=A(II,JJ)+SS(  KK)
  500 CONTINUE
  600 CONTINUE
C
C     DETERMINE IF STIFFNESS IS TO BE PLACED ON TAPE 7
C
      IF (MM.GT.1) GO TO 700
      DO 650 I=1,ND
      II=LM(I) -NSHIFT
      IF(II.GT.NE2B.AND.II.LE.NEBB) GO TO 660
  650 CONTINUE
      GO TO 700
  660 WRITE (7) LRD,ND,(LM(I),I=1,ND),(SS(I),I=1,LRD)
      NUM7=NUM7+1
C
  700 CONTINUE
      DO 710 L=1,NEQB
      ANORM=ANORM + A(L,1)
C##   IF (A(L,1).NE.0.) NDEG=NDEG + 1
      IF (A(L,1).NE.0.D0) NDEG=NDEG + 1
C##   IF (A(L,1).EQ.0.) A(L,1)=1.E+20
      IF (A(L,1).EQ.0.D0) A(L,1)=1.D+20
C##   IF (TMASS(L).NE.0.) NVV=NVV + 1
      IF (TMASS(L).NE.0.D0) NVV=NVV + 1
  710 CONTINUE
C
C     FOR STEP-BY-STEP ANALYSIS ADD THE MASS CONTRIBUTIONS TO
C     THE EQUATION DIAGONAL COEFFICIENTS
C
      IF(NDYN.NE.4) GO TO 716
      DO 714 I=1,NEQB
  714 A(I,1) = A(I,1) + A0* TMASS(I)
      WRITE (4) ((A(I,J),I=1,NEQB),J=1,MBAND)
      GO TO 718
  716 WRITE (4) ((A(I,J),I=1,NEQB),J=1,MBAND),((B(I,L),I=1,NEQB),L=1,LL)
  718 WRITE (9)  (TMASS(I),I=1,NEQB)
C
      IF(M.EQ.NBLOCK) GO TO 1000
      DO 720 L=K,NE2B
      ANORM=ANORM + A(L,1)
C##   IF (A(L,1).NE.0.) NDEG=NDEG + 1
      IF (A(L,1).NE.0.D0) NDEG=NDEG + 1
C##   IF (A(L,1).EQ.0.) A(L,1)=1.E+20
      IF (A(L,1).EQ.0.D0) A(L,1)=1.D+20
C##   IF (TMASS(L).NE.0.) NVV=NVV + 1
      IF (TMASS(L).NE.0.D0) NVV=NVV + 1
  720 CONTINUE
C
      IF(NDYN.NE.4) GO TO 726
      DO 724 I=K,NE2B
  724 A(I,1) = A(I,1) + A0* TMASS(I)
      WRITE (4) ((A(I,J),I=K,NE2B),J=1,MBAND)
      GO TO 728
  726 WRITE (4) ((A(I,J),I=K,NE2B),J=1,MBAND),((B(I,L),I=K,NE2B),L=1,LL)
  728 WRITE (9)  (TMASS(I),I=K,NE2B)
C
      IF (MM.EQ.MB) MM=0
      MM=MM+1
 1000 NSHIFT=NSHIFT+NE2B
      IF (NDEG.GT.0) GO TO 730
      WRITE (6,1010)
      STOP
C##30 ANORM=(ANORM/NDEG)*1.E-8
  730 ANORM=(ANORM/NDEG)*1.D-8
C
      RETURN
 1002 FORMAT (4F10.0)
 1004 FORMAT (5I5,3F10.0)
 1010 FORMAT (51H STRUCTURE WITH NO DEGREES OF FREEDOM  CHECK DATA   )
 2000 FORMAT (/// 10H STRUCTURE,13X,7HELEMENT,4X,4HLOAD,4X,
     1 11HMULTIPLIERS,/ 10H LOAD CASE,12X,1HA,9X,1HB,9X,1HC,9X,1HD,/ 1X)
 2002 FORMAT (I6,7X,4F10.3)
 2004 FORMAT (45H S T E P - B Y - S T E P   S O L U T I O N   ,
     1        37HC O N T R O L   I N F O R M A T I O N, ///
     2 5X, 35HNUMBER OF TIME VARYING FUNCTIONS  =, I5    //
     3 5X, 35HGROUND MOTION INDICATOR           =, I5    /
     4 8X, 10HEQ.0, NONE, /
     5 8X, 29HGT.0, READ ACCELERATION INPUT, //
     6 5X, 35HNUMBER OF ARRIVAL TIMES           =, I5    /
     7 8X, 26HEQ.0, ALL FUNCTIONS ARRIVE, /
     8 8X, 18H      AT TIME ZERO, //
     9 5X, 35HNUMBER OF SOLUTION TIME STEPS     =, I5    //
     A 5X, 35HOUTPUT (PRINT) INTERVAL           =, I5    //
     B 5X, 35HSOLUTION TIME INCREMENT           =, E14.4 //
     C 5X, 30HMASS-     PROPORTIONAL DAMPING, /
     D 5X, 35HCOEFFICIENT (ALPHA)               =, E14.4 //
     E 5X, 30HSTIFFNESS-PROPORTIONAL DAMPING, /
     F 5X, 35HCOEFFICIENT (BETA)                =, E14.4 /// 1X)
 3000 FORMAT (27H *** ERROR   ZERO TIME STEP, / 1X)
      END
       SUBROUTINE BANDET (A,B,V,MAXA,NN,NWA,RA,NSCH,DET,ISCALE,KK)
      IMPLICIT REAL*8(A-H,O-Z)
C
C     CALLED BY:  SECNTD
C
       COMMON /TAPES/NSTIF,NRED,NL,NR,NT,NMASS
       DIMENSION A(NWA),B(1),V(1),MAXA(1)
C
       NR=NN-1
       IF (KK-2) 100,700,800
C
C##0   TOL=1.0E+07
 100   TOL=1.0D+07
C##    RTOL=1.0E-10
       RTOL=1.0D-10
C     SCALE=2.0D0**200
      SCALE=2.D0**166
       NTF=3
       IS=1
 120   REWIND NSTIF
       READ (NSTIF) A
       DO 140 I=1,NN
 140   A(I)=A(I)-RA*B(I)
 160   IF (NWA.EQ.NN) GO TO 230
       DO 200 N=1,NR
       IH=N+NWA-NN
 210   IF (A(IH)) 220,215,220
 215   IH=IH-NN
       GO TO 210
 220   MAXA(N)=IH
       PIV=A(N)
      IF(PIV) 221,500,221
  500 IS = IS+1
      IF(IS.LE.NTF) GO TO 502
  501 WRITE (6,1000) NTF,RA
      STOP
C##02 RA = RA*(1.0-RTOL)
  502 RA = RA*(1.0D0-RTOL)
      GO TO 120
 221   IL=N+NN
       L=N
       DO 240 I=IL,IH,NN
       L=L+1
       C=A(I)
       IF (C) 225,240,225
 225   C=C/PIV
      IF (DABS(C) .LT. TOL) GO TO 235
 226   IS=IS+1
       IF (IS.LE.NTF) GO TO 245
      GO TO 501
C##5   RA=RA*(1.0-RTOL)
 245   RA=RA*(1.0D0-RTOL)
       GO TO 120
 235   J=L-I
       DO 260 K=I,IH,NN
 260   A(K+J)=A(K+J)-C*A(K)
       A(I)=C
 240   CONTINUE
 200   CONTINUE
C##0   IF (A(NN).NE.0.0) GO TO 280
 230   IF (A(NN).NE.0.0D0) GO TO 280
      AA=DABS(A(1))
       DO 290 I=2,NR
 290  AA=AA+DABS(A(I))
C##    A(NN)=-(AA/NR)*1.0E-16
       A(NN)=-(AA/NR)*1.0D-16
C
 280   NSCH=0
       ISC=0
C##    DET=1.0
       DET=1.0D0
       DO 300 I=1,NN
      IF (DABS(DET) .LT. SCALE) GO TO 320
       DET=DET/SCALE
       ISC=ISC+1
 320   DET=DET*A(I)
C##0   IF (A(I).LT.0.) NSCH=NSCH+1
 300   IF (A(I).LT.0.D0) NSCH=NSCH+1
C
       IF (ISCALE.LT.1000) GO TO 340
       ISCALE=ISC
       GO TO 900
 340   IF (ISC-ISCALE) 350,900,370
 350   DET=DET/SCALE
       GO TO 900
 370   DET=DET*SCALE
       GO TO 900
C
 700   IL=NN
       DO 400 N=1,NR
       C=V(N)
       V(N)=C/A(N)
       IF (NWA-NN) 410,400,410
 410   IL=IL+1
       IH=MAXA(N)
       K=N
       DO 420 I=IL,IH,NN
       K=K+1
 420   V(K)=V(K)-C*A(I)
 400   CONTINUE
       V(NN)=V(NN)/A(NN)
C
 800   IF (NWA-NN) 430,900,430
 430   N=NN
       DO 440 L=2,NN
       N=N-1
       IL=N+NN
       IH=MAXA(N)
       K=N
       DO 460 I=IL,IH,NN
       K=K+1
 460   V(N)=V(N)-A(I)*V(K)
 440   CONTINUE
 900   RETURN
C
 1000 FORMAT (37H ***ERROR   SOLUTION STOP IN *BANDET*, / 12X,
     1        1H(,I3,37H) TRIANGULAR FACTORIZATIONS ATTEMPTED, / 12X,
     2        16HCURRENT SHIFT = ,E20.14 / 1X)
C
       END
      SUBROUTINE BEAM
C
      IMPLICIT REAL*8(A-H,O-Z)
C
C     CALLS:  TEAM,STRSC
C     CALLED BY:  ELTYPE
C
      COMMON /ELPAR/ NPAR(14),NUMNP,MBAND,NELTYP,N1,N2,N3,N4,N5,MTOT,NEQ
      COMMON /JUNK/ LT,LH,L,IPAD,SIG(20),N6,N7,N8,N9,N10,IFILL(435)
      COMMON /EXTRA/ MODEX,NT8,N10SV,NT10,IFILL2(22)
C---      COMMON A(1)
      COMMON A(100000001)
C
      IF(NPAR(1).EQ.0) GO TO 500
      N5A=N5+NUMNP
      N6=N5+NPAR(5) + NUMNP
      N7=N6+NPAR(5)
      N8=N7+NPAR(5)
      N9=N8+12*NPAR(4)
      N10=N9+6*NPAR(3)
      N11=N10+NPAR(5)
      IF(N11.GT.MTOT) CALL ERROR(N11-MTOT)
C
      CALL TEAM(NPAR(2),NPAR(3),NPAR(4),NPAR(5),A(N1),A(N2),A(N3),
     1          A(N4),A(N5A),A(N6),A(N7),A(N8),A(N9),A(N10),
     2          NUMNP,MBAND)
C
      RETURN
C
  500 WRITE (6,2002)
      NUME=NPAR(2)
      DO 800 MM=1,NUME
      CALL STRSC (A(N1),A(N3),NEQ,0)
      WRITE (6,2001)
      DO 800 L=LT,LH
      CALL STRSC (A(N1),A(N3),NEQ,1)
      WRITE(6,3002) MM,L,(SIG(I),I=1,12)
C***  STRESS PORTHOLE
      IF(N10SV.EQ.1)
     *WRITE (NT10) MM,L,(SIG(I),I=1,12)
  800 CONTINUE
      RETURN
 2001 FORMAT (/)
 2002 FORMAT(/29H .....BEAM FORCES AND MOMENTS//
     . 10H BEAM LOAD,5X,5HAXIAL, 2(7X,5HSHEAR),5X,7HTORSION,
     .  2(5X,7HBENDING)/ 10H  NO.  NO. ,8X, 2HR1, 10X ,2HR2 ,10X,
     .  2HR3,10X,2HM1,10X,2HM2,10X,2HM3)
 3002 FORMAT (I5,I4,1PE11.3,5E12.3/8X,6E12.3/)
      END
      SUBROUTINE BENDDC (NEL,NI,NJ,X1,X2,X3,R,KODE,A,MODEX,THETA,TOL,PI)
      IMPLICIT REAL*8(A-H,O-Z)
C
C     CALLED BY:  PIPEK
C
C     COMPUTATION OF DIRECTION COSINE ARRAY FOR THE LOCAL AXES OF A
C     CIRCULAR BEND PIPE ELEMENT
C
C     NEL          =  ELEMENT NUMBER
C     NI           =  NODE NUMBER AT END I
C     NJ           =  NODE NUMBER AT END J
C     X1           =  GLOBAL COORDINATES OF END I
C     X2           =  GLOBAL COORDINATES OF END J
C     X3           =  GLOBAL COORDINATES OF THE THIRD POINT
C     KODE         =  CODE DEFINING THE THIRD POINT
C                     (EQ.1, TANGENT INTERSECTION POINT)
C                     (EQ.2, CENTER OF CURVATURE)
C     R            =  RADIUS OF THE BEND
C     A            =  MATRIX OF DIRECTION COSINES RELATING LOCAL TO THE
C                     GLOBAL SYSTEM.  A(I,J) IS THE PROJECTION ON THE
C                     I-TH GLOBAL AXIS OF A UNIT VECTOR IN THE LOCAL
C                     J-DIRECTION.
C     MODEX        =  EXECUTION MODE
C                     (EQ.0, SOLUTION)
C                     (EQ.1, DATA CHECK)
C     THETA        =  CENTRAL ANGLE SUBTENDED BY THE ARC OF THE BEND
C     TOL          =  DIMENSIONAL TOLERANCE USED FOR ERROR TESTING
C     PI           =  3.14159...
C
      DIMENSION X1(3),X2(3),X3(3),A(3,3),B(3)
C
      GO TO (10,110),KODE
C
C     TANGENT INTERSECTION IS THE THIRD POINT
C
C        1. LOCAL X-AXIS VECTOR
C
   10 A(1,1) = X3(1)-X1(1)
      A(2,1) = X3(2)-X1(2)
      A(3,1) = X3(3)-X1(3)
      XL1T = A(1,1)**2 + A(2,1)**2 + A(3,1)**2
      XL1T =DSQRT(XL1T)
C##   IF(XL1T.GT.1.0E-8) GO TO 20
      IF(XL1T.GT.1.0D-8) GO TO 20
      NN = NI
   15 WRITE (6,3000) NEL,NN
      MODEX = 1
      RETURN
C##20 DUM = 1.0/ XL1T
   20 DUM = 1.0D0/ XL1T
      DO 25 K=1,3
   25 A(K,1) = A(K,1)* DUM
C
C        2. VECTOR FROM TANGENT POINT TO NODE J
C
      DO 30 K=1,3
   30 B(K) = X2(K)-X3(K)
      XLT2 = B(1)**2 + B(2)**2 + B(3)**2
      XLT2 =DSQRT(XLT2)
C##   IF(XLT2.GT.1.0E-8) GO TO 40
      IF(XLT2.GT.1.0D-8) GO TO 40
      NN = NJ
      GO TO 15
C
C        3. COMPARE DISTANCES BETWEEN THE NODES AND THE COMMON TANGENT
C           INTERSECTION POINT
C
   40 DIF =DABS(XL1T-XLT2)
      IF(DIF.LE.TOL) GO TO 42
      WRITE (6,3010) NEL,TOL,XL1T,XLT2
      MODEX = 1
      RETURN
C
   42 CONTINUE
C
C        4. LOCAL Z-AXIS
C
      A(1,3) = A(2,1)*B(3) - A(3,1)*B(2)
      A(2,3) = A(3,1)*B(1) - A(1,1)*B(3)
      A(3,3) = A(1,1)*B(2) - A(2,1)*B(1)
C##   DUM = 0.0
      DUM = 0.0D0
      DO 44 K=1,3
   44 DUM = DUM + A(K,3)**2
      DUM =DSQRT(DUM)
C##   IF(DUM.GT.1.0E-8) GO TO 46
      IF(DUM.GT.1.0D-8) GO TO 46
      WRITE (6,3060) NEL
      MODEX = 1
      RETURN
C##46 DUM = 1.0/DUM
   46 DUM = 1.0D0/DUM
      DO 48 K=1,3
   48 A(K,3) = A(K,3) * DUM
C
C        5. LOCAL Y-AXIS
C
      A(1,2) = A(2,3)*A(3,1) - A(3,3)*A(2,1)
      A(2,2) = A(3,3)*A(1,1) - A(1,3)*A(3,1)
      A(3,2) = A(1,3)*A(2,1) - A(2,3)*A(1,1)
C
C        6. COMPUTE THE CENTRAL ANGLE
C
      DUM = XL1T/R
      THETA = 2.0D0*DATAN(DUM)
   50 CONTINUE
C##   IF(THETA.GT.1.0E-8 .AND. THETA.LE.PI) RETURN
      IF(THETA.GT.1.0D-8 .AND. THETA.LE.PI) RETURN
C##   DUM = THETA*180.0/PI
      DUM = THETA*180.0D0/PI
      WRITE (6,3020) DUM,NEL
      MODEX = 1
      RETURN
C
C     CENTER OF CURVATURE IS THE THIRD POINT
C
C        1. LOCAL Y-AXIS VECTOR
C
  110 A(1,2) = X3(1)-X1(1)
      A(2,2) = X3(2)-X1(2)
      A(3,2) = X3(3)-X1(3)
C##   D1C = 0.0
      D1C = 0.0D0
      DO 120 K=1,3
  120 D1C = D1C + A(K,2)**2
      D1C =DSQRT(D1C)
C##   IF(D1C.GT.1.0E-8) GO TO 130
      IF(D1C.GT.1.0D-8) GO TO 130
      NN = NI
  125 WRITE (6,3030) NEL,NN
      MODEX = 1
      RETURN
C##30 DUM = 1.0/ D1C
  130 DUM = 1.0D0/ D1C
      DO 135 K=1,3
  135 A(K,2) = A(K,2)* DUM
C
C        2. COMPUTE THE VECTOR FROM NODE J TO THE C.C.
C
      B(1) = X3(1)-X2(1)
      B(2) = X3(2)-X2(2)
      B(3) = X3(3)-X2(3)
C##   D2C = 0.0
      D2C = 0.0D0
      DO 140 K=1,3
  140 D2C = D2C + B(K)**2
      D2C =DSQRT(D2C)
C##   IF(D2C.GT.1.0E-8) GO TO 150
      IF(D2C.GT.1.0D-8) GO TO 150
      NN = NJ
      GO TO 125
  150 CONTINUE
C
C        3. COMPARE COMPUTED RADII VERSUS THE INPUT VALUE
C
      DIF =DABS(R-D1C)
      IF(DIF.LT.TOL) GO TO 165
      NN = NI
      RR = D1C
  160 WRITE (6,3040) NN,NEL,R,RR
      MODEX = 1
      RETURN
  165 DIF =DABS(R-D2C)
      IF(DIF.LT.TOL) GO TO 170
      NN = NJ
      RR = D2C
      GO TO 160
C
C        4. LOCAL Z-AXIS VECTOR
C
  170 A(1,3) = A(2,2)*B(3) - A(3,2)*B(2)
      A(2,3) = A(3,2)*B(1) - A(1,2)*B(3)
      A(3,3) = A(1,2)*B(2) - A(2,2)*B(1)
C##   DUM = 0.0
      DUM = 0.0D0
      DO 172 K=1,3
  172 DUM = DUM + A(K,3)**2
      DUM =DSQRT(DUM)
C##   IF(DUM.LT.1.0E-8) GO TO 177
      IF(DUM.LT.1.0D-8) GO TO 177
C##   DUM = 1.0/DUM
      DUM = 1.0D0/DUM
      DO 173 K=1,3
  173 A(K,3) = A(K,3) * DUM
C
C        5. TEST FOR NODES I AND J COINCIDENT
C
C##   CHORD = 0.0
      CHORD = 0.0D0
      DO 175 K=1,3
  175 CHORD = CHORD + (X2(K)-X1(K))**2
      CHORD =DSQRT(CHORD)
C##   IF(CHORD.GT.1.0E-8) GO TO 180
      IF(CHORD.GT.1.0D-8) GO TO 180
  177 WRITE (6,3050) NI,NJ,NEL
      MODEX = 1
      RETURN
C
C        6. LOCAL X-AXIS VECTOR
C
  180 A(1,1) = A(2,2)*A(3,3) - A(3,2)*A(2,3)
      A(2,1) = A(3,2)*A(1,3) - A(1,2)*A(3,3)
      A(3,1) = A(1,2)*A(2,3) - A(2,2)*A(1,3)
C
C        7. COMPUTE THE CENTRAL ANGLE
C
C##   DUM = 0.5*CHORD/R
      DUM = 0.5D0*CHORD/R
C--   THETA = 2.0D0*DARSIN(DUM)
      THETA = 2.0D0*DASIN(DUM)
      GO TO 50
C
 3000 FORMAT (25H ERROR***  BEND ELEMENT (,I4,19H) HAS ZERO DISTANCE,
     1 15H BETWEEN NODE (,I4,31H) AND THE TANGENT INTERSECTION., / 1X)
 3010 FORMAT (45H ERROR***  TANGENT LENGTHS FOR BEND ELEMENT (,I4,
     1 27H) ARE NOT EQUAL TO WITHIN (,E11.4, 2H)., /
     2 11X,23HI-NODE TANGENT LENGTH =,E20.8, /
     3 11X,23HJ-NODE TANGENT LENGTH =,E20.8, / 1X)
 3020 FORMAT (30H ERROR***  THE CENTRAL ANGLE (,F8.3,10H) FOR BEND,
     1 10H ELEMENT (,I4,18H) IS OUT OF RANGE., / 11X,
     2 38HTHETA MUST BE GT.0 AND LT.180 DEGREES., / 1X)
 3030 FORMAT (25H ERROR***  BEND ELEMENT (,I4,19H) HAS ZERO DISTANCE,
     1 15H BETWEEN NODE (,I4,30H) AND THE CENTER OF CURVATURE.,/ 1X)
 3040 FORMAT (36H ERROR***  COMPUTED RADIUS TO NODE (,I4,10H) FOR BEND,
     1 10H ELEMENT (,I4,38H) IS DISCREPANT FROM THE INPUT RADUIS., /
     2 11X, 17HRADIUS INPUT    =,E20.8 /
     3 11X, 17HRADIUS COMPUTED =,E20.8 / 1X)
 3050 FORMAT (44H ERROR***  ZERO CHORD LENGTH BETWEEN NODES (,I4,
     1 7H) AND (,I4,19H) IN BEND ELEMENT (,I4,2H)., / 1X)
 3060 FORMAT (51H ERROR***  TANGENT INTERSECTION POINT FOR ELEMENT (,
     1        I4,18H) IS ON THE CHORD., / 1X)
C
      END
C
C     CALLS:  PINVER
C     CALLED BY:  PIPEK
C
C     COMPUTATION OF THE ELEMENT STIFFNESS AND LOAD MATRICES FOR A
C     CIRCULARLY CURVED PIPE BEND ELEMENT.
C
C
C     ALFAV        =  SHAPE FACTOR FOR SHEAR DISTORTION
C                     (GT.99.99, NEGLECT)
C     NOAX         =  CODE FOR NEGLECTING AXIAL DEFORMATIONS
C                     (EQ.1, NEGLECT)
C     E            =  YOUNG*S MODULUS
C     XNU          =  POISSON*S RATIO
C     XKP          =  PRESSURE DEPENDENT FLEXIBILITY FACTOR
C     AREA         =  SECTION AREA
C     XMI          =  MOMENT OF INERTIA
C     T            =  ANGLE OF THE BEND, THETA
C     ST           =  SIN(THETA)
C     CT           =  COS(THETA)
C     NODE         =  NODE NUMBER AT END J OF THE BEND
C     NEL          =  PIPE ELEMENT NUMBER
C     MODEX        =  EXECUTION MODE
C                     (EQ.1, DATA CHECK)
C     F(6,6)       =  FLEXIBILITY MATRIX AT NODE J
C     R            =  RADIUS OF THE BEND
C     THERM        =  THERMAL EXPANSION COEFFICIENT
C     P            =  INTERNAL PIPE PRESSURE
C     WALL         =  PIPE WALL THICKNESS
C     DOUT         =  OUTSIDE DIAMETER OF THE PIPE
C     B            =  FREE END DEFLECTIONS AT NODE J DUE TO
C                     (1)  UNIFORM LOAD IN THE X(I) DIRECTION
C                     (2)  UNIFORM LOAD IN THE Y(I) DIRECTION
C                     (3)  UNIFORM LOAD IN THE Z(I) DIRECTION
C                     (4)  UNIFORM THERMAL EXPANSION (DT=1)
C                     (5)  P,   INTERNAL PRESSURE
C     H            =  FORCE TRANSFORMATION RELATING REACTIONS AT NODE I
C                     DUE TO UNIT LOADS AT NODE J
C     S            =  LOCAL BEND ELEMENT STIFFNESS MATRIX
C     FEF          =  FIXED END FORCES (ACTING ON THE NODES) DUE TO
C                     (1)  UNIFORM LOAD IN THE X(I) DIRECTION
C                     (2)  UNIFORM LOAD IN THE Y(I) DIRECTION
C                     (3)  UNIFORM LOAD IN THE Z(I) DIRECTION
C                     (4)  UNIFORM THERMAL EXPANSION (DT=1)
C                     (5)  P,   INTERNAL PRESSURE
C     XM           =  LUMPED MASS MATRIX
C     SA           =  STRESS-DISPLACEMENT TRANSFORMATION RELATING THE
C                     12 GLOBAL COMPONENTS OF DISPLACEMENT TO THE 6
C                     LOCAL COMPONENTS OF MEMBER LOADS LOCATED AT NODE
C                     I, MIDPOINT OF THE ARC AND AT NODE J.
C     FEFC         =  FIXED-END FORCE CORRECTIONS TO THE MEMBER LOADS
C                     DUE TO THE FIVE (5) TYPES OF ELEMENT LOADS
C     XMAS         =  MASS   PER UNIT LENGTH OF THE SECTION
C     DC           =  ARRAY OF DIRECTION COSINES WHICH TRANSFORMS LOCAL
C                     VECTORS TO GLOBAL VECTORS
      SUBROUTINE BENDKS
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON /PIPEC/ALFAV,E,XNU,XKP,T,NOAX,NODE,NEL,
     1              MODEX,R,THERM,P,AREA,XMI,WALL,DOUT,XMAS
      COMMON /EM/    IXX(14),S(12,12),RF(12,4),XM(12),SA(18,12),
     1               SF(18,4),FEF(12,5),FEFC(18,5),F(6,6),B(6,6),
     2               H(6,6),DC(3,3),IFILL2(4162)
      COMMON /ELPAR/ NPAR(14),IFILL1(10)
C
      DIMENSION      COL(6)
C
C     SET THE FACTOR FOR AXIAL DEFORMATIONS
C
C##   AXIAL = 1.0
      AXIAL = 1.0D0
C##   IF(NOAX.EQ.1) AXIAL = 0.0
      IF(NOAX.EQ.1) AXIAL = 0.0D0
C
C     SET THE FACTOR FOR SHEAR DEFORMATIONS (EQ.0, NEGLECT)
C
      XKAP = ALFAV
C##   IF(ALFAV.GT.99.99) XKAP = 0.0
      IF(ALFAV.GT.99.99D0) XKAP = 0.0D0
C
C     SET THE FLEXIBILITY FACTORS
C
      XK0 = XKP
      XK1 = XKP
C
C     COMPUTE THE MATERIAL FACTORS
C
C##   RE = 1.0/E
      RE = 1.0D0/E
C##   XNU1 = 1.0+XNU
      XNU1 = 1.0D0+XNU
C
C     COMPUTE SECTION PROPERTY CONSTANTS
C
      RA = AXIAL*R*RE/AREA
      RV = XKAP*XNU1*R*RE/AREA
C##   RT = 0.5*XNU1*R*RE/XMI
      RT = 0.5D0*XNU1*R*RE/XMI
C##   RB0 = 0.5*XK0*R*RE/XMI
      RB0 = 0.5D0*XK0*R*RE/XMI
      RB1 = XK1*R*RE/XMI
      R2 = R**2
C
C     COMPUTE COMMON TRIGONOMETRIC CONSTANTS
C
      ST =DSIN(T)
      CT =DCOS(T)
C##   S2T =DSIN(2.0*T)
      S2T =DSIN(2.0D0*T)
C##   C2T =DCOS(2.0*T)
      C2T =DCOS(2.0D0*T)
      T2 = T**2
C
C     FORM THE NODE FLEXIBILITY MATRIX AT NODE J REFERENCED TO THE
C     LOCAL (X,Y,Z) COORDINATE SYSTEM AT NODE I.
C
C     X - DIRECTION...  IN-PLANE TANGENT TO THE BEND AT NODE I AND
C                       DIRECTED TOWARD NODE J
C     Y - DIRECTION...  IN-PLANE AND DIRECTED RADIALLY INWARD TO THE
C                       CENTER OF CURVATURE
C     Z - DIRECTION...  OUT OF PLANE AND ORTHOGONAL TO X AND Y
C
      DO 50 I=1,6
      DO 50 K=I,6
C##   F(I,K) = 0.0
      F(I,K) = 0.0D0
   50 CONTINUE
C
C     A X I A L
C
C##   F(1,1) = F(1,1) + 0.25*RA*(2.0*T+S2T)
      F(1,1) = F(1,1) + 0.25D0*RA*(2.0D0*T+S2T)
C##   F(2,2) = F(2,2) + 0.25*RA*(2.0*T-S2T)
      F(2,2) = F(2,2) + 0.25D0*RA*(2.0D0*T-S2T)
C     N O T E   (COEFFICIENT CHANGE)
C##   F(1,2) = F(1,2) + 0.50*RA* ST**2
      F(1,2) = F(1,2) + 0.50D0*RA* ST**2
C
C     S H E A R
C
C##   F(1,1) = F(1,1) + 0.5*RV*(2.0*T-S2T)
      F(1,1) = F(1,1) + 0.5D0*RV*(2.0D0*T-S2T)
C##   F(2,2) = F(2,2) + 0.5*RV*(2.0*T+S2T)
      F(2,2) = F(2,2) + 0.5D0*RV*(2.0D0*T+S2T)
C##   F(3,3) = F(3,3) + 2.0*RV* T
      F(3,3) = F(3,3) + 2.0D0*RV* T
C     N O T E   (SIGN CHANGE)
      F(1,2) = F(1,2) -     RV* ST**2
C
C     T O R S I O N
C
C##   F(3,3) = F(3,3) + 0.5*RT*R2*(6.0*T+S2T-8.0*ST)
      F(3,3) = F(3,3) + 0.5D0*RT*R2*(6.0D0*T+S2T-8.0D0*ST)
C##   F(4,4) = F(4,4) + 0.5*RT*   (2.0*T+S2T)
      F(4,4) = F(4,4) + 0.5D0*RT* (2.0D0*T+S2T)
C##   F(5,5) = F(5,5) + 0.5*RT*   (2.0*T-S2T)
      F(5,5) = F(5,5) + 0.5D0*RT* (2.0D0*T-S2T)
      F(3,4) = F(3,4) +     RT*R *(ST-T*CT)
C##   F(3,5) = F(3,5) +     RT*R *(2.0-2.0*CT-T*ST)
      F(3,5) = F(3,5) +     RT*R *(2.0D0-2.0D0*CT-T*ST)
C##   F(4,5) = F(4,5) + 0.5*RT*   (1.0-C2T)
      F(4,5) = F(4,5) + 0.5D0*RT* (1.0D0-C2T)
C
C     B E N D I N G
C
C##   F(1,1) = F(1,1) + 0.25*RB1*R2*(2.0*T*(2.0+C2T)-3.0*S2T)
      F(1,1) = F(1,1) + 0.25D0*RB1*R2*(2.0D0*T*(2.0D0+C2T)-3.0D0*S2T)
C##   F(2,2) = F(2,2) + 0.25*RB1*R2*(2.0*T*(2.0-C2T)+3.0*S2T-8.0*ST)
      F(2,2) = F(2,2) + 0.25D0*RB1*R2*(2.0D0*T*(2.0D0-C2T)+3.0D0*S2T-8.0
     *D0*ST)
C##   F(3,3) = F(3,3) + 0.50*RB0*R2*(2.0*T-S2T)
      F(3,3) = F(3,3) + 0.50D0*RB0*R2*(2.0D0*T-S2T)
C##   F(4,4) = F(4,4) + 0.50*RB0*   (2.0*T-S2T)
      F(4,4) = F(4,4) + 0.50D0*RB0* (2.0D0*T-S2T)
C##   F(5,5) = F(5,5) + 0.50*RB0*   (2.0*T+S2T)
      F(5,5) = F(5,5) + 0.50D0*RB0* (2.0D0*T+S2T)
      F(6,6) = F(6,6) +      RB1*T
C##   F(1,2) = F(1,2) + 0.25*RB1*R2*(1.0+3.0*C2T+2.0*T*S2T-4.0*CT)
      F(1,2) = F(1,2) + 0.25D0*RB1*R2*(1.0D0+3.0D0*C2T+2.0D0*T*S2T-4.0D0
     **CT)
      F(1,6) = F(1,6) -      RB1*R *(ST-T*CT)
C##   F(2,6) = F(2,6) +      RB1*R *(T*ST+CT-1.0)
      F(2,6) = F(2,6) +      RB1*R *(T*ST+CT-1.0D0)
      F(3,4) = F(3,4) +      RB0*R *(ST-T*CT)
      F(3,5) = F(3,5) -      RB0*R *T*ST
C##   F(4,5) = F(4,5) - 0.50*RB0*   (1.0-C2T)
      F(4,5) = F(4,5) - 0.50D0*RB0* (1.0D0-C2T)
C
      DO 60 I=1,6
      DO 60 K=I,6
      F(K,I) = F(I,K)
   60 CONTINUE
C**** PRINT THE NODE FLEXIBILITY MATRIX
      IF(NPAR(10).LT.1) GO TO 6700
      WRITE (6,4000)
      WRITE (6,4010) ((F(I,K),K=1,6),I=1,6)
 6700 CONTINUE
C****
C
C     FORM THE NODE STIFFNESS MATRIX
C
      CALL PINVER (F,6,6,NODE,NEL,MODEX)
C**** PRINT THE NODE STIFFNESS MATRIX
      IF(NPAR(10).LT.1) GO TO 6701
      WRITE (6,4020)
      WRITE (6,4030) ((F(I,K),K=1,6),I=1,6)
 6701 CONTINUE
C****
C
C     COMPUTE THE DEFLECTIONS/ROTATIONS (MEASURED IN THE X,Y,Z SYSTEM
C     AT NODE I) AT NODE J DUE TO UNIFORM LOADS IN EACH OF THE X,Y,Z
C     DIRECTIONS (AT I).  THE UNIFORM LOADS ARE DIRECTION INVARIANT
C     WITH POSITION ALONG THE ARC, AND NODE I IS FIXED WHILE NODE J IS
C     COMPLETELY FREE.
C
      DO 70 I=1,6
      DO 70 K=1,3
C##   B(I,K) = 0.0
      B(I,K) = 0.0D0
   70 CONTINUE
C
C     A X I A L
C
C##   RA = 0.125*RA*R
      RA = 0.125D0*RA*R
C##   B(1,1) = B(1,1) +     RA*(2.0*T2-C2T+1.0)
      B(1,1) = B(1,1) +     RA*(2.0D0*T2-C2T+1.0D0)
C##   B(2,2) = B(2,2) +     RA*(2.0*T2+C2T-1.0)
      B(2,2) = B(2,2) +     RA*(2.0D0*T2+C2T-1.0D0)
C     N O T E   (COEFFICIENT CHANGE)
C##   B(1,2) = B(1,2) +     RA*(2.0*T-S2T)
      B(1,2) = B(1,2) +     RA*(2.0D0*T-S2T)
C     N O T E   (COEFFICIENT CHANGE)
C##   B(2,1) = B(2,1) +     RA*(2.0*T-S2T)
      B(2,1) = B(2,1) +     RA*(2.0D0*T-S2T)
C
C     S H E A R
C
C##   RV = 0.25*RV*R
      RV = 0.25D0*RV*R
C##   B(1,1) = B(1,1) +     RV*(2.0*T2+C2T-1.0)
      B(1,1) = B(1,1) +     RV*(2.0D0*T2+C2T-1.0D0)
C##   B(2,2) = B(2,2) +     RV*(2.0*T2-C2T+1.0)
      B(2,2) = B(2,2) +     RV*(2.0D0*T2-C2T+1.0D0)
C##   B(3,3) = B(3,3) + 4.0*RV*T2
      B(3,3) = B(3,3) + 4.0D0*RV*T2
C     N O T E   (SIGN CHANGE)
C##   B(1,2) = B(1,2) -     RV*(2.0*T-S2T)
      B(1,2) = B(1,2) -     RV*(2.0D0*T-S2T)
C     N O T E   (SIGN CHANGE)
C##   B(2,1) = B(2,1) -     RV*(2.0*T-S2T)
      B(2,1) = B(2,1) -     RV*(2.0D0*T-S2T)
C
C     T O R S I O N
C
      RT = RT*R2
C##   B(3,3) = B(3,3) + 0.5*RT*R*(1.0+2.0*T2-4.0*T*ST-C2T)
      B(3,3) = B(3,3) + 0.5D0*RT*R*(1.0D0+2.0D0*T2-4.0D0*T*ST-C2T)
C##   B(4,3) = B(4,3) +     RT*  (2.0-2.0*CT-T*ST)
      B(4,3) = B(4,3) +     RT*  (2.0D0-2.0D0*CT-T*ST)
C##   B(5,3) = B(5,3) +     RT*  (T*(2.0+CT)-3.0*ST)
      B(5,3) = B(5,3) +     RT*  (T*(2.0D0+CT)-3.0D0*ST)
C
C     B E N D I N G
C
      RB0 = RB0*R2
      RB1 = RB1*R2
C##   B(1,1) = B(1,1) + 0.125*RB1*R*(7.0+2.0*T2+9.0*C2T+4.0*T*S2T
C##  1                               -16.0*CT)
      B(1,1) = B(1,1) + 0.125D0*RB1*R*(7.0D0+2.0D0*T2+9.0D0*C2T+4.0D0*T*
     *S2T                            -16.0D0*CT)
C##   B(2,2) = B(2,2) + 0.125*RB1*R*(1.0+2.0*T2-9.0*C2T-4.0*T*S2T
C##  1                               +8.0*(CT-T*ST))
      B(2,2) = B(2,2) + 0.125D0*RB1*R*(1.0D0+2.0D0*T2-9.0D0*C2T-4.0D0*T*
     *S2T                            +8.0D0*(CT-T*ST))
C##   B(3,3) = B(3,3) + 0.500*RB0*R*(3.0+C2T-4.0*CT)
      B(3,3) = B(3,3) + 0.500D0*RB0*R*(3.0D0+C2T-4.0D0*CT)
C##   B(1,2) = B(1,2) + 0.125*RB1*R*(9.0*S2T-4.0*T*(C2T+2.0*CT)-6.0*T)
      B(1,2) = B(1,2) + 0.125D0*RB1*R*(9.0D0*S2T-4.0D0*T*(C2T+2.0D0*CT)-
     *6.0D0*T)
C##   B(2,1) = B(2,1) + 0.125*RB1*R*(9.0*S2T-4.0*T*C2T-24.0*ST+10.0*T)
      B(2,1) = B(2,1) + 0.125D0*RB1*R*(9.0D0*S2T-4.0D0*T*C2T-24.0D0*ST+1
     *0.0D0*T)
C##   B(4,3) = B(4,3) +       RB0*  (2.0-2.0*CT-T*ST)
      B(4,3) = B(4,3) +       RB0*  (2.0D0-2.0D0*CT-T*ST)
      B(5,3) = B(5,3) -       RB0*  (ST-T*CT)
C##   B(6,1) = B(6,1) -       RB1*  (2.0-2.0*CT-T*ST)
      B(6,1) = B(6,1) -       RB1*  (2.0D0-2.0D0*CT-T*ST)
C##   B(6,2) = B(6,2) +       RB1*  (2.0*ST-T-T*CT)
      B(6,2) = B(6,2) +       RB1*  (2.0D0*ST-T-T*CT)
C
C     COMPUTE THE FREE NODE DEFLECTIONS AT END J DUE TO A UNIFORM
C     THERMAL EXPANSION
C
      DO 80 I=1,6
C##   B(I,4) = 0.0
      B(I,4) = 0.0D0
   80 CONTINUE
C
      DUM = R*THERM
      B(1,4) = DUM*ST
C##   B(2,4) = DUM*(1.0-CT)
      B(2,4) = DUM*(1.0D0-CT)
C
C     COMPUTE THE FREE NODE DEFLECTIONS AT END J DUE TO PRESSURE
C
      DO 90 I=1,6
C##   B(I,5) = 0.0
      B(I,5) = 0.0D0
   90 CONTINUE
C
C     COMPUTE THE ANGLE CHANGE AND END DISPLACEMENTS AT THE FREE END
C     OF THE BEND DUE TO INTERNAL PRESSURE, P.
C
C##   RM = (DOUT-WALL)*0.5
      RM = (DOUT-WALL)*0.5D0
      KK = 1
      GO TO (92,94),KK
C
C     MEL REPORT 10-66, EQUATION (3-29).
C
   92 CONTINUE
C##   DUM = 3.14159265*RM**4*P*T
      DUM = 3.14159265D0*RM**4*P*T
C##   DUM = 0.5*DUM*RE/XMI
      DUM = 0.5D0*DUM*RE/XMI
      DU2 = RM/R
C##   BTA = DUM*(2.0-2.0*XNU + (3.0+1.5*XNU)*DU2**2)
      BTA = DUM*(2.0D0-2.0D0*XNU + (3.0D0+1.5D0*XNU)*DU2**2)
      GO TO 96
C
C     C. S. PARKER, EQUATION (10), 2-28-69.
C
   94 CONTINUE
      DU2 = R/RM
C##   DUM = P*RM*0.5*RE/WALL
      DUM = P*RM*0.5D0*RE/WALL
C##   DU3 = 1.0 + DUM*(1.0-XNU*(2.0*DU2-1.0)/(DU2-1.0))
      DU3 = 1.0D0 + DUM*(1.0D0-XNU*(2.0D0*DU2-1.0D0)/(DU2-1.0D0))
C##   BTA = DU3/(1.0 + DUM*(2.0-XNU))
      BTA = DU3/(1.0D0 + DUM*(2.0D0-XNU))
C##   BTA = T*(1.0-BTA)
      BTA = T*(1.0D0-BTA)
C
   96 CONTINUE
      DUM = R/T*BTA
      B(1,5) = DUM*(ST-T*CT)
C##   B(2,5) = DUM*(1.0-CT-T*ST)
      B(2,5) = DUM*(1.0D0-CT-T*ST)
      B(6,5) =-BTA
C
C     AXIAL GROWTH DUE TO PRESSURE.  MEL REPORT 10-66, EQUATION (3-28).
C
C##   DUM = 0.5* P* RM* RE* (1.0-2.0*XNU)* R/ WALL
      DUM = 0.5D0* P* RM* RE* (1.0D0-2.0D0*XNU)* R/ WALL
      B(1,5) = B(1,5) + DUM* ST
C##   B(2,5) = B(2,5) + DUM* (1.0-CT)
      B(2,5) = B(2,5) + DUM* (1.0D0-CT)
C**** PRINT THE FREE END DEFLECTIONS
      IF(NPAR(10).LT.1) GO TO 6702
      WRITE (6,4050)
      WRITE (6,4060) ((B(I,K),K=1,5),I=1,6)
 6702 CONTINUE
C****
C
C     SET UP THE FORCE TRANSFORMATION RELATING REACTIONS AT NODE I
C     ACTING ON THE MEMBER END DUE TO UNIT LOADS APPLIED TO THE MEMBER
C     END AT NODE J.
C
      DO 100 I=1,6
      DO 100 K=1,6
C##   H(I,K) = 0.0
      H(I,K) = 0.0D0
  100 CONTINUE
C
      DO 105 K=1,6
C##   H(K,K) =-1.0
      H(K,K) =-1.0D0
  105 CONTINUE
C
C##   H(4,3) =-R*(1.0-CT)
      H(4,3) =-R*(1.0D0-CT)
      H(5,3) = R* ST
      H(6,1) =-H(4,3)
      H(6,2) =-H(5,3)
C
C     FORM THE UPPER TRIANGULAR PORTION OF THE LOCAL ELEMENT STIFFNESS
C     MATRIX FOR THE BEND
C
      DO 110 K=1,6
      DO 110 I=K,6
      S(K+6,I+6) = F(K,I)
  110 CONTINUE
C
      DO 130 IR=1,6
      DO 130 IC=1,6
C##   S(IR,IC+6) = 0.0
      S(IR,IC+6) = 0.0D0
      DO 120 IN=1,6
      S(IR,IC+6) = S(IR,IC+6) + H(IR,IN)* F(IN,IC)
  120 CONTINUE
  130 CONTINUE
C
      DO 150 IR=1,6
      DO 150 IC=IR,6
C##   S(IR,IC) = 0.0
      S(IR,IC) = 0.0D0
      DO 140 IN=1,6
      S(IR,IC) = S(IR,IC) + S(IR,IN+6)* H(IC,IN)
  140 CONTINUE
  150 CONTINUE
C
C     REFLECT FOR SYMMETRY
C
      DO 160 I=1,12
      DO 160 K=I,12
      S(K,I) = S(I,K)
  160 CONTINUE
C**** PRINT THE BEND LOCAL STIFFNESS MATRIX
      IF(NPAR(10).LT.1) GO TO 6703
      WRITE (6,4500)
      WRITE (6,4510) ((S(I,J),J=1,6 ),I=1,12)
      WRITE (6,4510) ((S(I,J),J=7,12),I=1,12)
 6703 CONTINUE
C****
C
C     COMPUTE THE RESTRAINED NODE FORCES ACTING ON THE NODES OF THE
C     BEND DUE TO THE MEMBER LOADINGS
C
      DO 180 I=1,5
      DO 180 J=1,12
C##   FEF(J,I) = 0.0
      FEF(J,I) = 0.0D0
      DO 170 K=1,6
      FEF(J,I) = FEF(J,I) - S(J,K+6)* B(K,I)
  170 CONTINUE
  180 CONTINUE
C
C     FOR THE DISTRIBUTED LOADS SUPERIMPOSE THE CANTILEVER REACTIONS
C     ACTING ON THE ELEMENT AT NODE I.
C
      FEF(1,1) = FEF(1,1) - R*T
      FEF(6,1) = FEF(6,1) + R2*(T-ST)
C
      FEF(2,2) = FEF(2,2) - R*T
C##   FEF(6,2) = FEF(6,2) - R2*(1.0-CT)
      FEF(6,2) = FEF(6,2) - R2*(1.0D0-CT)
C
      FEF(3,3) = FEF(3,3) - R*T
      FEF(4,3) = FEF(4,3) - R2*(T-ST)
C##   FEF(5,3) = FEF(5,3) + R2*(1.0-CT)
      FEF(5,3) = FEF(5,3) + R2*(1.0D0-CT)
C**** PRINT THE FIXED END QUANTITIES
      IF(NPAR(10).LT.1) GO TO 6704
      WRITE (6,4600)
      WRITE (6,4610) ((FEF(I,J),J=1,5),I=1,12)
 6704 CONTINUE
C****
C
C     FORM THE LUMPED MASS MATRIX
C
C##   DUM = 0.5*R*T*XMAS
      DUM = 0.5D0*R*T*XMAS
      DO 200 K=1,3
      XM(K)   = DUM
      XM(K+6) = DUM
C##   XM(K+3) = 0.0
      XM(K+3) = 0.0D0
C##   XM(K+9) = 0.0
      XM(K+9) = 0.0D0
  200 CONTINUE
C
C     COMPUTE THE FIXED-NODE CORRECTIONS TO THE MEMBER LOADS RESULTING
C     FROM ELEMENT LOADINGS.  FORCES ACT ON THE SEGMENT BETWEEN THE
C     POINT WHERE EVALUATED AND NODE I.
C
C        1. AT NODE I (ACTING ON NODE I)
C
      DO 210 I=1,5
      DO 210 K=1,6
      FEFC(K,I) = -FEF(K,I)
  210 CONTINUE
C
C        2. AT NODE J (ROTATE IN-PLANE FROCES AN AMOUNT THETA)
C
      DO 220 I=1,5
      DO 215 K=1,4,3
      FEFC(K+12,I) =  CT* FEF(K+6,I) + ST* FEF(K+7,I)
      FEFC(K+13,I) = -ST* FEF(K+6,I) + CT* FEF(K+7,I)
      FEFC(K+14,I) =      FEF(K+8,I)
  215 CONTINUE
  220 CONTINUE
C
C        3. AT THE MIDPOINT OF THE ARC BETWEEN NODES I AND J.
C
C           A. TRANSFER FORCES AT NODE J TO THE MIDPOINT AND ROTATE
C              AN AMOUNT THETA/2
C
C##   S12T =DSIN(0.5*T)
      S12T =DSIN(0.5D0*T)
C##   C12T =DCOS(0.5*T)
      C12T =DCOS(0.5D0*T)
      DX   = R*(ST-S12T)
      DY   = R*(C12T-CT)
C
      DO 230 I=1,5
      XM10 = FEF(10,I) + FEF(9,I)* DY
      XM11 = FEF(11,I) - FEF(9,I)* DX
      FEFC( 7,I) =  C12T* FEF(7,I) + S12T* FEF(8,I)
      FEFC( 8,I) = -S12T* FEF(7,I) + C12T* FEF(8,I)
      FEFC( 9,I) =  FEF(9,I)
      FEFC(10,I) =  C12T* XM10     + S12T* XM11
      FEFC(11,I) = -S12T* XM10     + C12T* XM11
  230 FEFC(12,I) =  FEF(12,I) - FEF(7,I)* DY + FEF(8,I)* DX
C
C           B. FOR THE DISTRIBUTED LOADS SUPERIMPOSE THE RESULTANT
C              OF THE APPLIED LOADING TRANSFERRED TO THE MIDPOINT OF
C              THE ARC AND ROTATE AN AMOUNT THETA/2 (IN-PLANE)
C
C##   DDX = R*(2.0*(C12T-CT)/T - S12T)
      DDX = R*(2.0D0*(C12T-CT)/T - S12T)
C##   DDY = R*(2.0*(S12T-ST)/T + C12T)
      DDY = R*(2.0D0*(S12T-ST)/T + C12T)
C##   DUM = R*T*0.5
      DUM = R*T*0.5D0
C
      FEFC( 7,1) = FEFC( 7,1) + C12T* DUM
      FEFC( 8,1) = FEFC( 8,1) - S12T* DUM
      FEFC(12,1) = FEFC(12,1) - DDY * DUM
C
      FEFC( 7,2) = FEFC( 7,2) + S12T* DUM
      FEFC( 8,2) = FEFC( 8,2) + C12T* DUM
      FEFC(12,2) = FEFC(12,2) + DDX * DUM
C
      FEFC( 9,3) = FEFC( 9,3) + DUM
      XM10 =  DDY* DUM
      XM11 = -DDX* DUM
      FEFC(10,3) = FEFC(10,3) + C12T* XM10 + S12T* XM11
      FEFC(11,3) = FEFC(11,3) - S12T* XM10 + C12T* XM11
C**** PRINT THE FIXED-END CORRECTIONS
      IF(NPAR(10).LT.1) GO TO 6705
      WRITE (6,4650)
      WRITE (6,4660) ((FEFC(I,J),J=1,5),I=1,18)
 6705 CONTINUE
C****
C
C     FORM THE TRANSFORMATION RELATING GLOBAL DISPLACEMENTS AND MEMBER
C     FORCES AT NODE I, MIDPOINT AND NODE J.
C
C        1. STRESS RESULTANTS AT NODE I
C
      DO 260 K1=1,10,3
      NRS = K1-1
      DO 260 K2=1,10,3
      NCS = K2-1
      DO 250 IR=1,3
      NR = NRS+IR
      DO 250 IC=1,3
      NC = NCS+IC
C##   SA(NR,NC) = 0.0
      SA(NR,NC) = 0.0D0
      DO 240 IN=1,3
      N = NCS+IN
      SA(NR,NC) = SA(NR,NC) - S(NR,N)* DC(IC,IN)
  240 CONTINUE
  250 CONTINUE
  260 CONTINUE
C
C        2. STRESS RESULTANTS AT NODE J
C
      H(1,1) = CT
      H(1,2) = ST
      H(2,1) =-ST
      H(2,2) = CT
C##   H(3,3) = 1.0
      H(3,3) = 1.0D0
C
      DO 290 K1=7,10,3
      NRS = K1-1
      DO 280 IR=1,3
      NR = NRS+IR
      DO 280 IC=1,12
C##   SA(NR+6,IC) = 0.0
      SA(NR+6,IC) = 0.0D0
      DO 270 IN=1,3
      N = NRS+IN
      SA(NR+6,IC) = SA(NR+6,IC) - H(IR,IN)* SA(N,IC)
  270 CONTINUE
  280 CONTINUE
  290 CONTINUE
C
C        3. STRESS RESULTANTS AT THE MIDPOINT OF THE ARC
C
      H(1,1) =  C12T
      H(1,2) =  S12T
      H(2,1) = -S12T
      H(2,2) =  C12T
C##   H(3,3) =  1.0
      H(3,3) =  1.0D0
      DO 300 I=1,3
      DO 300 K=1,3
  300 H(I+3,K+3) = H(I,K)
      H(4,3) =  DY* C12T - DX* S12T
      H(5,3) = -DY* S12T - DX* C12T
      H(6,1) = -DY
      H(6,2) =  DX
C
      DO 320 IC=1,12
      DO 310 N=1,6
  310 COL(N) = SA(N+6,IC)
      DO 320 IR=1,6
C##   SA(IR+6,IC) = 0.0
      SA(IR+6,IC) = 0.0D0
      DO 315 IN=1,6
      SA(IR+6,IC) = SA(IR+6,IC) - H(IR,IN)* COL(IN)
  315 CONTINUE
  320 CONTINUE
C**** PRINT THE STRESS DISPLACEMENT TRANSFORMATION
      IF(NPAR(10).LT.1) GO TO 6706
      WRITE (6,4700)
      WRITE (6,4710) ((SA(I,J),J=1, 6),I=1,18)
      WRITE (6,4710) ((SA(I,J),J=7,12),I=1,18)
 6706 CONTINUE
C****
C
 4000 FORMAT (/// 24H NODE FLEXIBILITY MATRIX, // 1X)
 4010 FORMAT ( 1X / (6E20.8) )
 4020 FORMAT (/// 22H NODE STIFFNESS MATRIX, // 1X)
 4030 FORMAT ( 1X / (6E20.8) )
 4050 FORMAT (/// 42H FREE NODE DISPLACEMENTS  (5 MEMBER LOADS), // 1X)
 4060 FORMAT (1X / (5E20.8) )
 4500 FORMAT (23H LOCAL STIFFNESS MATRIX, // 1X)
 4510 FORMAT (// (/6E15.6) )
 4600 FORMAT (// 17H FIXED END FORCES, // 1X)
 4610 FORMAT (5E20.8)
 4650 FORMAT (// 43H STRESS CORRECTIONS DUE TO FIXED END FORCES, // 1X)
 4660 FORMAT (5E20.8)
 4700 FORMAT (//35H STRESS-DISPLACEMENT TRANSFORMATION, / 1X)
 4710 FORMAT (/// (6E20.8) )
C
      RETURN
      END
      SUBROUTINE BOUND
      IMPLICIT REAL*8(A-H,O-Z)
C
C     CALLS:  CLAMP,STRSC
C     CALLED BY:  ELTYPE
C
C---      COMMON A(1)
      COMMON A(100000001)
      COMMON /ELPAR/ NPAR(14),NUMNP,MBAND,NELTYP,N1,N2,N3,N4,N5,MTOT,NEQ
      COMMON /JUNK/  LT,LH,L,IPAD,SIG(20),IFILL(440)
      COMMON /EXTRA/ MODEX,NT8,N10SV,NT10,IFILL2(22)
C
      IF (NPAR(1).EQ.0) GO TO 500
C
      CALL CLAMP (NPAR(2),A(N1),A(N2),A(N3),A(N4),NUMNP,MBAND)
C
      RETURN
C
  500 WRITE (6,2002)
      NUME=NPAR(2)
      DO 800 MM=1,NUME
      CALL STRSC (A(N1),A(N3),NEQ,0)
      WRITE (6,2001)
      DO 800 L=LT,LH
      CALL STRSC (A(N1),A(N3),NEQ,1)
      WRITE (6,3002) MM,L,(SIG(I),I=1,2)
C***  STRESS PORTHOLE
      IF(N10SV.EQ.1)
     *WRITE (NT10) MM,L,SIG(1),SIG(2)
  800 CONTINUE
      RETURN
C
 2001 FORMAT (/)
 2002 FORMAT (48H B O U N D A R Y   E L E M E N T   F O R C E S /,
     1        14H M O M E N T S, // 8H ELEMENT,3X,4HLOAD,14X,5HFORCE,
     2        9X,6HMOMENT, / 8H  NUMBER,3X,4HCASE, // 1X)
 3002 FORMAT (I8,I7,4X,2E15.5)
      END
      SUBROUTINE BRICK8 (S,STR,NBRK8,NMAT,NLD,ID,X,Y,Z,T,EE,ENU,RHO,
     .                   ALPT,KTYPE,PR,YREF,NFACE,NUMNP)
      IMPLICIT REAL*8(A-H,O-Z)
C
C     CALLS:  DERIV,LOAD,LOSTR,CALBAN
C     CALLED BY:  THREED
C
C     STIFFNESS SUBROUTINE FOR 24 D.F. ISOPARAMETRIC HEXAHEDRON
C     LINEAR ELASTIC ISOTROPIC MATERIAL
C     'NINT*NINT*NINT' GAUSSIAN INTEGRATION RULE USED (NINT=1,2,3,4)
C
      DIMENSION KTYPE(1),PR(1),YREF(1),NFACE(1)
      DIMENSION T(1)
      DIMENSION X(1),Y(1),Z(1),ID(NUMNP,6)
      COMMON/EM/LM(24),ND,NS, SS(24,24),RF(24,4),XM(24),SA(12,24),
     .    SF(12,4),IFILL2(3604)
C---      EQUIVALENCE (IS1,SF(4)) , (IS2,SF(6))
      EQUIVALENCE (IS1,SF(4,1)) , (IS2,SF(6,1))
      DIMENSION EE(1),ENU(1),RHO(1),ALPT(1)
      COMMON /GASS/ XK(4,4),WGT(4,4),IPERM(3)
      COMMON /JUNK/ E1,E2,E3,DET,MLD(4),KLD(4),MULT(4),NP(8),INP(8),
     .              A(3,3),P(3,11),B(3,3),XX(8,3),Q(11),DL(8),
     .              TT(24),XLF(4),YLF(4),ZLF(4),TLF(4),PLF(4),
     .              REFT,INEL,ININT,IMAT,IINC,TTEMP,NEL,ML,NINT,MAT,
     .              INC,IPAD,TAG,TEMP,SKIP,I,J,K,L,FAC,CC1,CC2,CC3,CC4,
     .              G,DEN,FACT,GT,GG,C1,C2,C3,C,K1,K2,IFILL1(118)
      COMMON /ELPAR/ NPAR(14),NUMN ,MBAND,NELTYP,N1,N2,N3,N4,N5,MTOT,NEQ
      COMMON /EXTRA/ MODEX,NT8,IFILL3(24)
      DIMENSION S(33,33),STR(12,33)
      DIMENSION E(3,3)
      DIMENSION IS(2),ISP(2)
C
C     DATA STAR/'*'/,BLNK/' '/
      DATA STAR,BLNK /1H*,1H /
      DIMENSION STPTS(7,3)
C##   DATA STPTS / 0. , 1. ,-1. , 0. , 0. , 0. , 0. ,
C##  .             0. , 0. , 0. , 1. ,-1. , 0. , 0. ,
C##  .             0. , 0. , 0. , 0. , 0. , 1. ,-1. /
      DATA STPTS / 0.D0 , 1.D0 ,-1.D0 , 0.D0 , 0.D0 , 0.D0 , 0.D0 ,
     *             0.D0 , 0.D0 , 0.D0 , 1.D0 ,-1.D0 , 0.D0 , 0.D0 ,
     *             0.D0 , 0.D0 , 0.D0 , 0.D0 , 0.D0 , 1.D0 ,-1.D0 /
      XK(1,1)=0.0D0
      XK(2,1)=0.0D0
      XK(3,1)=0.0D0
      XK(4,1)=0.0D0
      XK(1,2)=-.5773502691896D0
      XK(2,2)=-XK(1,2)
      XK(3,2)=0.0D0
      XK(4,2)=0.0D0
      XK(1,3)=-.7745966692415D0
      XK(2,3)=0.0D0
      XK(3,3)=-XK(1,3)
      XK(4,3)=0.0D0
      XK(1,4)=-.8611363115941D0
      XK(2,4)=-.3399810435849D0
      XK(3,4)=-XK(2,4)
      XK(4,4)=-XK(1,4)
      WGT(1,1)=2.0D0
      WGT(2,1)=0.0D0
      WGT(3,1)=0.0D0
      WGT(4,1)=0.0D0
      WGT(1,2)=1.0D0
      WGT(2,2)=1.0D0
      WGT(3,2)=0.0D0
      WGT(4,2)=0.0D0
      WGT(1,3)=.55555555555556D0
      WGT(2,3)=.8888888888889D0
      WGT(3,3)=.5555555555556D0
      WGT(4,3)=0.0D0
      WGT(1,4)=.3478548451375D0
      WGT(2,4)=.6521451548625D0
      WGT(3,4)=WGT(2,4)
      WGT(4,4)=WGT(1,4)
      IPERM(1)=2
      IPERM(2)=3
      IPERM(3)=1
C
C
C
C     ZERO EM
C
      WRITE (6,3000) NBRK8,NMAT,NLD
      DO 9 I=1,1058
    9 LM(I)=0
C
C     MATERIAL PROPERTIES
C
      WRITE (6,1300)
      DO 1 I=1,NMAT
      READ  (5,1001) N,EE(N),ENU(N),RHO(N),ALPT(N)
    1 WRITE (6,2001) N,EE(N),ENU(N),RHO(N),ALPT(N)
C***  DATA PORTHOLE SAVE
      IF(MODEX.EQ.1)
     *WRITE (NT8) (EE(I),ENU(I),RHO(I),ALPT(I),I=1,NMAT)
C
C     ELEMENT DISTRIBUTED LOAD CARDS
C
      IF(NLD) 23,23,15
   15 WRITE (6,1302)
      DO 16 I=1,NLD
      READ  (5,1002) N,KTYPE(N),PR(N),YREF(N),NFACE(N)
   16 WRITE (6,2002) N,KTYPE(N),PR(N),YREF(N),NFACE(N)
C***  DATA PORTHOLE SAVE
      IF(MODEX.EQ.1)
     *WRITE (NT8) (KTYPE(N),PR(N),YREF(N),NFACE(N),N=1,NLD)
C
   23 READ  (5,1003) GRAV,PLF,TLF,XLF,YLF,ZLF
      WRITE (6,2003) GRAV,PLF,TLF,XLF,YLF,ZLF
C##   IF(GRAV.EQ.0.) GRAV=1.E+10
      IF(GRAV.EQ.0.D0) GRAV=1.D+10
C***  DATA PORTHOLE SAVE
      IF(MODEX.EQ.1)
     *WRITE (NT8) GRAV,PLF,TLF,XLF,YLF,ZLF
C
      WRITE (6,1301)
      NEL=0
   30 READ  (5,1000) INEL,INP,ININT,IMAT,IINC,MLD,ISP,TTEMP
      IF(IINC.EQ.0) IINC=1
      IF(IMAT.EQ.0) IMAT=1
   40 NEL=NEL+1
      ML=INEL-NEL
      IF(ML) 50,55,60
   50 WRITE (6,4003) INEL
      STOP
   55 DO 56 I=1,8
   56 NP(I)=INP(I)
      DO 39 I=1,4
   39 MULT(I)=1
      NINT=ININT
      MAT=IMAT
      INC=IINC
      TAG=STAR
      REFT=TTEMP
      IS(1)=ISP(1)
      IS(2)=ISP(2)
C##   SKIP=99999.
      SKIP=99999.D0
      IF(NINT) 33,33,57
   33 NINT=IABS(NINT)
C##   SKIP=1.
      SKIP=1.D0
C##   IF(NINT.EQ.0) SKIP=0.
      IF(NINT.EQ.0) SKIP=0.D0
   57 CONTINUE
      DO 59 I=1,4
      KLD(I)=IABS(MLD(I))
      IF(MLD(I)) 58,58,59
   58 MULT(I)=0
   59 CONTINUE
      GO TO 62
C
   60 DO 61 I=1,8
   61 NP(I)=NP(I)+INC
      TAG=BLNK
      DO 64 I=1,4
   64 KLD(I)=KLD(I)*MULT(I)
C
   62 IF(MODEX.EQ.1) GO TO 540
C##   TEMP = 0.0
      TEMP = 0.0D0
      DO 10 I=1,8
      K=NP(I)
      TEMP=TEMP+T(K)
      XX(I,1)=X(K)
      XX(I,2)=Y(K)
   10 XX(I,3)=Z(K)
C##   TEMP=TEMP*0.125
      TEMP=TEMP*0.125D0
      K=MAT
C##   FAC = EE(K)/((1.-2.*ENU(K))*(1.+ENU(K)))
      FAC = EE(K)/((1.D0-2.D0*ENU(K))*(1.D0+ENU(K)))
C##   FACT=FAC*ALPT(K)*(TEMP-REFT)*(1.+ENU(K))
      FACT=FAC*ALPT(K)*(TEMP-REFT)*(1.D0+ENU(K))
      IF(SKIP) 70,70,63
C##63 SKIP=SKIP-1.
   63 SKIP=SKIP-1.D0
C##   CC1=1.-ENU(K)
      CC1=1.D0-ENU(K)
      CC2=ENU(K)
C##   CC3=.5-ENU(K)
      CC3=.5D0-ENU(K)
C
      DO 100 I=1,33
      DO 100 J=1,33
  100 S(I,J)=0.0D0
      DO 110 I=1,24
C##10 TT(I)=0.
  110 TT(I)=0.D0
      DO 120 I=1,8
C##20 DL(I)=0.
  120 DL(I)=0.D0
C##   VOLUME = 0.0
      VOLUME = 0.0D0
C
C     LOOP OVER NINT**3 INTEGRATION POINTS
C
      DO 300  LX = 1,NINT
      E1=XK(LX,NINT)
      DO 300  LY = 1,NINT
      E2=XK(LY,NINT)
      DO 300  LZ = 1,NINT
      E3=XK(LZ,NINT)
C
      CALL DERIV(1,SA)
C
      GT= WGT(LX,NINT)*WGT(LY,NINT)*WGT(LZ,NINT)*DET
      VOLUME = VOLUME + GT
      GG=GT*RHO(MAT)
      G=GT*FAC
      C1=G*CC1
      C2=G*CC2
      C3=G*CC3
C
      L=0
      DO 310 I=1,8
      DL(I)=DL(I) + GG*Q(I)
      DO 310 K=1,3
      L=L+1
  310 TT(L)=TT(L) + GT*SA(I,K)
C
C     ADD CONTRIBUTION TO STIFFNESS MATRIX
C
      DO 300 I=1,11
      K3 = 3*I
      K2 = K3 - 1
      K1 = K2 - 1
      UI=SA(I,1)
      VI=SA(I,2)
      WI=SA(I,3)
      DO 300 J=I,11
      L3 = 3*J
      L2 = L3 - 1
      L1 = L2 - 1
      UJ=SA(J,1)
      VJ=SA(J,2)
      WJ=SA(J,3)
      UU=UI*UJ
      VV=VI*VJ
      WW=WI*WJ
      UV=UI*VJ
      VU=VI*UJ
      UW=UI*WJ
      WU=WI*UJ
      VW=VI*WJ
      WV=WI*VJ
      S(K1,L1) = S(K1,L1) + C1*UU + C3*(VV+WW)
      S(K2,L2) = S(K2,L2) + C1*VV + C3*(WW+UU)
      S(K3,L3) = S(K3,L3) + C1*WW + C3*(UU+VV)
      S(K1,L2) = S(K1,L2) + C2*UV + C3*VU
      S(K1,L3) = S(K1,L3) + C2*UW + C3*WU
      S(K2,L3) = S(K2,L3) + C2*VW + C3*WV
      IF (I.EQ.J)  GO TO 300
      S(K2,L1) = S(K2,L1) + C2*VU + C3*UV
      S(K3,L1) = S(K3,L1) + C2*WU + C3*UW
      S(K3,L2) = S(K3,L2) + C2*WV + C3*VW
  300 CONTINUE
C
C     FORM STRAIN MATRIX
C
      NSS=2
      IF(IS(2).EQ.0) NSS=1
      DO 305 I=1,12
      DO 305 J=1,33
C##05 STR(I,J)=0.
  305 STR(I,J)=0.D0
      DO 405 L=1,NSS
      LL=IS(L)+1
      E1=STPTS(LL,1)
      E2=STPTS(LL,2)
      E3=STPTS(LL,3)
      CALL DERIV(2,SA)
      L3=6*L-6
      DO 402 K=1,11
      K3=3*K
      K2=K3-1
      K1=K2-1
      STR(L3+1,K1) = SA(K,1)
      STR(L3+2,K2) = SA(K,2)
      STR(L3+3,K3) = SA(K,3)
      STR(L3+4,K1) = SA(K,2)
      STR(L3+4,K2) = SA(K,1)
      STR(L3+5,K2) = SA(K,3)
      STR(L3+5,K3) = SA(K,2)
      STR(L3+6,K1) = SA(K,3)
  402 STR(L3+6,K3) = SA(K,1)
  405 CONTINUE
      NS=6*NSS
C
C     STATIC CONDENSATION
C
      DO 710 M=1,9
      MN=34-M
      MO=MN-1
C      STIFFNESS MATRIX - S
      SP=S(MN,MN)
      DO 650 I=1,MO
  650 S(MN,I)=S(I,MN)/SP
      DO 700 K=1,MO
      SP=S(MN,K)
      DO 700 J=1,K
  700 S(J,K)=S(J,K) - SP*S(J,MN)
C      DERIVATIVE MATRIX - STR
      DO 710 J=1,NS
      SP=STR(J,MN)
C##   IF(SP.EQ.0.) GO TO 710
      IF(SP.EQ.0.D0) GO TO 710
      DO 705 K=1,MO
  705 STR(J,K)=STR(J,K) - SP*S(MN,K)
  710 CONTINUE
C
      DO 760 I=1,24
      DO 760 J=I,24
      SS(I,J)=S(I,J)
  760 SS(J,I)=SS(I,J)
C
C     STRAIN TO STRESS MATRIX
C
      E(1,1)=CC1*FAC
      E(2,2)=E(1,1)
      E(3,3)=E(1,1)
      E(1,2)=CC2*FAC
      E(1,3)=E(1,2)
      E(2,3)=E(1,2)
      E(2,1)=E(1,2)
      E(3,1)=E(1,2)
      E(3,2)=E(1,2)
C
      DO 900 I=1,NSS
      II=I*6-6
      DO 850 J=1,3
      DO 850 K=1,24
C##   SP=0.
      SP=0.D0
      DO 840 L=1,3
  840 SP = SP + E(J,L)*STR(II+L,K)
      SA(II+J,K)=SP
      JJ=II+3+J
  850 SA(JJ,K)=CC3*FAC*STR(JJ,K)
C
C
      DO 860 J=1,3
      JJ=J+3
      DO 860 K=1,4
      SF(II+J,K)=-FACT*TLF(K)
C##60 SF(II+JJ,K)=0.
  860 SF(II+JJ,K)=0.D0
C
      IF (IS(I).LE.0) GO TO 900
      LL=IS(I)+1
      E1=STPTS(LL,1)
      E2=STPTS(LL,2)
      E3=STPTS(LL,3)
      CALL DERIV (4,SA)
      CALL LOSTR (IS,A,B,SA,SF,I)
C
  900 CONTINUE
C
 70   CONTINUE
C
C     DISTRIBUTED LOAD
C
      DO 410 J=1,24
      DO 410 I=1,4
C##10 RF(J,I)=0.
  410 RF(J,I)=0.D0
      CALL LOAD (KTYPE,PR,YREF,NFACE)
C
C     SELF WQT.
C
      DO 460 II=1,8
      K=3*II
      J=K-1
      I=J-1
      DO 460 L=1,4
      RF(I,L) = RF(I,L)*PLF(L) + DL(II)*XLF(L)
      RF(J,L) = RF(J,L)*PLF(L) + DL(II)*YLF(L)
  460 RF(K,L) = RF(K,L)*PLF(L) + DL(II)*ZLF(L)
C
C     THERMAL LOADS
C
      DO 470 I=1,24
      GT=TT(I)*FACT
      DO 470 J=1,4
  470 RF(I,J)=RF(I,J) + GT*TLF(J)
C
C     MASS ARRAY
C
      L=0
C##   DUM=VOLUME*RHO(MAT)*.125/GRAV
      DUM=VOLUME*RHO(MAT)*.125D0/GRAV
      DO 465 I=1,8
      DO 465 J=1,3
      L=L+1
  465 XM(L) = DUM
C
  540 IJ = 0
      DO 550 I=1,8
      II=NP(I)
      DO 550 J=1,3
      IJ=IJ+1
  550 LM(IJ)=ID(II,J)
      ND=24
C
      IS1=IS(1)
      IS2=IS(2)
      NDM=24
      CALL CALBAN (MBAND,NDIF,LM,XM,SS,RF,ND,NDM,NS)
      IF(MODEX.EQ.1) GO TO 560
      WRITE (1) ND,NS,(LM(I),I=1,ND),((SA(I,J),I=1,NS),J=1,ND),
     1 ((SF(I,J),I=1,NS),J=1,4)
  560 IF(MODEX.EQ.1)
     *WRITE (NT8) NEL,NP,NINT,MAT,KLD,REFT,IS
      WRITE (6,2000) NEL,NP,NINT,MAT,TAG,KLD,REFT,IS,NDIF
C
C     CHECK IF LAST ELEMENT
C
      IF(NBRK8-NEL) 50,600,590
  590 IF(ML) 30,30,40
C
C
  600 RETURN
C
C
 1000 FORMAT (12I5,4I2,2I1,F10.2)
 1001 FORMAT (I5,4F10.0)
 1002 FORMAT (2I5,2F10.2,I5)
 1003 FORMAT(F10.2/(4F10.2))
 2000 FORMAT(I6,1X,8I5,I9,I12,8X,A1,3X,4I5,F9.1,5X,2I3,I8)
 2001 FORMAT(1X,I5,4E15.4)
 2002 FORMAT (I5,I9,2F13.3,I12)
 2003 FORMAT (/////
     . 35H .....ACCELERATION DUE TO GRAVITY = ,F10.2////
     . 38H LOAD FACTORS FOR 4 ELEMENT LOAD CASES  //
     . 46X,17HELEMENT LOAD CASE /
     . 36X,1HA,9X,1HB,9X,1HC,9X,1HD  /
     . 30H PRESSURE 
Comments: