[编程申请]：一杆系结构计算程序_FOR

拉登

C       ********************************************
C       * SEISMIC ACTION OF MULTISTORY SHEAR MODEL *
C       *   FRAME BY RESPONSE SPECTRUM TECHNIQUE   *
C       ********************************************
C                       1992.11.2  

        PROGRAM SEIS
        DIMENSION M(50),K(50),A(50,50),S(50,50),T(50),
     &  AF(50),R(50),OM(50),GM(50),F(50,50)
        INTEGER XD,CD
        REAL M,K
        OPEN(1,FILE='R.DAT')
        OPEN(2,FILE='W.DAT',STATUS='NEW')
        READ(1,*) N,XD,JY,LD,CD,NX
        READ(1,*) (M(I),I=1,N)
        READ(1,*) (K(I),I=1,N)
        WRITE(2,5) N,XD,JY,LD,CD,NX
5       FORMAT(5X,'N=',I2,3X,'XD=',I2,3X,'JY=',I2,
     &  3X,'LD=',I2,3X,'CD=',I2,3X,'NX=',I2)
        WRITE(2,10) (M(I),I=1,N)
10      FORMAT(12X,'LUMPED-MASS OF STORY'/(2X,5F12.3))
        WRITE(2,15) (K(I),I=1,N)
15      FORMAT(12X,'SHEAR STIFFNESS OF STORY'/2X,5F12.3))
        CALL FMA(N,M,K,A,R)
        CALL JACOBI(N,A,S,1.0E-6)
        CALL TVM(N,A,S,T,R,OM)
        CALL ALFA(XD,JY,LD,CD,NX,T,AF,N)
        CALL FIJ(N,NX,M,S,AF,GM,F)
        CLOSE(1)
        CLOSE(2)
        END

C       FORM MATRIX A
        SUBROUTINE FMA(N,M,K,A,R)
        DIMENSION R(N),M(N),K(N),A(N,N)
        REAL M,K
        DO 20 I=1,N
        R(I)=SQRT(M(I))
        DO 10 J=1,N
        A(I,J)=0.0
10      CONTINUE
20      CONTINUE
        DO 30 I=1,N-1
        A(I,I)=(K(I)+K(I+1))/M(I)
        A(I,I+1)=-K(I+1)/R(I)/R(I+1)
        A(I+1,I)=A(I,I+1)
30      CONTINUE
        A(N,N)=K(N)/M(N)
        END

C       JACOBI METHOD FOR FINDING EIGENVALUES AND
C       EIGENVECTORS OF SYMMERICAL REAL MATRIX
        SUBROUTINE JACOBI(N,A,S,EPS)
        DIMENSION A(N,N),S(N,N)
        INTEGER P,Q
        DO 10 I=1,N
        DO 10 J=1,N
        S(I,J)=0.0
10      S(I,I)=1.0
        G=0.0
        DO 20 I=2,N
        DO 20 J=1,I-1
        G=G+2.0*A(I,J)*A(I,J)
20      CONTINUE
        T1=SQRT(G)
        T2=EPS*T1/N
        T3=T1
        L=0
30      T3=T3/N
40      DO 80 Q=2,N
        DO 80 P=1,Q-1
        IF(ABS(A(P,Q)).GE.T3) THEN
        L=1
        V1=A(P,P)
        V2=A(P,Q)
        V3=A(Q,Q)
        U=0.5*(V1-V3)
        IF(U.EQ.0.0) G=1.0
        IF(ABS(U).GE.1.0E-10) G=-SIGN(1.0,U)*V2/
     &  SQRT(V2*V2+U*U)
        ST=G/SQRT(2.0*(1.0+SQRT(1.0-G*G)))
        CT=SQRT(1.0-ST*ST)
        DO 90 I=1,N
        G=A(I,P)*CT-A(I,Q)*ST
        A(I,Q)=A(I,P)*ST+A(I,Q)*CT
        A(I,P)=G
        G=S(I,P)*CT-S(I,Q)*ST
        S(I,Q)=S(I,P)*ST+S(I,Q)*CT
        S(I,P)=G
90      CONTINUE
        DO 100 I=1,N
        A(P,I)=A(I,P)
        A(Q,I)=A(I,Q)
100     CONTINUE
        A(P,P)=V1*CT*CT+V3*ST*ST-2.0*V2*ST*CT
        A(Q,Q)=V1*ST*ST+V3*CT*CT+2.0*V2*ST*CT
        A(P,Q)=0.0
        A(Q,P)=0.0
        END IF
