! EOF：经验正交函数分解! “实现时空场的分离”具体可以参考黄嘉佑老师的《气象统计分析与预报方法》
! 这是来自“lijianping”老师的EOF分解程序，在此对李建平老师表示感谢O(∩_∩)O~
! 相关参数说明
! m：格点数
! n：时间长度
! mnl：模态数
! ks：[1]-1表示对原场EOF；
! [2]0表示多距平场EOF；
! [3]1表示对归一化场EOF。
! "test.txt" ：需要分解的时空场
! "er.txt"：存放解释方差 [1]：er(mnl,1)特征值—— [2]：er(mnl,2)累计特征值—— [3]：er(mnl,3)解释方差—— [4]：er(mnl,4)累计解释方差!【20120922补充：解释方差就是对特征值的归一化】
! "ecof.txt"：存放时间系数------每一列是一个模态
! "egvt.txt"：存放各个模态的空间向量场------每一列是一个模态
integer,parameter :: nx=46
integer,parameter :: ny=46
integer,parameter :: n = 456
integer,parameter :: m= 1665
integer,parameter :: mnl= 456
real, dimension(m,n):: x
real, dimension(m,mnl):: egvt
real, dimension(mnl,n):: ecof
real, dimension(mnl,4):: er
real, dimension(m) :: lonin,latin
integer,dimension(n) :: kk
integer::k,ks=1
integer :: istatus
open(1,file="D:\method\EOF\sst.mnmean.v4_.txt",iostat=istatus)
do k=1,n
do j=1,m
read(1,'(i5,2f10.4,f10.1)') kk(k),latin(j),lonin(j), x(j,k)
end do; end do
call eof(m,n,mnl,x,ks,er,egvt,ecof)
open(2,file="D:\method\EOF\sst.mnmean.er456+1.txt")
do i=1,mnl
write(2,"(4f30.8)") (er(i,j),j=1,4)
enddo
close(2)
open(3,file="D:\method\EOF\sst.mnmean.ecof456+1.txt")
do j=1,n
write(3,"(<mnl>f20.6)") (ecof(i,j),i=1,mnl)
enddo
close(3)
open(4,file="D:\method\EOF\sst.mnmean.egvt456+1.txt")
do j=1,m
write(4,*) (latin(j),lonin(j), egvt(j,i),i=1,mnl)
end do
close(4)
end
!***************************************************************************
subroutine eof(m,n,mnl,f,ks,er,egvt,ecof)
implicit none
integer*4 :: m,n,mnl,ks
real*4 :: f(m,n),er(mnl,4),egvt(m,mnl),ecof(mnl,n)
real*4,allocatable :: cov(:,:),s(:,:),d(:)
call transf(m,n,f,ks)
allocate(cov(mnl,mnl))
call crossproduct(m,n,mnl,f,cov)
allocate(s(mnl,mnl))
allocate(d(mnl))
call jacobi(mnl,cov,s,d,0.00001)
call arrang(mnl,d,s,er)
call tcoeff(m,n,mnl,f,s,er,egvt,ecof)
return
end
subroutine transf(m,n,f,ks)
implicit none
integer*4 :: m,n,ks
real*4 :: f(m,n)
real*4,allocatable :: fw(:),wn(:)
integer*4 :: i0,i,j
real*4 :: af,sf,vf
allocate(fw(n))
allocate(wn(m))
i0=0
do i=1,m
do j=1,n
fw(j)=f(i,j)
enddo
call meanvar(n,fw,af,sf,vf)
if(sf.eq.0.)then
i0=i0+1
wn(i0)=i
endif
enddo
if(i0.ne.0)then
write(*,*)'************************* fault ***********************************'
write(*,*)'The program cannot go on because the original field has invalid data.'
write(*,*)'There are totally ',i0,'gridpionts with invalid data.'
write(*,*)'The array wn(m) stores the positions of those invalid grid-points.'
write(*,*)'You must pick off those invalid data from the orignal field---$ '
write(*,*)'$---and then reinput a new field to calculate its eofs.'
write(*,*)'************************ fault ************************************'
endif
if(ks.eq.-1)return
if(ks.eq.0)then
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
do i=1,m
do j=1,n
fw(j)=f(i,j)
enddo
call meanvar(n,fw,af,sf,vf)
if(sf==0.0)then
do j=1,n
f(i,j)=0.0
enddo
else
do j=1,n
f(i,j)=(f(i,j)-af)/sf
enddo
endif
enddo
endif
return
end
subroutine crossproduct(m,n,mnl,f,cov)
implicit none
integer*4 :: m,n,mnl
real*4 :: f(m,n),cov(mnl,mnl)
integer*4 :: i,j,is,js
if((n-m)<0)then
do i=1,mnl
do 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)
enddo
enddo
else
do i=1,mnl
do 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)
enddo
enddo
endif
return
end
subroutine jacobi(m,a,s,d,eps)
implicit none
integer*4 :: m
real*4 :: a(m,m),s(m,m),d(m),eps
integer*4 :: i,j,i1,l,iq,iq1,ip
real*4 :: g,s1,s2,s3,v1,v2,v3,u,st,ct
s=0.0
do i=1,m
s(i,i)=1.0
enddo
g=0.
do i=2,m
i1=i-1
do j=1,i1
g=g+2.0*a(i,j)*a(i,j)
enddo
enddo
s1=sqrt(g)
s2=eps/float(m)*s1
s3=s1
l=0
50 s3=s3/float(m)
60 do iq=2,m
iq1=iq-1
do ip=1,iq1
