dat文件进行eof分析出现断点，跪求大神帮忙！

islandowner

如题，上代码
!此程序为EOF程序
! 运行时要改动前面的空间、时间格点以及文件路径，ks和kvt根据自己的需要进行改动
!程序中自动去除缺省值并将其写回生成数据(生成数据中缺省值为-9999.0)
!对程序中data_in到F的传递 进行调整后 此程序也可用于s-eof和mv-eof
!
PROGRAM EOF
    IMPLICIT NONE
    INTEGER,PARAMETER :: nt=612,nx=180,ny=120          ! you need change,NT为时间长度                                             ***************
    INTEGER,PARAMETER  ::M=nt,KS=0,KVT=1              !kvt为输出的模态数                                                         ***************
                                                      ! KS的设置： ks>0 计算前先将数据标准化 ，
                      ! ks=0时取距平，ks<0时不进行这一步处理  
    INTEGER :: i,j ,MNH,N ,K,IM , m1
    REAL, allocatable,dimension(:,:,:)::DATA_IN
    REAL, allocatable,dimension(:,:)::F,S,ER,A,S1,F1
    CHARACTER(LEN=20) :: NOW  , TRACK
    REAL :: land(nx,ny), D,AVE,PT(NX,NY,kvt)  ,ran1
     
     TRACK=''                                         !输出的目标文件夹，默认为程序所在文件夹
      call  time(now)
      print*, now   
!!1111111111读入数据并去掉缺省值11111111111111         
   ALLOCATE(DATA_IN(NX,NY,NT))
   OPEN(1,file='ssta.dat',access='direct',recl=nx*ny*nt)  !****修改路径                                                        ***************
      READ(1,rec=1) (((data_in(I,J,K),I=1,nx),J=1,ny),K=1,nt)   
      CLOSE(1)   
                                       !注意数据排列顺序

!************做纬度加权平均，中、高纬度使用，热带或小范围不必******（未验证）
      !do j=1,ny
      !z(j)=0.+(real(j)-1.)*2.5/180.*3.1415926575       !使用时需要改动格距和起始纬度
      !data_in(:,j,:)=data_in(:,j,:)*sqrt(cos(z(j)))
      !enddo   
   land=0.0
      N=NX*NY
      DO I=1,nx
      DO J=1,ny
   DO K=1,nt
        IF(abs(data_in(I,J,K))>1000)then             !判断缺省值 (注意条件)
       land(I,J)=-9999.0
    N=N-1
    EXIT
     ENDIF
   ENDDO
      ENDDO
      ENDDO
   ALLOCATE(F(1:N,1:M))
      im=0
      DO I=1,nx
      DO J=1,ny
        IF(land(I,J)/=-9999.0)then
      im=im+1
      F(IM,1:m)=data_in(I,J,1:m)      
        ENDIF
      ENDDO
      ENDDO     
      print*, '空间点数' , nx*ny, '非缺省值空间点数:',im,N
    DEALLOCATE(DATA_IN)
           MNH=min(N,M)
       ALLOCATE( A(MNH,MNH))
       ALLOCATE(S(MNH,MNH))
    ALLOCATE(ER(mnh,6))
      
!222222222222222222计算过程22222222222222222222222

      CALL TRANSF(N,M,F,KS)                           !根据KS的设置，-1时跳出，0时距平，1时标准化
   print*,"**"
      CALL FORMA(N,M,MNH,F,A)                         !求协方差矩阵A
   print*,"***"   
      CALL JCB(MNH,A,S,0.0000001)                       !雅可比过关法求特征值特征向量
    print*,"****"                               !最后这个EPS的值控制计算精度，越小精度越高  
      
      CALL ARRANG(MNH,A,ER,S)                         !按照特征值大小排序
   print*,"*****"
   DEALLOCATE( A)
      CALL TCOEFF(KVT,N,M,MNH,S,F,ER)                 !给出时间序列和标准化的空间场
   print*,"******"
      
  ALLOCATE(S1(MNH,MNH))
  ALLOCATE(F1(N,M))
        
