这个eof程序为什么出不了结果？？

Susie

      PARAMETER(M=52,N=19,MNH=13,KS=1,KV=8,KVT=3)!N为站点数，若N>M，则MNH=M
      DIMENSION F(N,M),A(MNH,MNH),S(MNH,MNH),ER(MNH,4),t(n,m)
        dimension err(kvt,2)
      DIMENSION DF(N),V(MNH),AVF(N)
        real EGVT(N,KVT),ECOF(M,KVT)
      dimension H(N),U(N),VR(N),WT(N),WK(N)
      dimension WBT(M),WBK(M),R(KVT),ST(KVT)
        dimension fc(kvt)
        INTEGER iSTa(n)
CCCCCCCCCCCCCCCCCCCCCC   input data CCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        open(2,file='D:\lw\eof\winter.txt')
        do i=1,n
        do it=1,m
      READ(2,*)f(i,it)
        enddo
        enddo   
      print*,'read OK OK oK'
        close(2)

CCCCCCCCCCCCCCCC INPUT DATA CCCCCCCCCCCCCCCCCCC
      CALL TRANSF(N,M,F,AVF,DF,KS)
      CALL FORMA(N,M,MNH,F,A)
      CALL JCB(MNH,A,S,0.00001)
          write(*,*)'call 3 ok'

      CALL ARRANG(KV,MNH,A,ER,S,TR0)
        call jy(KV,M,ER)
        write(*,*)'jian yan OK'
      CALL TCOEFF(KVT,KV,N,M,MNH,S,F,V,ER)
C      CALL OUTER(KV,KVT,ER,R)
      CALL OUTER(KV,KVT,MNH,ER,R)
      CALL OUTVT(KVT,N,M,MNH,S,F,EGVT,ECOF)
cccccccccccccccccc output unrotatted time cooficience cccccccccc
       do 131 js=1,kvt
         fc(js)=0.0
       do 132 j=1,M
132    fc(js)=fc(js)+ecof(j,js)*ecof(j,js)/float(M)
        fc(js)=sqrt(fc(js))
131     print*,'fc= ',js,'  ',fc(js)


      open(16,file='D:\lw\eof\winteruntmc.grd',form='binary')
       do 348 j=1,m
348   write(16)(ecof(j,is)/fc(is),is=1,kvt)
       close(16)
cccccccccccccc output unrotatted eigenvectors cccccccccccccc

      open(25,file='D:\lw\eof\winterunvec.grd',form='binary')
       write(25) ((fc(js)*egvt(i,js),i=1,n),js=1,kvt)
       close(25)
       print*,' EOF OK !'
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      CALL rotation(EGVT,ECOF,R,H,ST,N,M,KVT,U,VR,WT,WK,WBT,WBK)
      call rarrange(KVT,N,M,TR0,EGVT,ECOF,ERr)
cccccccccccccccccc output rotatted time cooficience cccccccccc
       do 133 js=1,kvt
         fc(js)=0.0
       do 134 j=1,M
134    fc(js)=fc(js)+ecof(j,js)*ecof(j,js)/float(M)
        fc(js)=sqrt(fc(js))
133     print*,'fc= ',js,'  ',fc(js)


      open(26,file='D:\lw\eof\winterrtmc.grd',form='binary')
      do 449 j=1,M
449  write(26)(ecof(j,is)/fc(is),is=1,kvt)
      close(26)
cccccccccccccc output rotatted eigenvectors cccccccccccccc

        OPEN(25,file='D:\lw\eof\winterrvec.grd',
     &form='binary')
          write(25) ((fc(js)*egvt(i,js),i=1,n),js=1,kvt)
        close(25)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      STOP
      END
!----------------------将资料标准化或距平----------------------!
!ks=1将资料标准化,即均值为零,方差为1;ks=0将资料化为距平值；ks=-1原始资料不变化
!avf(i)为第i个格点的多年均值，df(i)为方差；m为时间长度，n为格点数。
      SUBROUTINE TRANSF(N,M,F,AVF,DF,KS)
      DIMENSION F(N,M),AVF(N),DF(N)
      DO  I=1,N
      AVF(I)=0.0
      DF(I)=0.0
      enddo
      IF(KS.eq.0)then
      DO  I=1,N
       DO J=1,M
       AVF(I)=AVF(I)+F(I,J)
         enddo
       AVF(I)=AVF(I)/M
       DO  J=1,M
       F(I,J)=F(I,J)-AVF(I)
       enddo
        enddo
      elseif(KS.EQ.-1) then
      RETURN
      ELSE
        DO  I=1,N
      DO  J=1,M
      AVF(I)=AVF(I)+F(I,J)
        enddo
      AVF(I)=AVF(I)/M
      DO j=1,M
      F(I,J)=F(I,J)-AVF(I)
        enddo
        enddo
      DO  I=1,N
      DO  J=1,M
      DF(I)=DF(I)+F(I,J)*F(I,J)
        enddo
      DF(I)=SQRT(DF(I)/M)
      DO  J=1,M
      F(I,J)=F(I,J)/DF(I)
      enddo
        enddo
      ENDIF
      RETURN
      END
