我用compaq visual fortran运行下列程序，报错，麻烦高手帮我调试一下，我都不知道错在哪了。。。
C  ***************************************************************
C  *              PROGRAM OF DISK STRESS ANANLYSIS               *
C  ***************************************************************
  DIMENSION R(21),H(21),SIGR(2,21),ALF(21),E(21),T(21),AN(2,2,19),
  *AM(2,2,20),XM(20),A(20),B(20),C(2,20),CL(2,20),AMN(2,2,20),
  *AI(2,2,20),BI(2,20),AIP(2,2,20)BIP(2,20),RE(22),HE(22),DD(2,21),
  *SIGE(2,21),SIGP(2,21),HH(21)
  DOUBLE PRECISION R,H,T,ALF,E,SIGR,AN,AM,XM,A,B,C,ALFA,EE,CL,F,
  *AMN,AI,BI,V,GAMA,NN,DELTA,XO,YO,HH
C  CHARACTER*(69)ZS
  OPEN(5,FILE='CIN.DATA')
  READ(5,*)N
C  READ(5,*)ZS
  READ(5,*)(R(I),HH(I),E(I),T(I),ALF(I),I=1,N),V,GAMA,NN
  CLOSE(5)
  WRITE(6,400)
  DO 10 I=1,N
10  WRITE(6,420) R(I),HH(I),E(I),T(I),ALF(I)
  WRITE(6,430) V,GAMA,NN
  WRITE(6,440)
  N1=N-1
  DO 20 I=1,N1
  XM(I)=R(1)/R(I+1)
  H(I)=(HH(I)+HH(I+1))/2
  A(I)=0.5*(1.0+XM(I)*XM(I))
  B(I)=0.5*(1.0-XM(I)*XM(I))
  C(1,I)=0.125*(2.0*(1.0+V)*XM(I)*XM(I)+(1.0-V)*XM(I)**4-3.0-V)
  C(2,I)=0.125*(2.0*(1.0+V)*XM(I)*XM(I)-(1.0-V)*XM(I)**4-1.0-3.0-V)
  CL(1,I)=C(1,I)*GAMA*(2.0*3.1415926*NN/60.0)**2*R(I+1)*R(I+1)
20  CL(2,I)=C(2,I)*GAMA*(2.0*3.1415926*NN/60.0)**2*R(I+1)*R(I+1)
  DO 25 I1=1,2
25  BI(I1,1)=CL(I1,1)
  AI(1,1,1)=A(1)
  AI(1,2,1)=B(1)
  AI(2,1,1)=B(1)
  AI(2,2,1)=A(1)
  I=1
100  AM(1,1,I+1)=A(I+I)
  AM(1,2,I+1)=B(I+I)
  AM(2,1,I+1)=B(I+I)
  AM(2,2,I+1)=A(I+I)
  AN(1,1,I)=H(I)/H(I+1)
  AN(1,2,I)=0.0
  AN(2,1,I)=V*(H(I)/H(I+1)-E(I+1)/E(I))
  AN(2,2,I)=E(I+1)/E(I)
  DD(2,I)=E(I+1)*(ALF(I)*T(I)-ALF(I+1)*T(T+1))
  DD(1,I)=0.0
  DO 30 I1=1,2
  BIP(I1,I)=0.0
  DO 30 I2=1,2
  AIP(I1,I2,I)=0.0
  DO 30 K=1,2
30  AIP (I1,I2,I)=AIP(I1,I2,I)+AN(I1,K,I)*AI(K,I2,I)
  DO 40 I1=1,2
  DO 35 I2=1,2
35  BIP(I1,I)=BIP(I1,I)+AN(I1,I2,I)*BI(I2,I)
40  BIP(I1,I)=BIP(I1,I)+DD(I1,I)
  DO 42 J1=1,2
  BI(J1,I+1)0.0
  DO 42 J2=1,2
  AI(J1,J2,I+1)=0.0
  DO 42 K=1,2
42  AI(J1,J2,I+1)=AI(J1,J2,I+I)+AM(J1,K,I+1)*AIP(K,J2,I)
  DO 46 J1=1,2
  DO 44 J2=1,2
44  BI(J1,I+1)=BI(J1,I+1)+AM(J1,J2,I+1)*BIP(J2,I)
46  BI(J1,I+1)=BI(J1,I+1)+CL(J1,I+1)
  I=I+1
  IF(I.LT.N1)GOTO 100
  IF(R(1).EQ.0.0) GOTO 300
  READ (5,*) XO,YO
C  
  SIGR(1,1)=XO
  SIGR(1,N)=YO
  SIGR(2,1)=(YO-BI(1,N1)-AI(1,1,N1)*XO)/AI(1,2,N1)
  GOTO 160
C 300  READ(5,*) XO
  SIGR(1,1)=(XO-BI(1,N1))/(AI(1,1,N1)+AI(1,2,N1))
  SIGR(2,1)=SIGR(1,1)
160  DO 60 I=1,N1
  DO 60 L1=1,2
  SIGR(L1,I+1)=BIP(L1,I)
60  SIGP(L1,I)=BI(L1,I)
  DO 80 I=1,N1
  DO 80 L1=1,2
  DO 70 L2=1,2
  IF (I.EQ.N1) GOTO 70
  SIGR(L1,I+1)=SIGR(L1,I+1)+AIP(L1,L2,I)*SIGR(L2,1)
70  SIGP(L1,I)=SIGP(L1,I)+AI(L1,L2,I)*SIGR(L2,1)
  SIGE(L1,I+1)=(SIGR(L1,I)+SIGP(L1,I))*0.5
  RE(I+1)=(R(I)+R(I+1))*0.5
80  HE(I+1)=H(I)
  DO 90 L2=1,2
  SIGE(L2,1)=SIGR(L2,1)
90  SIGE(L2,N+1)=SIGP(L2,N1)
  RE(1)=R(1)
  RE(N+1)=R(N)
  HE(1)=HH(1)
  HE(N+1)=HH(N)
  WRITE(6,650)
  WRITE(6,500)
  WRITE(6,650)
  IF(R(1).NE.0.0) GOTO 120
  WRITE(6,670)N,XO
  GOTO 121
120  WRITE(6,660) N,XO,YO
121  WRITE(6,850)
  WRITE(6,860)
  DO 600 I=1,N+1
  WRITE(6,800) RE(I),HE(I),SIGE(1,I),SIGE(2,I)
600  CONTINUE
400  FORMAT(7X,'R',13X,'HH',14X,'E',14X,'T',11X,'ALF',/)
420  FORMAT(1X,5(D12.7,4X))
430  FORMAT(1X,/,7X,'V',29X,'GAMA',26X,'NN',//,1X,3(D12.7,20X))
440  FORMAT(1X,/)
650  FORMAT(5X,55H****************************************************
  ****)
500  FORMAT(5X,55H*  STRESS DISTRIBUTION ALONG R DIRECTION)
670  FORMAT(14X,1H(,2X,2HN=,I2,13X,3HXO=,F10.2,3X,1H))
660  FORMAT(7X,1H(,2HN=,I2,13X,3HXO=,F10.2,4X,3HYO=,F10.2,2X,1H))
850  FORMAT(1X,6HR(cm),4X,6HH(cm.),7X,16HSIG--R(kg/cm*cm),8X,
  *16HSIG--Q(kg/cm*cm))
860  FORMAT(1X)
800  FORMAT(1X,F5.2,4X,F6.2,6X,D18.9,6X,D18.9)
  STOP
  END

  QQ

小白程序应先检查R是不是大于0