关于EOF

吴阿好大人

      program main
parameter(m=37,n=51,mnl=n,np=20,ll=m)
      parameter(ks=1,undef=32767.0)
!-----input array
      dimension f0(m,n),F(M,N)
!-----work array
      dimension gvt(ll,mnl),rgvt(ll,np),cof(mnl,n),rcof(np,n)
!-----output arrays
      dimension er(mnl,4),egvt(m,mnl),ecof(mnl,n)
      dimension rer(np,2),regvt(m,np),recof(np,n)
       integer sta(m)
  real jd(m),wd(m)
character aa
c-----read data
      open(1,file='1.txt')
read(1,*)aa
  do I=1,M
c   read(1,101)(f0(i,j),J=1,N)
   read(1,*)sta(i),jd(i),wd(i),(f0(i,j),j=1,n)
  end do
!101  FORMAT(33F10.9)
close(1)
c-----Remove the terrain or missing value.
      l=0
do j=1,m
   if(f0(j,1).ne.undef)then
          l=l+1
     do k=1,n
            f(l,k)=f0(j,k)
     enddo
   endif
enddo
      write(*,*)'ll=',l
c-----call the subroutine
      call reof(ll,n,mnl,np,f,ks,er,gvt,ecof,rer,rgvt,recof)
c-----Add the terrain or missing value.
      l=0
do i=1,m
   if(f0(i,1).ne.undef)then
     l=l+1
     do k=1,n
       egvt(i,k)=gvt(l,k)
     end do
  do k=1,np
       regvt(i,k)=rgvt(l,k)
     end do
   else
     do k=1,n
       egvt(i,k)=undef
     end do
     do k=1,np
       regvt(i,k)=undef
     end do
   endif
enddo
!-----output the result        
c-----output the error
      OPEN(10,file='er.dat')
    !  write(10,*)'error=',sum-tr
      write(10,*)'EOF lanmda (eigenvalues) from big to small'
write(10,*)(er(i,1),i=1,mnl)
write(*,*)'lamda='
      write(*,*)(er(i,1),i=1,mnl)
write(10,*)'EOF accumulated eigenvalues from big to small'
write(10,*)(er(i,2),i=1,mnl)
write(10,*)'EOF explained variances'
write(10,*)(er(i,3),i=1,mnl)
write(10,*)'EOF accumulated explained variances'
write(10,*)(er(i,4),i=1,mnl)
write(10,*)'REOF explained variances'
write(10,*)(rer(i,1),i=1,np)
write(10,*)'REOF accumulated explained variances'
write(10,*)(rer(i,2),i=1,np)
c-----output eigenvactors matrix of EOF
      OPEN(11,FILE='egvt.dat',form='binary')
do i=1,m
write(11)(egvt(i,j),j=1,np)
end do
close(11)
      OPEN(11,FILE='egvt.txt')
do i=1,m
write(11,102)sta(i),jd(i),wd(i),(egvt(i,j),j=1,np)
enddo
102   format(i7,2x,2(f7.2,2x),20f8.2)
close(11)
c-----output time coefficients matrix of EOF
      open(14,file='ecof.dat',form='binary')
      do i=1,n
write(14)(ecof(j,i),j=1,np)
end do
close(14)
      open(14,file='ecof.txt')
      do i=1,n
write(14,103)(ecof(j,i),j=1,np)
end do
103 format(20f8.2)
close(14)
c-----output loading vectors of REOF
      OPEN(21,FILE='regvt.dat',form='binary')
      do i=1,m
    write(21)(regvt(i,j),j=1,np)
end do
close(21)
      OPEN(21,FILE='regvt.txt')
      do i=1,m
write(21,102)sta(i),jd(i),wd(i),(regvt(i,j),j=1,20)
      end do
!103   format(i7,2x,2(f7.2,2x),20f10.2)
close(21)
c-----output time coefficients matrix of REOF
      OPEN(22,FILE='recof.dat',form='binary')
      do j=1,n
      write(22)(recof(i,j),i=1,np)
end do
close(22)
      OPEN(22,FILE='recof.txt')
      do j=1,n
      write(22,103)(recof(i,j),i=1,np)
end do
!104   format(20f10.4)
close(22)
      stop
end
c-----
      subroutine reof(m,n,mnl,np,f,ks,er,egvt,ecof,rer,regvt,recof)
c-----input array
      dimension f(m,n)
c-----work arrays
      dimension cov(mnl,mnl),s(mnl,mnl),d(mnl)
