用Fortran做eof原始和据平贡献差不多，哪里错了？

weijiade99

C        THIS IS A PROGRAM FOR CACULATING EMPIRICAL ORTHOGONAL FUNCTION
        PROGRAM EOF
        PARAMETER(N=50,M=160,JOB=1)
        DIMENSION F(N,M),T(N,M),IX(M),V(M,M),A(M,M),
     &        D(M),W(M),H(N,M),B(M),Z(M),H1(M),V1(M,M),FM(M)
C        ***********************************************
C        * N:   SAMPLE SIZE                            *
C        * M:   NUMBER OF STATION OR GRID POINT        *
C        * JOB: CONTRAL PARAMETER;                     *
C       *      WHEN JOB=0 ORIGINAL DATA;              *
C        *      WHEN JOB=1 DEPATURE DATA;              *
C        *      WHEN JOB=2 STANDARDIZE DATA.           *
C        * F(N,M): ORIGINAL DATA FIELD                 *
C        * D(M): EIGENVALUE                            *
C        * V(M,M): CHARACTERISTIC VECTOR               *
C        * T(N,M): TIME COE FFICIENT                   *
C        ***********************************************
        OPEN(11,FILE='F:\zhang3\qw.txt')
        READ(11,*)((F(I,J),J=1,M),I=1,N)
        OPEN(6,FILE='F:\zhang3\14.txt',STATUS='NEW')
*************************
      open(7,FILE='F:\zhang3\er.grd',form='binary')
        open(8,FILE='F:\zhang3\et.grd',form='binary')
