$NOFLOATCALLS $NODEBUG C ................................................................... PROGRAM SHAKE91 CHARACTER*32 FIN,FOUT,PUNCH COMMON X(102520) COMMON /TIME/ T(9) COMMON /WGK/ WW, GT C ................................................................... WRITE(*,100) WRITE(*,101) 100 FORMAT( + 2X,'***********************************************************'/ + 2X,'* SHAKE -- A COMPUTER PROGRAM FOR EARTHQUAKE RESPONSE *'/ + 2X,'* ANALYSIS OF HORIZONTALLY LAYERED SITES *'/ + 2X,'* by: Per B. Schnabel & John Lysmer -- 1970 *'/ + 2X,'* ------------------------------------------------------- *'/ + 2X,'* shake85 IBM-PC version of SHAKE *'/ + 2X,'* by: S.S. (Willie) Lai, January 1985 *'/ + 2X,'* ------------------------------------------------------- *'/ + 2X,'* shake88 : New modulus reduction curves for clays added*'/ + 2X,'* using results from Sun et al (1988) *'/ + 2X,'* by: J. I. Sun & Ramin Golesorkhi *'/ + 2X,'* February 26, 1988 *'/ + 2X,'* ------------------------------------------------------- *'/ + 2X,'* SHAKE90/91: Adjust last iteration; Input now is either *'/ + 2X,'* Gmax or max Vs; up to 13 material types can *'/ + 2X,'* be specified by user; up to 50 Layers can *'/ + 2X,'* be specified; object motion can be read in *'/ + 2X,'* from a separate file and can have user *'/ + 2X,'* specified format; Different periods for *'/ + 2X,'* response spectral calculations; options *'/ + 2X,'* are renumbered; and general cleanup *'/ + 2X,'* by: J. I. Sun, I. M. Idriss & P. Dirrim *'/ + 2X,'* June 1990 - February 1991 *'/ + 2X,'* ------------------------------------------------------- *'/ + 2X,'* SHAKE91 : General cleanup and finalization of input/ *'/ + 2X,'* output format ... etc *'/ + 2X,'* by: I. M. Idriss *'/ + 2X,'* December 1991 *') 101 FORMAT(2X,'*----------------------------------------------------' + '---- *'/ + 2X,'* SHAKEVT : 1) Array sizes increased to read in ', + 'accel *'/ + 2X,'* time histories digitized at 0.005 sec. *'/ + 2X,'* 2) Algorithm used to compute shear stress *'/ + 2X,'* corrected to compute shear stress time *'/ + 2x,'* history for visco-elastic media. *'/ + 2x,'* 3) Option 8 changed to compute normalized *'/ + 2x,'* dissipated energy profile. *'/ + 2X,'* by: R. Green *'/ + 2X,'* February 2000 *'/ + 2X,'***********************************************************') C WRITE(*,200) 200 FORMAT(4X,'Name of Input File ='\) READ(*,10) FIN C WRITE(*,300) 300 FORMAT(4X,'Name of Output File #1 (input, peak values .. etc) ='\) READ(*,10) FOUT C WRITE(*,400) 400 FORMAT(4X,'Name of Output File #2 (time histories .. etc) ='\) READ(*,10) PUNCH C 10 FORMAT(A32) C OPEN(5,FILE=FIN,STATUS='OLD') OPEN(6,FILE=FOUT,STATUS='NEW') OPEN(7,FILE=PUNCH,STATUS='NEW') C WRITE(6,500) FIN,FOUT,PUNCH WRITE(7,500) FIN,FOUT,PUNCH 500 FORMAT(5X,'NAME OF INPUT FILE: ',A32,/ 1 5X,'NAME OF FIRST OUTPUT FILE: ',A32,/ 2 5X,'NAME OF SECOND OUTPUT FILE: ',A32,/) C WRITE(6,100) WRITE(6,101) WW = .0624 GT = 32.2 MAMAX=16400 C ...................................................................... C NAX = MAMAX + 5 NAA = NAX + 3*(MAMAX + 4) NS = NAA + 2*MAMAX NINV = NS + NAX/8 + 1 NTOT = NINV + NAX/8 + 1 WRITE(6,2000) MAMAX, NTOT C CALL SHAKIT(X(1), X(NAX), X(NAA), X(NS), X(NINV)) C WRITE(*,1500) WRITE(6,1750) STOP C **************************************************** 1000 FORMAT(I5, F10.0) 1500 FORMAT(//,3X,'I have enjoyed the brief time that we have spent ' 1 'together.'/) 1750 FORMAT(/,3X,'Analysis Complete: Program Terminated Normally') 2000 FORMAT( 45H MAX. NUMBER OF TERMS IN FOURIER TRANSFORM = I10/ 1 45H NECESSARY LENGTH OF BLANK COMMON X = I10) END C******************************************************************** SUBROUTINE EARTHQ(X,AX,S,INV) C*********************************************************************** C C THIS ROUTINE READS THE MOTION IN THE TIME DOMAIN, ADDS TRAILING C ZEROS, SCALES THE VALUES, FIND MAXIMUM VALUE AND VARIOUS PARAMETERS C AND TRANSFER THE MOTION INTO THE FREQUENCY DOMAIN. C C CODED BY PBS SEPT. 1970 C C X = INPUT MOTION C AX = TEMPORARY STORAGE OF X C TITLE = IDENTIFICATION FOR MOTION C DT = TIME STEP BETWEEN VALUES IN TIME DOMAIN C NV = NUMBER OF ACC. VALUES TO BE READ C MA = LENGTH OF MOTION INCLUDING TRAILING ZEROS C MMA = LENGTH OF SIGNIFICANT PART OF MOTION C XF = MULTIPLICATION FACTOR FOR ACCELERATION VALUES C DF = FREQUENCY STEPS IN FREQ. DOMAIN C C*********************************************************************** C CHARACTER*6 TITLE CHARACTER*30 FINPEQ CHARACTER*80 HEAD CHARACTER*12 FMAT COMPLEX X, AX DIMENSION XR(8),X(300),AX(3,270),S(70),INV(70) COMMON /EQ/ MFOLD,MA2,TITLE(5),DT, MA , MMA, DF, MX COMMON/FRCUT/ NCUT,NZERO COMMON /JISCK/ JIS,FINPEQ C PI2 = 6.283185307 READ(5,1001) NV, MA, DT, FINPEQ, FMAT READ(5,1004) XF,XMAX,FMAX,NHEAD, NPL IF (FMAX.LT. .001) FMAX = 100000. IF (FMAT .EQ. ' ') FMAT = '(8F9.6,I7)' MA2=2 2 IF(MA2.GE.MA) GO TO 3 MA2=MA2*2 GO TO 2 3 MA=MA2 C .................................................................... WRITE(6,2012) FINPEQ, NV, MA, NHEAD, NPL, DT, FMAT WRITE(7,2013) FINPEQ, NV, MA, DT WRITE (*,2026) FINPEQ,FMAT OPEN(8,FILE=FINPEQ,STATUS='OLD') WRITE (6,2021) DO 4 I=1,NHEAD READ(8,2022) HEAD 4 WRITE(6,2022) HEAD C WRITE (6,2023) MMA = NV + NV/10 IF (MMA.GT.MA) MMA=MA MA2 = MA + 2 MFOLD = MA2/2 MFOLD = MFOLD + 1 DF = 1./(MA*DT) FMA = FLOAT(MA) MX = (ALOG10(FMA)/ALOG10(2.))- 1 1 NMX = 2**(MX+1) IF (MA .LE. NMX) GO TO 11 MX = MX + 1 GO TO 1 11 NCARDS = (NV-1)/NPL + 1 JL=NPL*NCARDS-NV NV = NV + 1 N = 0 LC = 0 WRITE(6,2024) DO 31 I = 1,NCARDS LC = LC + 1 READ(8,FMAT) (XR(J), J=1,NPL) IF(I.NE.NCARDS) GO TO 6 IF(JL.EQ.0) GO TO 6 JL=NPL+1-JL DO 5 J=JL,NPL 5 XR(J)=0. C--------------------------------------------------------------- C ONLY PRINT OUT A FEW OF THE FIRST AND THE LAST LINES of Input Motion C--------------------------------------------------------------- 6 ICHECK = NCARDS - I IF (I .LE. 5 .OR. ICHECK .LT. 5) WRITE(6,2008) I,(XR(J), J=1,NPL) IF (I .EQ. 10) WRITE (6,2009) 2009 FORMAT(3X,'........ INPUT MOTION READ NOT ECHOED...........') C ENDIF C C FIND MAX. INPUT ACC. (XMAX) C 311 DO 31 J = 1,NPL,2 N = N + 1 X(N) = CMPLX(XR(J),XR(J+1)) 31 CONTINUE CLOSE (8) N = N + 1 DO 32 I = N,MFOLD 32 X(I) = 0. CALL XMX(X,MA,XM,NXMAX) IF (XMAX.LT. .000001) GO TO 300 XF = XMAX/XM 300 DO 30 I = 1,N 30 X(I) = X(I)*XF XMAX = XM*XF TMAX = FLOAT(NXMAX-1)*DT WRITE(6,2014) XM,TMAX,XF,XMAX C CALL RFFT(X,MX,INV,S,IFERR,1) C C REMOVE FREQUENCIES ABOVE FMAX AND FIND MAX. ACC. OF NEW MOTION C FREQ = 0. SXX = 0. SFX = 0. NCUT=0 DO 33 I = 1,MFOLD IF(FREQ.LE.FMAX) GO TO 34 NCUT=NCUT+1 X(I)=0.0 34 CONTINUE XA = CABS(X(I)) SXX= SXX + XA*XA SFX = SFX + FREQ*XA*XA AX(1,I) = X(I) FREQ = FREQ + DF 33 CONTINUE SFX = SFX/SXX NCUT=MFOLD-NCUT NZERO=NCUT+1 WRITE(6,2005) SFX IF (FMAX.GT.FREQ) RETURN CALL RFSN(X,MX,INV,S,IFERR,-2) CALL XMX(X,MA,XM,NXMAX) DO 72 I = 1,MFOLD 72 X(I) = AX(1,I) WRITE(6,2001) XM,FMAX C 1001 FORMAT(2I5, F10.3, A30, A12) 1002 FORMAT(8F9.5,I7) 1003 FORMAT(8F10.0) 1004 FORMAT(3F10.0,2I5) 2001 FORMAT(21H MAX ACCELERATION = F10.5, 22H FOR FREQUENCIES REMOV 19HED ABOVE F10.2, 7H C/SEC.) 2003 FORMAT(17H ACC. CARD NO. I4,16H OUT OF SEQUENCE ) 2005 FORMAT(25H MEAN SQUARE FREQUENCY = F10.2, 7H C/SEC. ) C2008 FORMAT(2X, I5, 5X, 8F15.6) 2008 FORMAT(1X,I5,1X, 8F9.6) 2012 FORMAT(/1X, ' FILE NAME FOR INPUT MOTION = ', A30,/ + 1X, ' NO. OF INPUT ACC. POINTS = ',I5,/ + 1X, ' NO. OF POINTS USED IN FFT = ',I5/ + 1X, ' NO. OF HEADING LINES = ',I5/ + 1X, ' NO. OF POINTS PER LINE = ',I5/ + 1X, ' TIME STEP FOR INPUT MOTION = ',F6.4/ + 1X, ' FORMAT FOR OF TIME HISTORY = ', A12, /) 2013 FORMAT(/1X, ' FILE NAME FOR INPUT MOTION = ', A30,/ + 1X, ' NO. OF INPUT ACC. POINTS = ',I5,/ + 1X, ' NO. OF POINTS USED IN FFT = ',I5/ + 1X, ' TIME STEP FOR INPUT MOTION = ',F6.4/) 2014 FORMAT(/23H MAXIMUM ACCELERATION = F9.5/ 1 23H AT TIME = F6.2, 4H SEC/ 1 44H THE VALUES WILL BE MULTIPLIED BY A FACTOR = F7.3/ 3 44H TO GIVE NEW MAXIMUM ACCELERATION = F9.5 ) 2021 FORMAT (/1X,'***** H E A D E R ') 2022 FORMAT (A80) 2023 FORMAT (1X,'*************************************************' 1 ,'*****************') 2024 FORMAT (' ** FIRST & LAST 5 LINES OF INPUT MOTION *****'/) 2025 FORMAT (' *****************************************************'/) 2026 FORMAT(/,1X, ' READING INPUT MOTION FROM ----> ',A30/ + 1X, ' FORMAT OF INPUT MOTION USED --> ',A12) RETURN END