c-----output arrays
      dimension er(mnl,4),egvt(m,mnl),ecof(mnl,n)
      dimension rer(np,2),regvt(m,np),recof(np,n)
c---- Preprocessing data
      call transf(m,n,f,ks)
c---- Crossed product matrix of the data f
      call crossproduct(m,n,mnl,f,cov,sum)
c---- Eigenvalues and eigenvectors by the Jacobi method
      eps=1e-7
      call jcb(mnl,cov,s,d,eps)
c---- Specified eigenvalues
      call arrang(m,n,mnl,d,s,er,tr)
write(*,*)'error=',sum-tr
c---- Normalized eignvectors and their time coefficients
      call tcoeff(m,n,mnl,f,s,er,egvt,ecof)
write(*,*)'EOF accomplished'
c-----linear transposed (rotate)
      call rotaor(m,n,mnl,np,egvt,ecof,tr,rer,regvt,recof)
write(*,*)'REOF accomplished'
c---- Specified eignvectors and time coefficients with
c     explained variances of REOF (if you need)
      call arrange2(m,n,np,regvt,recof,rer)
      return
      end
      subroutine transf(m,n,f,ks)
!-----input array
      dimension f(m,n)
!-----work array
      dimension fw(n)
      if(ks.eq.-1)return
      if(ks.eq.0)then                !anomaly of f
        do i=1,m
          do j=1,n
            fw(j)=f(i,j)
          enddo
          call meanvar(n,fw,af,sf,vf)
          do j=1,n
            f(i,j)=f(i,j)-af
          enddo  
        enddo
        return
      endif
      if(ks.eq.1)then                 !normalizing f
        do i=1,m
          do j=1,n
            fw(j)=f(i,j)
          enddo
          call meanvar(n,fw,af,sf,vf)
          do j=1,n
            f(i,j)=(f(i,j)-af)/sf
          enddo
        enddo
      endif
      return
      end
c-----
      subroutine crossproduct(m,n,mnl,f,cov,sum)
!-----input array
      dimension f(m,n)
!-----work array
dimension cov(mnl,mnl)
      if(n-m) 10,50,50
  10  do 20 i=1,mnl
      do 20 j=i,mnl
        cov(i,j)=0.0
        do is=1,m
          cov(i,j)=cov(i,j)+f(is,i)*f(is,j)
        enddo
        cov(j,i)=cov(i,j)
  20  continue
      sum=0.
do i=1,mnl
sum=sum+cov(i,i)
end do
sum=sum/(mnl*1.)
      return
  50  do 60 i=1,mnl
      do 60 j=i,mnl
        cov(i,j)=0.0
        do js=1,n
          cov(i,j)=cov(i,j)+f(i,js)*f(j,js)
        enddo
        cov(j,i)=cov(i,j)
  60  continue
      sum=0.
do i=1,mnl
sum=sum+cov(i,i)
end do
sum=sum/(mnl*1.)
      return
      end
c-----
      subroutine jcb(m,a,s,d,eps)
!-----input array
      dimension a(m,m)
!-----output arrays
dimension s(m,m),d(m)
      do 30 i=1,m
      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,m
        i1=i-1
        do 40 j=1,i1
  40      g=g+2.*a(i,j)*a(i,j)
      s1=sqrt(g)
      s2=eps/float(m)*s1
      s3=s1
      l=0
  50  s3=s3/float(m)
  60  do 130 iq=2,m
        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,m
          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,m
          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
      do 160 i=1,m
        d(i)=a(i,i)
  160 continue
      return
      end
c-----
      subroutine arrang(m,n,mnl,d,s,er,tr)
!-----input arrays
      dimension d(mnl),s(mnl,mnl)
!-----output array
dimension er(mnl,4)
tr=0.0
      do 10 i=1,mnl
        tr=tr+d(i)/(mnl*1.)
        er(i,1)=d(i)/(mnl*1.)
  10  continue
      mnl1=mnl-1
      do 20 k1=mnl1,1,-1
      do 20 k2=k1,mnl1
        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 15 i=1,mnl
            c=s(i,k2+1)
            s(i,k2+1)=s(i,k2)
            s(i,k2)=c
  15      continue
        endif
  20  continue
      er(1,2)=er(1,1)
      do 30 i=2,mnl
        er(i,2)=er(i-1,2)+er(i,1)
  30  continue
      do 40 i=1,mnl
        er(i,3)=er(i,1)*100./tr
        er(i,4)=er(i,2)*100./tr
  40  continue
      return
      end
c-----
      subroutine tcoeff(m,n,mnl,f,s,er,egvt,ecof)
