eof运行不出来？？什么原因

云一林

我参考了以下帖子的程序：http://bbs.06climate.com/forum.php?mod=viewthread&tid=13501&extra=&highlight=eof&page=1，参考的是带north检验的eof程序
以下为eof修改部分：
      PROGRAM EOF
!********************************************************************************
!C     THIS PROGRAM USES EOF FOR ANALYSING TIME SERIES
!C     OF METEOROLOGICAL FIELD
!C     M:LENGTH OF TIME SERIES
!C     N:NUMBER OF GRID-POINTS
!C     KS=-1:SELF; KS=0:DEPATURE; KS=1:STANDERDLIZED DEPATURE
!C     KV:NUMBER OF EIGENVALUES WILL BE OUTPUT
!C     KVT:NUMBER OF EIGENVECTORS AND TIME SERIES WILL BE OUTPUT
!C     MNH=MIN(M,N)
!C     Evf=EIGENVACTORS, tcF=TIME COEFFICIENTS FOR EGVT.
!C     ER(KV,1)=LAMDA,LAMDA EIGENVALUE
!C     ER(KV,2)=ACCUMULATE LAMDA
!C     ER(KV,3)=THE SUM OF COMPONENTS VECTORS PROJECTED ONTO EIGENVACTOR.
!C     ER(KV,4)=ACCUMULATE ER(KV,3)
!********************************************************************************
!C     Both of grid data and station data are used in this program.
!C     M is the length of the time series, N is the number of grid, N=x*y
!C     MNH=MIN(M,N)
!C     KV(space pattern),KVT(time pattern) Generally, KV=KVT
!C     ff is undef, this parameter is based on origion data
!C     INPUT DATA---- modify  x and y  
!C     tt(x,y,M)                                          
!C     KV>=3
!C     modify one input path and three output paths  
!C     特别声明该eof程序在输出方差贡献的时候对方差贡献进行了north检验，其中样本时间长度需要手动更改
!*******************************************************************************
      PARAMETER(N=76,M=57)此程序中N为站点76，时间为57年
      PARAMETER(MNH=57,KS=1,KV=4,KVT=4)因为57比76小，所以MNH取57，数据已经手动距平，所以取KS=1做标准化eof，只想要4个模态，所以KV,KVT取4
      PARAMETER(ff=-9.99E+08)
      DIMENSION F(N,M),A(MNH,MNH),S(MNH,MNH),ER(mnh,4),DF(N),V(MNH),AVF(N),evf(N,KVT),tCF(M,KVT),nf(N),north(mnh)
!****  INPUT DATA ***********************************************
        open(11,file='1.txt',form='binary')读入数据
        do i=1,n
    do it=1,m
        read(11)f(i,it)
    enddo
    enddo
    close(11)
!*************************************************************
      CALL Test1(n,m,ff,f,nf)
           write(*,*)'ok2'
      CALL TRANSF(N,M,F,nf,AVF,DF,KS)
           write(*,*)'ok3'
      CALL FORMA(N,M,MNH,F,A)
           write(*,*)'ok4'
      CALL JCB(MNH,A,S,0.00001)
           write(*,*)'ok5'
      CALL ARRANG(KV,MNH,A,ER,S)
           write(*,*)'ok6'
      CALL TCOEFF(KVT,KV,N,M,MNH,S,F,V,evf,tcf,ER)
           write(*,*)'ok7'
      call test3(N,ff,nf,evf,kvt)
           write(*,*)'ok8'

!*************************************************************
!****   output the space field   ******
     open(21,file='Veof-f.txt',form='binary')将数据输出为txt格式
       do 668 kk=1,kvt
668    write(21)(evf(i,kk),i=1,n)
      close(21)

!*************************************************************
!****   output the time series   ******
      open(21,file='Teof.txt',form='binary')将数据输出为txt格式
      do 345 it=1,m
345        write(21)(tcf(it,iik),iik=1,kvt)
      close(21)

