[数值问题] 程序的变量值变化的很奇怪，求助

PROGRAM PFSAP!ANALYSIS PROGRAM FOR PLANE FRAME
INTEGER NO,NJ,N,NE,NM,NPJ,NPF,I,J,IE,IP,KK,I1
REAL K(200,200),KE(6,6),AKE(6,6),X(100),Y(100),AL(100),EAI(3,100),PJ(100),PF(2,100),R(6,6),P(100),FF(6),FE(6),ADE(6),DE(6),RT(6,6),AFE(6),F(3),D(100)
INTEGER JE(2,100),JN(3,100),JPJ(100),JPF(2,100),M(6),JEAI(100)
CHARACTER *80 TITLE
1 OPEN (6,FILE='PFSAP.IN')
OPEN (8,FILE='PFSAP.OUT')
READ(6,'(I3)')NO
IF(NO.EQ.0) STOP
WRITE(8,"('(NO.='I3,')'/)")NO
CALL PUTIN(NJ,N,NE,NM,NPJ,NPF,JN,X,Y,JE,JEAI,EAI,JPJ,PJ,JPF,PF)
DO I=1,N
   P(I)=0.0
   DO J=1,N
      K(I,J)=0
   END DO
END DO
DO IE=1,NE
   CALL MKE(KE,IE,JE,JEAI,EAI,X,Y,AL)
   CALL MR(R,IE,JE,X,Y)
   CALL MAKE(KE,R,AKE)
   CALL CALM(M,IE,JN,JE)
   CALL MK(K,AKE,M)
END DO
DO IP=1,NPF
   CALL MR(R,JPF(1,IP),JE,X,Y)
   CALL TRAN(R,RT)
   CALL PE(FE,IP,JPF,PF,AL)
   CALL MULV6(RT,FE,AFE)
   CALL CALM(M,JPF(1,IP),JN,JE)
   CALL MF(P,AFE,M)
END DO
DO I=1,NPJ
   P(JPJ(I))=P(JPJ(I))+PJ(I)
END DO
CALL SLOV(K,P,D,N)
WRITE(8,*) "/10X,5('* *'),'RESULTS OF CALCULATION',5('* *')"
WRITE(8,40)
40 FORMAT(/5X,'NO.N',4X,'X-DISPLACEMENT',2X,'Y-DISPLEMENT',3X,'ANG.ROT.(RAD)')
DO KK=1,NJ
   DO I=1,3
      F(I)=0.0
      I1=JN(I,KK)
      IF(I1.GT.0) F(I)=D(I1)
   END DO
WRITE(8,70)KK,F(1),F(2),F(3)
70 FORMAT(I8,2X,3G16.5)
END DO
WRITE(8,80)
80 FORMAT(/'NO.E',5X,'N(1)',8X,'Q(1)',8X,'M(1)',8X,'N(2)',8X,'Q(2)',8X,'M(2)')
DO IE=1,NE
   CALL MADE(IE,JN,JE,D,ADE)
   CALL MKE(KE,IE,JE,JEAI,EAI,X,Y,AL)
   CALL MR(R,IE,JE,X,Y)
   CALL MULV6(R,ADE,DE)
   CALL MULV6(KE,DE,FF)
     DO IP=1,NPF
        IF(JPF(1,IP).EQ.IE)THEN
          CALL PE(FE,IP,JPF,PF,AL)
          DO I=1,6
             FF(I)=FF(I)-FE(I)
          END DO
        END IF
     END DO
    WRITE(8,110)IE,(FF(I),I=1,6)
110    FORMAT(I5,2X,6G12.5)
END DO
GOTO 1
END PROGRAM PFSAP
SUBROUTINE PUTIN(NJ,N,NE,NM,NPJ,NPF,JN,X,Y,JE,JEAI,EAI,JPJ,PJ,JPF,PF)
REAL X(100),Y(100),EAI(3,100),PJ(100),PF(2,100)
INTEGER JE(2,100),JN(3,100),JPJ(100),JPF(2,100),JEAI(100)
CHARACTER*80 TITLE
READ(6,*)TITLE
WRITE(8,*)TITLE
READ(6,*)NJ,N,NE,NM,NPJ,NPF
WRITE(8,3)NJ,N,NE,NM,NPJ,NPF
3 FORMAT('NJ='I2,5X,'N='I2,5X,'NE='I2,5X,'NM='I2,5X,'NPJ='I2,5X,'NPF='I2)
WRITE(8,5)
5 FORMAT(/3X,'NO. (1)  (2)  (3)',10X,'X',8X,'Y')
READ(6,10)((JN(J,I),J=1,3),X(I),Y(I),I=1,NJ)
10 FORMAT(2(3I5,2G16.4))
WRITE(8,11)(I,(JN(J,I),J=1,3),X(I),Y(I),I=1,NJ)
11 FORMAT(2X,I2,3I5,4X,2F10.3)
WRITE(8,30)
30 FORMAT(/10X,'ELEMENT NO. NODE-1  NODE-2  MATERIALS')
READ(6,40)(JE(1,I),JE(2,I),JEAI(I),I=1,NE)
40 FORMAT(5(3I5))
DO I=1,NE
   WRITE(8,'(14X,I2,3(7X,I3))')I,JE(1,I),JE(2,I),JEAI(I)