!-----input arrays
      dimension f(m,n),s(mnl,mnl),er(mnl,4)
!-----work arrays
      dimension v(m),v1(n)
!-----output arrays
dimension egvt(m,mnl),ecof(mnl,n)
      do j=1,mnl
        do i=1,m
          egvt(i,j)=0.
        enddo
        do i=1,n
          ecof(j,i)=0.
        enddo
      enddo
c-----Normalizing the input eignvectors (c=sqrt(lamda))
      do 10 j=1,mnl
        c=0.
        do i=1,mnl
          c=c+s(i,j)*s(i,j)
        enddo
        c=sqrt(c)
        do i=1,mnl
          s(i,j)=s(i,j)/c
        enddo
  10  continue
c-----(m<n)
      if(m.le.n) then      
   do js=1,mnl
          do i=1,m
            egvt(i,js)=s(i,js)
          enddo
        enddo
        do 30 j=1,n
          do i=1,m
            v(i)=f(i,j)
          enddo
          do is=1,mnl
            do i=1,m
              ecof(is,j)=ecof(is,j)+v(i)*egvt(i,is)
            enddo
          enddo
  30    continue
      else
c-----(m>n)
        do 40 i=1,m
          do j=1,n
            v1(j)=f(i,j)
          enddo
          do js=1,mnl
          do j=1,n
            egvt(i,js)=egvt(i,js)+v1(j)*s(j,js)
          enddo
          enddo
  40    continue
        do 50 js=1,mnl
          do j=1,n
            ecof(js,j)=s(j,js)*sqrt(abs(er(js,1)))
          enddo
          do i=1,m
            egvt(i,js)=egvt(i,js)/sqrt(abs(er(js,1)))
          enddo
  50    continue
      do js=1,mnl
   do i=1,m
          egvt(i,js)=(sqrt(abs(er(js,1)/(mnl*1.0))))*egvt(i,js)
   end do
   do i=1,n
     ecof(js,i)=(ecof(js,i))/(sqrt(abs(er(js,1)/(mnl*1.0))))
   end do
end do
end if
      return
      end
c-----
subroutine rotaor(m,n,mnl,np,egvt,ecof,tr,rer,regvt,recof)
c-----input arrays
dimension egvt(m,mnl),ecof(mnl,n)  
c-----work arrays
dimension h(m),st(np),vrlv(50),er(mnl,4),s(m,np)
c-----output arrays  
      dimension regvt(m,np),recof(np,n),rer(np,2)
      integer t,k
      write(*,*)'rotate beginning'
c-----Take out forward NP eigenvectors to rotate
      do j=1,np
        do i=1,m
     regvt(i,j)=egvt(i,j)
   end do
end do
c-----Take out forward NP time coefficients to rotate
      do i=1,np
        do j=1,n
  recof(i,j)=ecof(i,j)
        end do
end do
c-----compute the common degree
      do i=1,m
   h(i)=0.0
        do j=1,np
   h(i)=h(i)+regvt(i,j)**2
        end do
end do
      do i=1,m
   h(i)=sqrt(h(i))
      end do
c-----compute the variance of the variance contribution by per common degree
c-----(equal to the ratated one)
      vrlv0=0.0
      do k=1,np
   g1=0.0
        g2=0.0
        do i=1,m
  g1=g1+(regvt(i,k)**2/h(i)**2)**2/real(m)
  g2=g2+(regvt(i,k)**2/h(i)**2)/real(m)
        end do
   g2=g2**2
   vrlv0=vrlv0+g1-g2
      end do
c-----rotated
      lll=0
write(*,*)'rotate times'
800  continue
      do 505 t=1,np
      do 505 k=1,np
        if(t.ge.k) goto 505
  call rot(regvt,recof,h,t,k,m,n,np)
505  continue
      lll=lll+1
write(*,*)'LLL=',LLL
      vrlv(lll)=0.0
      do k=1,np
   g1=0.0
   g2=0.0
        do i=1,m
     g1=g1+(regvt(i,k)**2/h(i)**2)**2/real(m)
      g2=g2+(regvt(i,k)**2/h(i)**2)/real(m)
        end do
   g2=g2**2
   vrlv(lll)=vrlv(lll)+g1-g2
      end do
      if(lll.lt.2) goto 800
ci=(vrlv(lll)-vrlv(lll-1))/vrlv0
if(ci.lt.0.001) goto 600
if(lll.ge.40) goto 600
      goto 800
600  continue
c-----compute the explained variances
do j=1,np
   st(j)=0.0
        do i=1,m
   st(j)=st(j)+regvt(i,j)**2
        end do
      end do
      do j=1,np
   rer(j,1)=st(j)*100/tr