!*************************************************************
!****   output the VARiance  ******
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 在这里加入North检验，输出为文本格式文档 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      ee=sqrt(2/57.0)                  !此处需要手动修改样本（时间）长度，即修改395.0（务必为实型）样本时间长度为57，修改为57.0
      open(21,file='Var-east.txt')输出为txt格式
        do iiii=1,kv
                  if(er(iiii,3)-er(iiii+1,3)>=er(iiii,3)*ee) then
          write(21,106)(er(iiii,j),j=1,4),'pass'
          else
          write(21,106)(er(iiii,j),j=1,4),'notpass'
          end if
                 enddo
      close(21)
      STOP
106 format(4f18.4,a10)
      END


!*************************************************************
!******************************      
      SUBROUTINE Test1(n1,m,ff,f,nf)
      DIMENSION F(N1,M)
      DIMENSION nF(N1)
      do i=1,n1
      nf(i)=0.0
      enddo

      do i=1,n1
      do k=1,m
      if(f(i,k).eq.ff)then
      f(i,k)=0.0
      nf(i)=nf(i)+1
      endif
      enddo
      enddo
      return
      end

!*************************************************************
!***  THIS SUBROUTINE PROVIDES INITIAL F BY KS and kv.
!****************************** ******************************
      SUBROUTINE TRANSF(n1,m,f,nf,avf,df,ks)
      DIMENSION F(N1,M),AVF(N1)
      DIMENSION DF(N1)
      DIMENSION nF(N1)
      if(ks.eq.-1)then
      goto 30
      endif
      do i=1,n1
      avf(i)=0.0
      enddo
      if(ks.eq.0)then
      goto 5
      endif
      do i=1,n1
      df(i)=0.0
      enddo
!************************
5     continue
      DO 141 I=1,N1
      if(nf(i).ne.0) goto 141
      do 12 j=1,m
  12  AVF(I)=AVF(I)+F(I,J)/m
      DO 14 J=1,M
14    F(I,J)=F(I,J)-AVF(I)
141  CONTINUE
      IF(KS.EQ.0) THEN
      RETURN
      ELSE
      DO 241 I=1,N1
      if(nf(i).ne.0) goto 241
      DO 22 J=1,M
  22  DF(I)=DF(I)+F(I,J)*F(I,J)
      DF(I)=SQRT(DF(I)/M)
      DO 24 J=1,M
  24  F(I,J)=F(I,J)/DF(I)
241  CONTINUE
      ENDIF
  30  CONTINUE

      RETURN

      END
!*************************************************************
!***  THIS SUBROUTINE FORMS A BY F.
!****************************** ******************************
      SUBROUTINE FORMA(N,M,MNH,F,A)
      DIMENSION F(N,M),A(MNH,MNH)
      IF(M-N) 40,50,50
  40  DO 44 I=1,MNH
      DO 44 J=I,MNH
      A(I,J)=0.0
      DO 42 IS=1,N
  42  A(I,J)=A(I,J)+F(IS,I)*F(IS,J)
      A(J,I)=A(I,J)
  44  CONTINUE
      RETURN
  50  DO 54 I=1,MNH
      DO 54 J=I,MNH
      A(I,J)=0.0
      DO 52 JS=1,M
  52  A(I,J)=A(I,J)+F(I,JS)*F(J,JS)
      A(J,I)=A(I,J)
  54  CONTINUE
      RETURN
      END