END DO
READ(6,*)((EAI(I,J),I=1,3),J=1,NM)
WRITE(8,60)(J,(EAI(I,J),I=1,3),J=1,NM)
60 FORMAT(/3X,'NO.MAT',3X,'ELASTIK MOUDULUS',5X,'AREA',5X,'MOMENT OF INERTIA'/(I6,3G16.6))
IF(NPJ.NE.0)THEN
  WRITE(8,15) 
15 FORMAT(/20X,'NODEL LOADS')
  WRITE(8,12) 
12FORMAT(/16X,'NO.DISP.     VALUE')
  READ(6,70)(JPJ(I),PJ(I),I=1,NPJ)
70 FORMAT(5(I5,G16.4))
  DO I=1,NPJ
   WRITE(8,'(14X,I7,F16.3)')JPJ(I),PJ(I)
  END DO
ELSE
END IF
IF(NPF.NE.0) THEN
  WRITE(8,13) 
13 FORMAT(/20X,'NON-NODEL LOADS')
  WRITE(8,14) 
14 FORMAT(11X,'NO.E NO.LOAD.MODEL',8X,'A',9X,'C')
  READ(6,100)(JPF(1,I),JPF(2,I),PF(1,I),PF(2,I),I=1,NPF)
100 FORMAT(2(2I5,2G16.4))
  DO I=1,NPF
     WRITE(8,120)(JPF(J,I),J=1,2),PF(1,I),PF(2,I)
120 FORMAT(6X,2I8,10X,2F10.3)
  END DO
ELSE 
END IF
RETURN
END
SUBROUTINE MKE(KE,IE,JE,JEAI,EAI,X,Y,AL)
INTEGER I,J,II,JJ,MT
REAL L,A1,A2,A3,A4
REAL KE(6,6),X(100),Y(100),EAI(3,100),AL(100)
INTEGER JE(2,100),JEAI(100)
II=JE(1,IE)
JJ=JE(2,IE)
MT=JEAI(IE)
L=SQRT((X(JJ)-X(II))**2+(Y(JJ)-Y(II))**2)
AL(IE)=L
A1=EAI(1,MT)*EAI(2,MT)/L
A2=EAI(1,MT)*EAI(3,MT)/L**3
A3=EAI(I,MT)*EAI(3,MT)/L**2
A4=EAI(1,MT)*EAI(3,MT)/L
KE(1,1)=A1
KE(1,4)=-A1
KE(2,2)=12*A2
KE(2,3)=6*A3
KE(2,5)=-12*A2
KE(2,6)=6*A3
KE(3,3)=4*A4
KE(3,5)=-6*A3
KE(3,6)=2*A4
KE(4,4)=A1
KE(5,5)=12*A2
KE(5,6)=-6*A3
KE(6,6)=4*A4
DO I=1,6
   DO J=I,6
      KE(J,I)=KE(I,J)
   END DO
END DO
RETURN
END
SUBROUTINE MR(R,IE,JE,X,Y)
INTEGER I,J
INTEGER JE(2,100)
REAL L,CX,CY
REAL R(6,6),X(100),Y(100)
I=JE(1,IE)
J=JE(2,IE)
L=SQRT((X(J)-X(I))**2+(Y(J)-Y(I))**2)
CX=(X(J)-X(I))/L
CY=(Y(J)-Y(I))/L
DO J=1,6
   DO I=1,6
      R(I,J)=0.0
   END DO
END DO
DO I=1,4,3
   R(I,I)=CX
   R(I,I+1)=CY
   R(I+1,I)=-CY
   R(I+1,I+1)=CX
   R(I+2,I+2)=1.0
END DO
RETURN 
END
SUBROUTINE MAKE(KE,R,AKE)
REAL KE(6,6),R(6,6),RT(6,6),TMP(6,6),AKE(6,6)
CALL TRAN(R,RT)
CALL MULV(RT,KE,TMP)
CALL MULV(TMP,R,AKE)
RETURN
END
SUBROUTINE CALM(M,IE,JN,JE)
INTEGER I,IE
INTEGER M(6),JN(3,100),JE(2,100)
DO I=1,3
   M(I)=JN(I,JE(1,IE))
   M(I+3)=JN(I,JE(2,IE))