end do
c-----compute the accumulated explained variances
      ddd=0.0
      do k=1,np
        ddd=ddd+st(k)/tr
        rer(k,2)=ddd*100.
      enddo
      return
      end
c-----
      subroutine rot(a,b,h,t,k,m,n,np)
c-----input arrays
      dimension a(m,np),b(np,n),h(m)
c-----work arrays
      dimension u(m),vr(m),wt(m),wk(m),wbt(n),wbk(n)
c-----output arrays
c      dimension a(m,np),b(np,n)
      integer t,k
c-----compute fi
      do 25 i=1,m
      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,m
      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/(m*1.)
      tg2=c-(s**2-g**2)/(m*1.)
      fi=atan2(tg1,tg2)/4.0
c-----compute new a and b with fi
      do 75 i=1,m
      wt(i)=a(i,t)
75   wk(i)=a(i,k)
      do 84 kk=1,n
      wbt(kk)=b(t,kk)
84   wbk(kk)=b(k,kk)
      do 78 i=1,m
      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,n
      b(t,it)=wbt(it)*cos(fi)+wbk(it)*sin(fi)
89   b(k,it)=wbt(it)*(-sin(fi))+wbk(it)*cos(fi)
      return
      end
c-----
      subroutine arrange2(m,n,np,regvt,recof,rer)
c-----input and output arrays  
      dimension regvt(m,np),recof(np,n),rer(np,2)
      mnl1=np-1
      do 20 k1=mnl1,1,-1
      do 20 k2=k1,mnl1
        if(rer(k2,1).lt.rer(k2+1,1)) then
          c=rer(k2+1,1)
          rer(k2+1,1)=rer(k2,1)
          rer(k2,1)=c
          do i=1,m
            d=regvt(i,k2+1)
            regvt(i,k2+1)=regvt(i,k2)
            regvt(i,k2)=d
          end do
          do i=1,n
            e=recof(k2+1,i)
            recof(k2+1,i)=recof(k2,i)
            recof(k2,i)=e
          end do
        endif
  20  continue
      do i=1,np
   rer(i,2)=0.0
end do
      rer(1,2)=rer(1,1)
      do 30 i=2,np
        rer(i,2)=rer(i-1,2)+rer(i,1)
  30  continue
      return
      end

      subroutine meanvar(n,x,ax,sx,vx)
!-----input array
      dimension x(n)
      ax=0.
      vx=0.
      sx=0.
      do 10 i=1,n
        ax=ax+x(i)
  10  continue
      ax=ax/float(n)
      do 20 i=1,n
        vx=vx+(x(i)-ax)**2
  20  continue
      vx=vx/float(n)
      sx=sqrt(vx)
      return
      end
37个站点51年的资料做EOF分解，程序是老师给的，为什么会出错？但是把年份改成小于或者37年就不会出错，好痛苦。。。

river

错误提示数组溢出。看看你读取的顺序和你自己的数据排列是不是一致的

吴阿好大人

river 发表于 2015-4-15 17:31
错误提示数组溢出。看看你读取的顺序和你自己的数据排列是不是一致的

是一致的，就是这个程序最多只能处理37年数据，不知道为什么

river

吴阿好大人 发表于 2015-4-15 18:04
是一致的，就是这个程序最多只能处理37年数据，不知道为什么

那就不一定是程序的问题了，你在win7或者win8上运行CVF6.5/6.6可能就会出现各种问题。建议你以XP的兼容模式运行fortran再试试

默默dreamer

楼主的问题最后解决了么？我也碰到这个问题了

lqouc

默默dreamer 发表于 2015-4-17 10:03
楼主的问题最后解决了么？我也碰到这个问题了

parameter(m=37,n=51,mnl=n,np=20,ll=m)这一句写的有问题
根据我的经验，其中mnl=min(m,n)，以楼主的情况，应该是mnl=m，这里涉及一些矩阵的运算，一般的程序会对参量进行注释说明的，用之前最好还是仔细读一读。
你可以试一试这么改。

lqouc

楼主试一下楼上给出的方法。

你、是你的

楼主用的那个6.6版的fortran吗？我最近写论文也是用eof方法，正在纠结中

centaurus

你这格点数m小于样本量n了，如果程序没改，里面的时空变换那部就不对了，所以报错。要么改雅克比方法那的程序，要么时间取小于m的值

陨落の殇

把mnl改成20以下，mnl指代特征向量输出数，20个出图够用了