!*************************************************************
!****  THIS SUBROUTINE COMPUTS EIGENVALUES AND EIGENVECTORS OF A
!*************************************************************
      SUBROUTINE JCB(N,A,S,EPS)
      DIMENSION A(N,N),S(N,N)
      DO 30 I=1,N
      DO 30 J=1,I
      IF(I-J) 20,10,20
  10  S(I,J)=1.
      GO TO 30
  20  S(I,J)=0.
      S(J,I)=0.
  30  CONTINUE
      G=0.
      DO 40 I=2,N
      I1=I-1
      DO 40 J=1,I1
  40  G=G+2.*A(I,J)*A(I,J)
      S1=SQRT(G)
      S2=EPS/FLOAT(N)*S1
      S3=S1
      L=0
  50  S3=S3/FLOAT(N)
  60  DO 130 IQ=2,N
      IQ1=IQ-1
      DO 130 IP=1,IQ1
      IF(ABS(A(IP,IQ)).LT.S3) GOTO 130
      L=1
      V1=A(IP,IP)
      V2=A(IP,IQ)
      V3=A(IQ,IQ)
      U=0.5*(V1-V3)
      IF(U.EQ.0.0) G=1.
      IF(ABS(U).GE.1E-10) G=-SIGN(1.,U)*V2/SQRT(V2*V2+U*U)
      ST=G/SQRT(2.*(1.+SQRT(1.-G*G)))
      CT=SQRT(1.-ST*ST)
      DO 110 I=1,N
      G=A(I,IP)*CT-A(I,IQ)*ST
      A(I,IQ)=A(I,IP)*ST+A(I,IQ)*CT
      A(I,IP)=G
      G=S(I,IP)*CT-S(I,IQ)*ST
      S(I,IQ)=S(I,IP)*ST+S(I,IQ)*CT
  110 S(I,IP)=G
      DO 120 I=1,N
      A(IP,I)=A(I,IP)
  120 A(IQ,I)=A(I,IQ)
      G=2.*V2*ST*CT
      A(IP,IP)=V1*CT*CT+V3*ST*ST-G
      A(IQ,IQ)=V1*ST*ST+V3*CT*CT+G
      A(IP,IQ)=(V1-V3)*ST*CT+V2*(CT*CT-ST*ST)
      A(IQ,IP)=A(IP,IQ)
  130 CONTINUE
      IF(L-1) 150,140,150
  140 L=0
      GO TO 60
  150 IF(S3.GT.S2) GOTO 50
      RETURN
      END

!*************************************************************************
!*****   THIS SUBROUTINE PROVIDES A SERIES OF EIGENVALUES FROM MAX TO MIN.
!*************************************************************************
      SUBROUTINE ARRANG(KV,MNH,A,ER,S)
      DIMENSION A(MNH,MNH),ER(mnh,4),S(MNH,MNH)
      TR=0.0
      DO 200 I=1,MNH
      TR=TR+A(I,I)
  200 ER(I,1)=A(I,I)
      MNH1=MNH-1
      DO 210 K1=MNH1,1,-1
      DO 210 K2=K1,MNH1
      IF(ER(K2,1).LT.ER(K2+1,1)) THEN
      C=ER(K2+1,1)
      ER(K2+1,1)=ER(K2,1)
      ER(K2,1)=C
      DO 205 I=1,MNH
      C=S(I,K2+1)
      S(I,K2+1)=S(I,K2)
      S(I,K2)=C
  205 CONTINUE
      ENDIF
  210 CONTINUE
      ER(1,2)=ER(1,1)
      DO 220 I=2,KV
      ER(I,2)=ER(I-1,2)+ER(I,1)
  220 CONTINUE
      DO 230 I=1,KV
      ER(I,3)=ER(I,1)/TR
      ER(I,4)=ER(I,2)/TR
  230 CONTINUE
      WRITE(6,250) TR
  250 FORMAT(/5X,'TOTAL SQUARE ERROR=',F20.5)
      RETURN
      END