*****************************       
        CALL SEOF(M,N,JOB,F,H,T,IX,V,A,D,W,FM,B,Z,H1,V1)
        STOP
        END
        SUBROUTINE SEOF(N,LL,JOB,F,H,T,IX,V,A,D,W,FM,B,Z,H1,V1)
        DIMENSION F(LL,N),T(LL,N),IX(N),V(N,N),A(N,N),D(N),
     &        W(N),H(LL,N),B(N),Z(N),H1(N),V1(N,N),FM(N)
        IF(JOB.EQ.0)THEN
        GOTO 1
        ELSE IF(JOB.EQ.1)THEN
        GOTO 2
        ELSE
        GOTO 3
        END IF
  1        DO 10 I=1,LL
        DO 10 J=1,N
  10        H(I,J)=F(I,J)
        GOTO 111
  2        CALL DEP(N,LL,F,H,W)
        GOTO 111
  3        CALL NOR(N,LL,F,H,FM,W)
  111        CONTINUE
        DO 20 I=1,N
        DO 20 J=1,N
        A(I,J)=0.0
        DO 17 K=1,LL
  17        A(I,J)=A(I,J)+H(K,I)*H(K,J)
        A(J,I)=A(I,J)
  20        CONTINUE
        WRITE(6,200)
  200        FORMAT(30X,'---------------EOF---------------')         
        WRITE(6,900)
  900        FORMAT(//20X,'****** CORVARINCE ******')
        DO 800 I=1,N
        WRITE(6,801)I,(A(I,J),J=1,N)
  801        FORMAT(2X,I3/(2X,8F15.4))
  800        CONTINUE
        CALL JACOBI (N,.TRUE.,D,IRT,B,Z,V,A)
        DO 400 I=1,N
  400        IX(I)=I
        DO 401 I=1,N-1
        DO 402 J=I+1,N
        IF(ABS(D(I)).GE.ABS(D(J))) GOTO 402
        W1=D(I)
        D(I)=D(J)
        D(J)=W1
        K=IX(I)
        IX(I)=IX(J)
        IX(J)=K
  402        CONTINUE
  401        CONTINUE
        WRITE (6,104)
  104         FORMAT (//20X,'******* EIGENVALUE VALUE ******')
        WRITE(6,500)(IX(I),I=1,N)
  500        FORMAT(/20X,'MAX-MIN'/(12X,20I5))
        WRITE(6,501)(D(I),I=1,N)
  501        FORMAT(2X,8F15.5)
        DO 502 J=1,N
        ILW=IX(J)
        DO 503 I=1,N
  503        V1(I,J)=V(I,ILW)
  502        CONTINUE
        WRITE (6,105)
  105         FORMAT (//20X, '****** CHARACTERISTIC VECTOR ******')
        DO 600 I=1,N
        WRITE(6,601)(V1(I,J),J=1,10)
*************************************
        write(7) (V1(i,j),j=1,2)
***********************************
  601        FORMAT(10X,10F10.2)
  600        CONTINUE
        DO 25 L=1,LL
        DO 25 J=1,N
        T(L,J)=0.
        DO 31 I=1,N
  31        T(L,J)=T(L,J)+F(L,I)*V1(I,J)
  25        CONTINUE
        WRITE (6,106)
  106        FORMAT (///20X, '****** TIME COEFFICIENT ******')
        DO 700 I=1,LL
        WRITE(6,701)(T(I,J),J=1,10)
*****************************
        write(8)(T(I,J),J=1,2)
*****************************************
  701        FORMAT(10X,10F10.2)
  700        CONTINUE
        DO 705 I=1,LL
  705        CONTINUE
        AP=0.0
        DO 702 I=1,N
  702        AP=AP+D(I)
        DO 703 I=1,N
        AP1=0.0
        DO 704 J=1,I
  704        AP1=AP1+D(J)
  703        H1(I)=1.0*AP1/AP
        WRITE(6,107)
  107        FORMAT(//20X,'****** ACCUMULATE PROPORTION ****** ')
        WRITE(6,108)(H1(I),I=1,N)
  108        FORMAT(/10X,10F10.5)
        RETURN
        END
        SUBROUTINE DEP(N,LL,F,H,W)
        DIMENSION F(LL,N),H(LL,N),W(N)
        DO 10 J=1,N
        W(J)=0.0
        DO 20 I=1,LL
  20        W(J)=W(J)+F(I,J)
        W(J)=W(J)/FLOAT(LL)
  10        CONTINUE
        DO 30 J=1,N
        DO 30 I=1,LL
        H(I,J)=F(I,J)-W(J)
  30        CONTINUE
        RETURN
        END
        SUBROUTINE NOR(N,LL,F,H,FM,W)
        DIMENSION F(LL,N),H(LL,N),FM(N),W(N)
        DO 10 J=1,N
        W(J)=0.0
        DO 20 I=1,LL
  20        W(J)=W(J)+F(I,J)
        FM(J)=W(J)/FLOAT(LL)
        W(J)=0.0
        DO 30 I=1,LL
  30        W(J)=W(J)+(F(I,J)-FM(J))**2
        W(J)=SQRT(W(J)/FLOAT(LL))
        DO 40 I=1,LL
  40        H(I,J)=(F(I,J)-FM(J))/W(J)
  10        CONTINUE
        RETURN
        END
        SUBROUTINE JACOBI(N,EV,D,IRT,B,Z,V,A)
        DIMENSION D(N),B(N),Z(N),V(N,N),A(N,N)
        LOGICAL EV
        IF(. NOT. EV) GOTO 20
        DO 10 K=1,N
        DO 10 L=1,N
        IF (K-L) 5,6,5
   6        V (K,K)=1.0
        GOTO 10
   5        V(K,L)=0.0
  10        CONTINUE
  20        DO 30 K=1,N
        B(K)=A(K,K)
        D(K)=B(K)
  30        Z(K)=0.0
        IRT=0
        DO 200 I=1,50
        SM=0.0
        N1=N-1
        DO 40 K=1,N1
        K1=K+1
        DO 40 L=K1,N
  40        SM=SM+ABS(A(K,L))
        IF (SM)101,210,101
  101        IF(I-4)50,60,60
  50        TRESH=0.2*SM/(FLOAT(N)*FLOAT(N))
        GOTO 70
  60        TRESH=0.0
  70         DO 180 K=1,N1
        K1=K+1
        DO 180 L=K1,N
        G=100.0*ABS(A(K,L))
        IF (I. GT. 4. AND. ABS(D(K))+G. EQ. ABS(D(K)).
     *  AND. ABS(D(L))+G.EQ. ABS(D(L))) GOTO 170
        IF (ABS(A(K,L)). LE . TRESH) GOTO 180
        H=D(L)-D(K)
        IF (ABS(H)+G.EQ.ABS(H)) GOTO 80
        THETA=0.5*H/A(K,L)
        T=1.0/(ABS(THETA)+SQRT(1.0+THETA*THETA))
        IF (THETA. LT. 0.0) T=-T
        GOTO 90
  80        T=A(K,L)/H
  90    C=1.0/SQRT(1.0+T*T)
        S=T*C
        H=T*A(K,L)
        Z(K)=Z(K)-H
        Z(L)=Z(L)+H
        D(K)=D(K)-H
        D(L)=D(L)+H
        A(K,L)=0.0
        KM1=K-1
        IF (KM1) 120,120,100
100        DO 110 J=1,KM1
        G=A(J,K)
        H=A(J,L)
        A(J,K)=C*G-S*H
110        A(J,L)=S*G+C*H
120        L1=L-1
        IF (L1-K1)140,300,300
300        DO 130 J=K1,L1
        G=A(K,J)
        H=A(J,L)
        A(K,J)=C*G-S*H
130         A(J,L)=S*G+C*H
140        L1=L+1
        IF (L1. GT. N) GOTO 150
        DO 145 J=L1,N
        G=A(K,J)
        H=A(L,J)
        A(K,J)=C*G-S*H
145        A(L,J)=S*G+C*H
150        IF (.NOT. EV) GOTO 160
        DO 155 J=1,N
        G=V(J,K)
        H=V(J,L)
        V(J,K)=C*G-S*H
155        V(J,L)=S*G+C*H
160        IRT=IRT+1
        GOTO 180
170        A(K,L)=0.0
180        CONTINUE
        DO 200 K=1,N
        D(K)=B(K)+Z(K)
        B(K)=D(K)
200        Z(K)=0.0
210        RETURN
        END

weijiade99

C:\Users\Administrator\Desktop

weijiade99

原数据http://bbs.06climate.com/forum.php?mod=attachment&aid=ODUyNzV8M2VjZDk4NjcxZjBhMWVhY2ExMjg1NTZjYzY3YTk0MzZ8MTc2MTkyODQwNg%3D%3D&request=yes&_f=.jpg