! ------------- 计算相关矩阵A(mnh,mnh)------------!
      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

!----------计算相关矩阵A(mnh,mnh)的特征值和特征向量---------!
      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/real(N)*S1
      S3=S1
      L=0
  50  S3=S3/real(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

!---------------特征值从大到小排序---------------!
      SUBROUTINE ARRANG(KV,MNH,A,ER,S,TR)
c      DIMENSION A(MNH,MNH),ER(KV,4),S(MNH,MNH)
      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
      RETURN
      END

!--------------------计算标准化特征向量及其时间系数-----------------------!
      SUBROUTINE TCOEFF(KVT,KV,N,M,MNH,S,F,V,ER)
C     THIS SUBROUTINE PROVIDES STANDARD EIGENVECTORS (M.GE.N,SAVED IN S;
C          M.LT.N,SAVED IN F) AND ITS TIME COEFFICENTS SERIES (M.GE.N,
C          SAVED IN F; M.LT.N,SAVED IN S)
c      DIMENSION S(MNH,MNH),F(N,M),V(MNH),ER(KV,4)
      DIMENSION S(MNH,MNH),F(N,M),V(MNH),ER(MNH,4)
      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
  160 S(I,J)=S(I,J)/C
  360 CONTINUE
      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 380 IS=1,KVT
      DO 380 I=1,N
  380 F(IS,J)=F(IS,J)+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
      S(J,JS)=S(J,JS)*SQRT(ER(JS,1))
  420 CONTINUE
      DO 430 I=1,N
      F(I,JS)=F(I,JS)/SQRT(ER(JS,1))
  430 CONTINUE
      ENDIF
      RETURN
      END

C      SUBROUTINE OUTER(KV,KVT,ER,R)
      SUBROUTINE OUTER(KV,KVT,MNH,ER,R)
C     THIS SUBROUTINE PRINTS ARRAY ER
C     ER(KV,1) FOR  SEQUENCE OF EIGENVALUE FROM BIG TO SMALL
C     ER(KV,2) FOR  EIGENVALUE FROM BIG TO SMALL
C     ER(KV,3) FOR  SMALL LO=(LAMDA/TOTAL VARIANCE)
C     ER(KV,4) FOR  BIG LO=sum OF SMALL LO)
C      DIMENSION ER(KV,4),R(KVT)
      DIMENSION ER(MNH,4),R(KVT)
      open(26,file='D:\lw\eof\unrvalue.dat')
      WRITE(26,510)
  510 FORMAT(/30X,'EIGENVALUE AND ANALYSIS ERROR',/)
      WRITE(26,520)
  520 FORMAT(10X,1HH,8X,5HLAMDA,10X,6HSLAMDA,11X,2HPH,12X,3HSPH)
      WRITE(26,530) (IS,(ER(IS,J),J=1,4),IS=1,KV)
  530 FORMAT(1X,I10,4F15.5)
      WRITE(26,540)
  540 FORMAT(//)
      DO I=1,KVT
      R(I)=ER(I,1)
      ENDDO
      RETURN
      END
      SUBROUTINE OUTVT(KVT,N,M,MNH,S,F,EGVT,ECOF)
C     THIS SUBROUTINE PRINTS STANDARD EIGENVECTORS
C          AND ITS TIME-COEFFICENT SERIES
      DIMENSION F(N,M),S(MNH,MNH),EGVT(N,KVT),ECOF(M,KVT)
      DO 550 I=1,N
      IF(M.GE.N) THEN
      DO 11 JS=1,KVT
      EGVT(I,JS)=S(I,JS)
   11 CONTINUE
      ELSE
      DO 12 JS=1,KVT
      EGVT(I,JS)=F(I,JS)
   12 CONTINUE
      ENDIF
  550 CONTINUE

      DO 600 J=1,M
      IF(M.GE.N) THEN
      DO 13 IS=1,KVT
      ECOF(J,IS)=F(IS,J)
   13 CONTINUE
      ELSE
      DO 14 IS=1,KVT
      ECOF(J,IS)=S(J,IS)
   14 CONTINUE
      ENDIF
  600 CONTINUE
      RETURN
      END


      subroutine rotation(a,b,r,h,st,mm,nn,np,u,vr,wt,wk,wbt,wbk)
c-----------------------------------------------------------------c
c     a(mm,np): eigenvector in unrotated origin field             c
c     b(nn,np): time coefficient(PC) in unrotated origin field    c
c     r(np):    eigenvalue in unrotated origin field              c
c-----------------------------------------------------------------c

      dimension a(mm,np),b(nn,np),h(mm),u(mm),vr(mm),wt(mm),wk(mm)
      dimension wbt(nn),wbk(nn),r(np),st(np),vrlv(40)
      integer t,k
      do 333 j=1,np
      do 333 i=1,mm
      a(i,j)=a(i,j)*sqrt(r(j))
333  continue

      do 321 j=1,np
      do 321 it=1,nn
      b(it,j)=b(it,j)/sqrt(r(j))
321  continue

      do 20 i=1,mm
      h(i)=0.0
      do 21 j=1,np
21   h(i)=h(i)+a(i,j)**2
20   continue

      do 23 i=1,mm
23   h(i)=sqrt(h(i))

      lll=0
      vrlv0=0.0
      do 230  k=1,np
      g1=0.0
      g2=0.0
      do  434  i=1,mm
      g1=g1+(a(i,k)**2/h(i)**2)**2/real(mm)
434  g2=g2+(a(i,k)**2/h(i)**2)/real(mm)
      g2=g2**2
      vrlv0=vrlv0+g1-g2
230  continue
800  continue
      do 505 t=1,np
      do 505 k=1,np
      if(t.ge.k) goto 505
      call rot(a,b,h,t,k,mm,nn,np,u,vr,wt,wk,wbt,wbk)
505  continue

      lll=lll+1
      vrlv(lll)=0.0
      do 530  k=1,np
      g1=0.0
      g2=0.0
      do  534  i=1,mm
      g1=g1+(a(i,k)**2/h(i)**2)**2/real(mm)
534  g2=g2+(a(i,k)**2/h(i)**2)/real(mm)
      g2=g2**2
      vrlv(lll)=vrlv(lll)+g1-g2
530  continue
      if(lll.lt.2) goto 800
      ci=(vrlv(lll)-vrlv(lll-1))/vrlv0
      if(ci.lt.0.01) goto 600
      if(lll.ge.40) goto 600
      goto 800
600  continue
      write(6,538)
      write(6,549)
      write(6,539) vrlv0,(vrlv(k),k=1,lll)
538  format(//1x,'sum of variance of squared element of each colum of
     1 rotated loding martix ')
549  format(1x,'at each circulation')
539  format((1x,10f8.3))
      return
      end

      subroutine rot(a,b,h,t,k,mm,nn,np,u,vr,wt,wk,wbt,wbk)
      dimension a(mm,np),b(nn,np),h(mm),u(mm),vr(mm),wt(mm),wk(mm),
     1 wbt(nn),wbk(nn)
      integer t,k
c      write(6,*)((b(it,j),it=1,nn),j=1,np)
      do 25 i=1,mm
      u(i)=(a(i,t)/h(i))**2-(a(i,k)/h(i))**2
25   vr(i)=2.0*(a(i,t)/h(i))*(a(i,k)/h(i))
      c=0.0
      d=0.0
      s=0.0
      g=0.0
      do 30 i=1,mm
      c=c+(u(i)**2-vr(I)**2)
      d=d+2.0*u(i)*vr(i)
      s=s+u(i)
30   g=g+vr(i)
      tg1=d-2.0*s*g/mm
      tg2=c-(s**2-g**2)/mm
      fi=atan2(tg1,tg2)/4.0
      do 75 i=1,mm
      wt(i)=a(i,t)
75   wk(i)=a(i,k)
      do 84 kk=1,nn
      wbt(kk)=b(kk,t)
84   wbk(kk)=b(kk,k)
      do 78 i=1,mm
      a(i,t)=wt(i)*cos(fi)+wk(i)*sin(fi)
78   a(i,k)=wt(i)*(-sin(fi))+wk(i)*cos(fi)
      do 89 it=1,nn
      b(it,t)=wbt(it)*cos(fi)+wbk(it)*sin(fi)
89   b(it,k)=wbt(it)*(-sin(fi))+wbk(it)*cos(fi)
      return
      end

      SUBROUTINE rarrange(KVT,N,M,TR,a,b,ER)
C     THIS SUBROUTINE PROVIDES A SERIES OF EIGENVALUES
C          FROM MAX TO MIN
      DIMENSION a(N,KVT),b(M,KVT),ER(KVT,2)
      do 802 I=1,KVT
      ER(I,1)=0.0
      do 803 J=1,N
803  ER(I,1)=ER(I,1)+a(J,I)**2
802  continue
      KVT1=KVT-1
      DO 210 K1=KVT1,1,-1
      DO 210 K2=K1,KVT1
      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,N
      C=a(I,K2+1)
      a(I,K2+1)=a(I,K2)
      a(I,K2)=C
  205 CONTINUE
      DO 305 I=1,M
      C=b(I,K2+1)
      b(I,K2+1)=b(I,K2)
      b(I,K2)=C
  305 CONTINUE
      ENDIF
  210 CONTINUE
      DO 230 I=1,KVT
      ER(I,2)=ER(I,1)/TR
  230 CONTINUE
      ee=0.0
      do i=1,KVT
      ee=ee+ER(i,2)
      enddo
      open(16,file='D:\lw\eof\rvalue.dat')
        do I=1,KVT
      write(16,*) er(i,1),ER(I,2)      
        enddo
      write(16,*) ee
      RETURN
      END

      SUBROUTINE jy(KV,M,ER)
      DIMENSION ER(KV,4),wc(kv)
        open(1,file='D:\lw\eof\jy.dat')
        do i=1,kv
        wc(i)=er(i,1)*sqrt(2./real(m))
        enddo
        do i=1,kv-1
        a=er(i,1)-er(i+1,1)
        if(a.ge.wc(i))then
      b=i+1
        write(1,*)i,b,'you jia zhi'
        endif
        enddo
        close(1)
        end

言深深

额，程序太长，没法看···换一个别的有使用说明的吧···不熟别的，你的read部分就看不懂···所有数据（需要被分解的时空场）就一列还是怎么回事？

Susie

言深深 发表于 2013-5-8 16:07

                               
登录/注册后可看大图

额，程序太长，没法看···换一个别的有使用说明的吧···不熟别的，你的read部分就看不懂···所有数据 ...

PROGRAM EOFPW
C     M:LENTH 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特征值的个数
! KVT:NUMBER OF EIGENVECTORS AND TIME SERIES WILL BE OUTPUT输出旋转特征向量及其时间系数的个数
C     MNH=MIN(M,N)
C     EGVT==EIGENVACTORS旋转特征向量, ECOF==TIME COEFFICIENTS FOR EGVT旋转特征向量的时间系数.
C     ER(KV,1)====LAMDA 特征值 ， LAMDA======EIGENVALUE
C     ER(KV,2)====ACCUMULATE LAMDA 逐步累积特征值
! ER(KV,3)====THE sum OF COMPONENTS VECTORS PROJECTED ONTO EIGENVACTOR.各分量在特征向量里所占比重
C     ER(KV,4)====ACCUMULATE ER(KV,3)逐步累积分量在特征向量里所占比重

Susie

言深深 发表于 2013-5-8 16:07

                               
登录/注册后可看大图

额，程序太长，没法看···换一个别的有使用说明的吧···不熟别的，你的read部分就看不懂···所有数据 ...

学长你不要多看，帮我看看这个MNH该怎么设置。。。谢谢了

Susie

言深深 发表于 2013-5-8 16:07

                               
登录/注册后可看大图

额，程序太长，没法看···换一个别的有使用说明的吧···不熟别的，你的read部分就看不懂···所有数据 ...

是只有一列啊，我把季节平均单拿出来存在winter.txt里了，一共是52年19个站点的

言深深

Susie 发表于 2013-5-8 16:23

                               
登录/注册后可看大图

学长你不要多看，帮我看看这个MNH该怎么设置。。。谢谢了

你的MNH是19，取的是M，N中较小的一个···

言深深

Susie 发表于 2013-5-8 16:29

                               
登录/注册后可看大图

是只有一列啊，我把季节平均单拿出来存在winter.txt里了，一共是52年19个站点的

你是否确定你的数据read正确了，要求是一个时空场的，也就是read的矩阵里面包括的是很多个时刻的所有空间数据。

如果以上没有问题，那么你的程序是否运行结束，没有运行结束的话，错在哪里？运行结束了请查看'D:\lw\eof\文件夹下的生成文件，必须要有此文件才对！

Susie

言深深 发表于 2013-5-8 18:07

                               
登录/注册后可看大图

你是否确定你的数据read正确了，要求是一个时空场的，也就是read的矩阵里面包括的是很多个时刻的所有空间 ...

谢谢这位哥哥指导，解决了解决了，果然是read的问题

言深深

Susie 发表于 2013-5-9 21:14

                               
登录/注册后可看大图

谢谢这位哥哥指导，解决了解决了，果然是read的问题

解决了就好···

一闭一睁

没看懂啊！！！！！！