!*************************************************************************
!***  THIS SUBROUTINE PROVIDES STANDARD EIGENVECTORS (M.GE.N,SAVED IN S;
!***  M.LT.N,SAVED IN F) AND ITS TIME COEFFICENTS SERIES (M.GE.N,SAVED
!***  IN F; M.LT.N,SAVED IN S)
!*************************************************************************
      SUBROUTINE TCOEFF(KVT,KV,N,M,MNH,S,F,V,evf,tcf,ER)
      DIMENSION S(MNH,MNH),F(N,M),V(MNH),ER(mnh,4),evf(n,kvt),tcf(m,kvt)
      DO 360 J=1,KVT
      C=0.
      DO 350 I=1,MNH
  350 C=C+S(I,J)*S(I,J)
      C=SQRT(C)
      DO 160 I=1,MNH
      s(I,J)=S(I,J)/C
  160 evf(I,J)=S(I,J)/C
  360 CONTINUE
!****************************************************************
!********cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      IF(N.LE.M) THEN
      DO 390 J=1,M
      DO 370 I=1,N
      V(I)=F(I,J)
      F(I,J)=0.
  370 CONTINUE
      do 371 is=1,kvt
         tcf(j,is)=0.
371  continue
      DO 380 IS=1,KVT
      DO 380 I=1,N
      f(is,j)=F(IS,J)+V(I)*S(I,IS)
  380 tcf(j,is)=tcf(J,is)+V(I)*S(I,IS)
  390 CONTINUE

!!****************************************************************
      ELSE
      DO 410 I=1,N
      DO 400 J=1,M
      V(J)=F(I,J)
      F(I,J)=0.
  400 CONTINUE
      DO 410 JS=1,KVT
      DO 410 J=1,M
      f(I,JS)=F(I,JS)+V(J)*S(J,JS)
  410 CONTINUE
      DO 430 JS=1,KVT
      DO 420 J=1,M
      tcf(J,JS)=S(J,JS)*SQRT(ER(JS,1))
  420 CONTINUE
      DO 430 I=1,N
      evf(I,JS)=F(I,JS)/SQRT(ER(JS,1))
  430 CONTINUE
      t=0.0
      do 3650 i=1,m
3650   t=t+tcf(i,1)*tcf(i,2)
      t=0.0
      do 3651 i=1,n
3651   t=t+evf(i,1)*evf(i,2)
      ENDIF
      RETURN
      END

!*************************************************************
!***   this subroutine sent undefine value ff to evf
!*************************************************************
      SUBROUTINE test3(N1,ff,nf,evf,kvt)
      dimension nf(n1),evf(n1,kvt)
      do i=1,n1
      if(nf(i).ne.0)then
      do k=1,kvt
      evf(i,k)=ff
      enddo
      endif
      enddo
      return
      end
!*********************************************************
      SUBROUTINE OUTER(KV,ER,mnh)