END DO
RETURN
END
SUBROUTINE MK(K,AKE,M)
INTEGER I,J
REAL K(200,200),AKE(6,6)
INTEGER M(6)
DO I=1,6
   DO J=1,6
      IF(M(I).NE.0.AND.M(J).NE.0) K(M(I),M(J))=K(M(I),M(J))+AKE(I,J)
   END DO
END DO
RETURN
END
SUBROUTINE PE(FE,IP,JPF,PF,AL)
INTEGER I,IND
INTEGER JPF(2,100)
REAL L,A,C
REAL FE(6),PF(2,100),AL(100)
A=PF(1,IP)
C=PF(2,IP)
L=AL(JPF(1,IP))
IND=JPF(2,IP)
DO I=1,6
   FE(I)=0.0
END DO
SELECT CASE(IND)
  CASE(1)
   FE(2)=(7.0*A/20.0+3.0*C/20.0)*L
   FE(3)=(A/20.0+C/30.0)*L**2
   FE(5)=(3.0*A/20.0+7.0*C/20.0)*L
   FE(6)=-(A/30.0+C/20.0)*L**2
   RETURN
  CASE(2)
   FE(5)=A*C**3*(2*L-C)/(2*L**3)
   FE(6)=-A*C**3*(4.0*L-3.0*C)/(12.0*L**2)
   FE(2)=A*C-FE(5)
   FE(3)=A*C**2*(6.0*L**2-8.0*C*L+3.0*C**2)/(12.0*L**2)
   RETURN
  CASE(3)
   FE(2)=A*(L-C)**2*(L+2*C)/L**3
   FE(3)=A*C*(L-C)**2/L**2
   FE(5)=A-FE(2)
   FE(6)=-A*C**2*(L-C)/L**2
   RETURN
  CASE(4)
   FE(2)=-6.0*A*C*(L-C)/L**3
   FE(3)=A*(L-C)*(L-3.0*C)/L**2
   FE(5)=-FE(2)
   FE(6)=A*C*(3.0*C-2.0*L)/L**2
   RETURN
  CASE(5)
   FE(1)=A*(1-C/L)
   FE(4)=A*C/L
  CASE(6)
   FE(1)=C*L/2
   FE(4)=FE(1)
 END SELECT
END
SUBROUTINE MULV6(A,B,C)
INTEGER I,J
REAL C(6),A(6,6),B(6)
DO I=1,6
   C(I)=0.0
   DO J=1,6
      C(I)=C(I)+A(I,J)*B(J)
   END DO
END DO
RETURN
END
SUBROUTINE MF(P,AFE,M)
INTEGER I
INTEGER M(6)
REAL P(100),AFE(6)
DO I=1,6
   IF(M(I).NE.0)P(M(I))=AFE(I)+P(M(I))
END DO
RETURN
END
SUBROUTINE SLOV(K,P,D,N)
INTEGER I,J,KK,N
REAL K(200,200),P(100),D(100)
DO I=1,100
   D(I)=P(I)
END DO
DO KK=1,N-1
   DO I=KK+1,N
      C=-K(KK,I)/K(KK,KK)
      DO J=I,N
         K(I,J)=K(I,J)+C*K(KK,J)
      END DO
      D(I)=D(I)+C*D(KK)
   END DO
END DO
D(N)=D(N)/K(N,N)
DO I=N-1,1,-1
   DO J=I+1,N
      D(I)=D(I)-K(I,J)*D(J)
   END DO
   D(I)=D(I)/K(I,I)
END DO
RETURN
END
SUBROUTINE MADE(IE,JN,JE,D,ADE)
INTEGER I
REAL ADE(6),D(100)
INTEGER IE,JN(3,100),JE(2,100)
DO I=1,6
   ADE(I)=0
END DO
DO I=1,3
   IF(JN(I,JE(1,IE)).NE.0) ADE(I)=D(JN(I,JE(1,IE)))
   IF(JN(I,JE(2,IE)).NE.0) ADE(I+3)=D(JN(I,JE(2,IE)))
END DO
RETURN
END
SUBROUTINE TRAN(R,RT)
INTEGER I,J
REAL R(6,6),RT(6,6)
DO I=1,6
   DO J=1,6
      RT(I,J)=R(J,I)
   END DO
END DO
RETURN
END
SUBROUTINE MULV(A,B,C)
INTEGER I,J,K
REAL A(6,6),B(6,6),C(6,6)
DO I=1,6
   DO J=1,6
      C(I,J)=0.0
      DO K=1,6
         C(I,J)=C(I,J)+A(I,K)*B(K,J) 
      END DO
   END DO
END DO
RETURN
END