80      CONTINUE
        IF(L.EQ.1) THEN
        L=0
        GO TO 40
        ELSE IF(T3.GT.T2) THEN
        GO TO 30
        END IF
        END

C       CALCULAT AND OUTPUT PERIODS AND VIBRATION MODES
        SUBROUTINE TVM(N,A,S,T,R,OM)
        DIMENSION A(N,N),S(N,N),T(N),R(N),OM(N)
        DO 10 I=1,N
        OM(I)=SQRT(A(I,I))
        DO 10 J=1,N
        S(I,J)=S(I,J)/R(I)
10      CONTINUE
        DO 40 I=1,N-1
        DO 30 J=I+1,N
        IF(OM(I).GT.OM(J)) THEN
        G=OM(I)
        OM(I)=OM(J)
        OM(J)=G
        DO 20 L=1,N
        G=S(L,I)
        S(L,I)=S(L,J)
        S(L,J)=G
20      CONTINUE
        END IF
30      CONTINUE
40      CONTINUE
        DO 45 J=1,N
        C=S(1,J)
        DO 45 I=1,N
        S(I,J)=S(I,J)/C
45      CONTINUE
        DO 50 J=1,N
        T(J)=2.0*3.14159/OM(J)
50      CONTINUE
        WRITE(2,60) (T(J),J=1,N)
60      FORMAT(/12X,'VIBRATION PERIODS'/(2X,5F12.4))
        WRITE(2,70)
70      FORMAT(/12X,'VIBRATION MODES')
        DO 80 J=1,N
80      WRITE(2,90) J,(S(I,J),I=1,N)
90      FORMAT(5X,'MODE NUMBER:',I3/(2X,5F12.4))
        END

C       CALCULAT SEISMIC EFFECT COEFFICIENT
        SUBROUTINE ALFA(XD,JY,LD,CD,NX,T,AF,N)
        INTEGER XD,CD
        DIMENSION T(N),AF(NX)
        IF(XD.EQ.1) THEN
        IF(LD-8) 10,20,30
10      AM=0.08
        GO TO 70
20      AM=0.16
        GO TO 70
30      AM=0.32
        GO TO 70
        ELSE
        IF(LD-8) 40,50,60
40      AM=0.50
        GO TO 70
50      AM=0.90
        GO TO 70
60      AM=1.40
        END IF
70      IF(JY.EQ.1) THEN
        GO TO(80,90,100,110) CD
80      TG=0.2
        GO TO 160
90      TG=0.3
        GO TO 160
100     TG=0.4
        GO TO 160
110     TG=0.65
        GO TO 160
        ELSE
        GO TO(120,130,140,150) CD
120     TG=0.25
        GO TO 160
130     TG=0.4
        GO TO 160
140     TG=0.55
        GO TO 160
150     TG=0.85
        END IF
160     DO 170 J=1,NX
        IF(T(J).LE.0.1) THEN
        AF(J)=0.45*AM+5.5*T(J)*AM
        ELSE IF(T(J).LE.TG) THEN
        AF(J)=AM
        ELSE IF(T(J).LE.3.0) THEN
        AF(J)=(TG/T(J))**0.9*AM
        IF(AF(J).LT.0.2*AM) AF(J)=0.2*AM
        END IF
170     CONTINUE
        END

C       CALCULAT AND OUTPUT SEISMIC ACTION
        SUBROUTINE FIJ(N,NX,M,S,AF,GM,F)
        DIMENSION M(N),S(N,N),AF(N),GM(NX),F(N,NX)
        REAL M
        DO 10 J=1,NX
        S1=0.0
        S2=0.0
        DO 5 I=1,N
        S1=S1+M(I)*S(I,J)
        S2=S2+M(I)*S(I,J)*S(I,J)
5       CONTINUE
        GM(J)=S1/S2
10      CONTINUE
        DO 20 J=1,NX
        DO 20 I=1,N
        F(I,J)=AF(J)*GM(J)*S(I,J)*M(I)*9.8
20      CONTINUE
        WRITE(2,30)
30      FORMAT(/12X,'SEISMIC ACTION')
        DO 40 J=1,NX
40      WRITE(2,50) J,(F(I,J),I=1,N)
50      FORMAT(5X,'MODE NUMBER:',I2/(2X,5F12.3))
        END

zhengqg

楼主是一名在校学生吧,这个程序在学校时有一种似曾相识的感觉

拉登

是一名在校学生

liqinlove

编写得不错，谁敢用