!**     THIS SUBROUTINE PRINTS ARRAY ER
!**     ER(KV,1) FOR  SEQUENCE OF EIGENVALUE FROM BIG TO SMALL
!**     ER(KV,2) FOR  EIGENVALUE FROM BIG TO SMALL
!**     ER(KV,3) FOR  SMALL LO=(LAMDA/TOTAL VARIANCE)
!**     ER(KV,4) FOR  BIG LO=SUM OF SMALL LO)
      DIMENSION ER(mnh,4)
      WRITE(16,510)
  510 FORMAT(/30X,'EIGENVALUE AND ANALYSIS ERROR',/)
      WRITE(16,520)
  520 FORMAT(10X,1HH,8X,5HLAMDA,10X,6HSLAMDA,11X,2HPH,12X,3HSPH)
      WRITE(16,530) (IS,(ER(IS,J),J=1,4),IS=1,KV)
  530 FORMAT(1X,I10,4F15.5)
      WRITE(16,540)
  540 FORMAT(//)
      RETURN
      END

!*****************************************************************
      SUBROUTINE OUTVT(KVT,N,M,MNH,S,F,EGVT,ECOF)
!**     THIS SUBROUTINE PRINTS STANDARD EIGENVECTORS AND ITS TIME-COEFFICENT SERIES
!*****************************************************************
      DIMENSION F(N,M),S(MNH,MNH),EGVT(N,KVT),ECOF(M,KVT)
      WRITE(16,560)
  560 FORMAT(30X,'STANDARD EIGENVECTORS',/)
      WRITE(16,570) (IS,IS=1,KVT)
  570 FORMAT(3X,80I7)
      DO 550 I=1,N
      IF(M.GE.N) THEN
      WRITE(16,580) I,(S(I,JS),JS=1,KVT)
  580 FORMAT(1X,I3,80F7.3,/)
      DO 11 JS=1,KVT
      EGVT(I,JS)=S(I,JS)
   11 CONTINUE
      ELSE
      WRITE(16,590) I,(F(I,JS),JS=1,KVT)
  590 FORMAT(1X,I3,80F7.3)
      DO 12 JS=1,KVT
      EGVT(I,JS)=F(I,JS)
   12 CONTINUE
      ENDIF
  550 CONTINUE
      WRITE(16,720)
  720 FORMAT(//)
      WRITE(16,610)
  610 FORMAT(30X,'TIME-COEFFICENT SERIES OF S. E.'/)
      WRITE(16,620) (IS,IS=1,KVT)
  620 FORMAT(3X,80I7)
      DO 600 J=1,M
      IF(M.GE.N) THEN
      WRITE(16,630) J,(F(IS,J),IS=1,KVT)
  630 FORMAT(1X,I3,80F7.1)
      DO 13 IS=1,KVT
      ECOF(J,IS)=F(IS,J)
   13 CONTINUE
      ELSE
      WRITE(16,640) J,(S(J,IS),IS=1,KVT)
  640 FORMAT(1X,I3,80F7.1)
      DO 14 IS=1,KVT
      ECOF(J,IS)=S(J,IS)
   14 CONTINUE
      ENDIF
  600 CONTINUE
      RETURN
      END
输入数据txt格式：

62.51754386        -20.88245614        -42.18245614        36.11754386        -74.88245614        96.41754386        -12.08245614        37.01754386        -95.28245614        -93.48245614        -57.58245614        10.91754386        28.81754386        -132.3824561        14.31754386        114.1175439        -101.2824561        23.21754386        77.51754386        33.61754386        63.61754386        -45.48245614        -1.78245614        204.0175439        40.71754386        3.91754386        -40.68245614        -79.98245614        63.31754386        -1.28245614        -88.18245614        12.71754386        80.01754386        -36.08245614        25.21754386        33.01754386        -73.88245614        127.8175439        -24.48245614        -26.28245614        -89.28245614        -65.78245614        117.2175439        48.11754386        -47.88245614        -64.08245614        25.31754386        -65.48245614        -83.28245614        226.9175439        18.31754386        106.8175439        56.81754386        -52.68245614        -89.38245614        -106.7824561        -75.68245614
-99.09649123        71.60350877        109.7035088        29.70350877        154.3035088        44.40350877        -113.5964912        -8.896491228        9.103508772        84.70350877        -110.3964912        -167.7964912        -124.0964912        85.20350877        -11.09649123        20.10350877        -12.79649123        -135.0964912        56.60350877        70.00350877        24.40350877        51.20350877        4.003508772        -20.39649123        43.40350877        -96.79649123        5.503508772        0.003508772        -11.89649123        7.603508772        46.20350877        -73.89649123        106.3035088        -118.8964912        130.2035088        -60.39649123        -83.49649123        83.50350877        78.80350877        8.803508772        -46.69649123        -23.69649123        130.4035088        -21.79649123        137.5035088        -134.3964912        -133.4964912        8.903508772        19.20350877        -67.09649123        -38.89649123        97.30350877        -24.29649123        93.00350877        -42.89649123        -75.79649123        46.00350877
-4.759649123        -22.95964912        -8.059649123        -53.85964912        62.04035088        -14.65964912等等
就是76行57列
fortran运行结果：已上传图片，请问这样的运行是什么原因？fortran并没有报错，但是就是运行到ok4就不往下了，请问是我数据的放置有问题吗？