!33333********数据输出**********333333333                  !输出数据为标准化后的时间序列及相应的空间场                     
!求时间序列的标准差，时间序列除以标准差，空间乘以该标准差
   IF (M>=N) THEN     
      DO K=1,KVT
        AVE=SUM(F(K,1:M))/REAL(M)
        F1(K,1:M)=F(K,1:M)-AVE
        D=   SQRT(SUM(F1(K,1:M)*F1(K,1:M))/REAL(M))
        F(K,1:M)= F(K,1:M)/D       !时间
        S(K,1:N)= S(K,1:N)*D       !空间
      ENDDO
      m1=0
      DO i=1,nx
      DO j=1,ny
       IF(land(i,j).eq.0.0)then
      m1=m1+1   
         PT(i,j,1:kvt)=S(m1,1:kvt)   
       ELSE
        PT(i,j,1:kvt)=-9999.00
       ENDIF
      ENDDO
      ENDDO
     OPEN (1,file=TRIM(track)//'pt.dat',access='direct',recl=NX*NY)
     DO k=1,KVT
       WRITE(1,rec=k)((PT(i,j,k),i=1,nx),j=1,ny)
     ENDDO
     CLOSE(1)
     OPEN (2,FILE=TRIM(track)//'PC.DAT',ACCESS='DIRECT',RECL=M)
     DO K=1,KVT
     WRITE(2,REC=K) ((F(K,1:M)))
     ENDDO
  CLOSE(2)
ELSE
   DO K=1,KVT
     AVE=SUM(S(1:M,K))/REAL(M)
     S1(1:M,K)=S(1:M,K)-AVE  
     D= SQRT( SUM(S1(1:M,K)*S1(1:M,K))/REAL(M))  
     S(1:M,K)=S(1:M,K)/D     !时间
     F(1:N,K)=F(1:N,K)*D     !空间
   ENDDO
    m1=0
    DO i=1,nx
    DO j=1,ny
      IF(land(i,j).eq.0.0)then
      m1=m1+1   
        PT(i,j,1:kvt)=F(m1,1:kvt)   
    ELSE
        PT(i,j,1:kvt)=-9999.00
       ENDIF
    ENDDO
    ENDDO
     OPEN (1,file=TRIM(track)//'PT.dat',access='direct',recl=NX*NY)   
     DO k=1,KVT
      WRITE(1,rec=k) ((PT(i,j,k),i=1,nx),j=1,ny)
     ENDDO
     CLOSE(1)
     OPEN (2,FILE=TRIM(track)//'PC.DAT',ACCESS='DIRECT',RECL=M)
     DO K=1,KVT
     WRITE(2,REC=K) ((S(1:M,K)))
     ENDDO
  CLOSE(2)
  OPEN(3,FILE=TRIM(TRACK)//'PC10.TXT')
  DO K=1,M
   WRITE(3,'(8F16.4)') S(K,1:KVT)
  ENDDO
     ENDIF
      call  time(now)
      print*, now ,'OK!'
    ENDPROGRAM

!######################################!
!                                      !
!     以下为计算过程调用的5个子程序    !
!                                      !
!######################################!
   
     
!11111111111111111!根据KS的设置进行初步处理
      SUBROUTINE TRANSF(N,M,F,KS)
   IMPLICIT NONE
!     THIS SUBROUTINE PROVIDES INITIAL F BY KS
      INTEGER ::KS, I,M,N
   REAL ::F(N,M),AVF(N),DF(N)   
   
      AVF=0.0
      DF=0.0
      IF(KS>0 .or. KS ==0) then     !根据KS的设置，-1时跳出，0时距平，1时标准化
        DO  I=1,N
         AVF(I)=SUM(F(I,1:M)/M)
            F(I,1:M)=F(I,1:M)-AVF(I)
        ENDDO
     IF(KS==0)   RETURN
        DO  I=1,N
            DF(I)=SUM(F(I,1:M)*F(I,1:M))
            DF(I)=SQRT(DF(I)/M)  
            F(I,1:M)=F(I,1:M)/DF(I)
        ENDDO
      ENDIF
      RETURN
      END
!!22222222222222222222222222222222求协方差矩阵A
      SUBROUTINE FORMA(N,M,MNH,F,A)  
   IMPLICIT NONE
!     THIS SUBROUTINE FORMS A BY F  
       INTEGER :: I,J,M,N,MNH
    REAL  :: F(N,M),A(MNH,MNH)
  
      A=0.0  
   IF(M<N)  THEN  
         DO  I=1,M
         DO  J=I,M
         A(I,J)=SUM(F(1:N,I)*F(1:N,J))
         A(J,I)=A(I,J)
         ENDDO
   ENDDO
      ELSE
         DO  I=1,N
         DO  J=I,N
         A(I,J)=SUM(F(I,1:M)*F(J,1:M))
         A(J,I)=A(I,J)
         ENDDO
   ENDDO
      ENDIF
      RETURN
      END
!!333333333333333333333333333333333333333雅可比过关法求特征值特征向量
           SUBROUTINE JCB(N,A,S,EPS)
   IMPLICIT NONE
!     THIS SUBROUTINE COMPUTS EIGENVALUES
!     AND EIGENVECTORS OF A   RETUERN S
      INTEGER   :: I,J,K,N,L ,I1
      REAL   ::A(N,N),S(N,N)
   REAL  :: EPS,G,S1,S2,S3,V1,V2,V3,ST,CT,IP,IQ,U,IQ1
      S=0.
       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)
   print*,"999"  
      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)  
   !PRINT*,V2*V2+U*U,1.-G*G,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
!!444444444444444444444444444444444444444按照特征值大小排序
      SUBROUTINE ARRANG(MNH,A,ER,S)
   IMPLICIT NONE
!     THIS SUBROUTINE PROVIDES A SERIES OF EIGENVALUES
!          FROM MAX TO MIN
     INTEGER  ::  MNH,K1,K2,I ,MNH1
  REAL :: A(MNH,MNH),ER(mnh,6),S(MNH,MNH)
   
   REAL     ::  TR,C
      TR=0.0
   DO  I=1,MNH
      TR=TR+A(I,I)
      ER(I,1)=A(I,I)
      ENDDO

   MNH1=MNH-1
      DO  K1=MNH1,1,-1   
      DO  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 I=1,MNH  
            C=S(I,K2+1)
            S(I,K2+1)=S(I,K2)
            S(I,K2)=C
         ENDDO  
      ENDIF
      ENDDO
      ENDDO  

      ER(1,2)=ER(1,1)
      DO  I=2,mnh
      ER(I,2)=ER(I-1,2)+ER(I,1)
   enddo
      CONTINUE
     
   er(:,5)=er(:,1)*sqrt(2/real(mnh))
   er(:,6)=er(:,5)/TR*100
      ER(:,3)=ER(:,1)/TR*100
      ER(:,4)=ER(:,2)/TR*100

      
      OPEN(119,file="eigenvalues.dat")
      DO i=1,mnh
      WRITE(119,'(2f16.2,4f16.6)') er(i,1),er(i,2),er(i,3),er(i,4),er(i,5),er(i,6)
  
      ENDDO
   CLOSE(119)
      RETURN
      END
!!555555555555555555555555555求Y!给出时间序列和标准化的空间场
      SUBROUTINE TCOEFF(KVT,N,M,MNH,S,F,ER)
   IMPLICIT NONE
!  THIS SUBROUTINE PROVIDES EIGENVECTORS (M.GE.N,SAVED IN S;
!  M.LT.N,SAVED IN F) AND ITS STANDARD TIME COEFFICENTS SERIES (M.GE.N,
!  SAVED IN F; M.LT.N,SAVED IN S)
     INTEGER   :: i,j,k,M,N,JS ,MNH,IS,KVT
  REAL  :: S(MNH,MNH),F(N,M),V(MNH),ER(mnh,6)
     REAL  :: C
      DO  J=1,KVT                     
      C=0.
      C=SUM(S(:,J)*S(:,J))
      C=SQRT(C)
      S(:,J)=S(:,J)/C
      ENDDO
      IF(N.LE.M) THEN
         DO J=1,M
            V(1:N)=F(1:N,J)
            F(1:N,J)=0.         
            DO  IS=1,KVT
            F(IS,J)=SUM(V(1:N)*S(1:N,IS))   
      ENDDO
         ENDDO
       ELSE
         DO I=1,N      
           V(1:M)=F(I,1:M)
           F(I,1:M)=0.  
           DO JS=1,KVT        
           F(I,JS)=SUM(V(1:M)*S(1:M,JS))        
           ENDDO
         ENDDO
        DO  JS=1,KVT
          S(1:M,JS)=S(1:m,JS)*SQRT(ER(JS,1))  
          F(1:N,JS)=F(1:N,JS)/SQRT(ER(JS,1))         
          ENDDO
      ENDIF
     RETURN
      END

islandowner

{:eb303:}{:eb303:}