- 积分
- 269
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2012-7-12
- 最后登录
- 1970-1-1

|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
在兰溪之水的EOF程序基础上改的矢量EOF分解,flux2.dat是放了两个变量的数据,73*17,结果算出来的egvt.2里面的值为零,请问是哪里出错了呢?
program main
character*80,parameter :: datafname='E:\ep-flux\flux2.dat'
parameter(n=30,mx=73,my=34,m=mx*my,mnl=n)
parameter(undef=-9.99e+08,mg1=mx*my,mg2=mx*my)
parameter(ks=0,km=0,nd=1,std=1e-6)
!-----input array
dimension f0(n,m)
!-----work arrays
dimension f1(n,mg1),f2(n,mg1),g(mg1),h1(mg1,n),h2(mg1,n)
dimension f(mg2,n),gvt(mg2,mnl),cof(mnl,n)
!-----output arrays
dimension er(mnl,4),egvt(m,mnl),ecof(mnl,n)
c-----Read data.
open(20,file=datafname,form='unformatted',
& access='direct',recl=m)
do j=1,n
read(20,rec=j)(f0(j,i),i=1,m)
end do
close(20)
write(*,*)'Read data OK!'
c-----Remove the terrain or missing value.
l1=0
do j=1,m
if(f0(1,j).ne.undef)then
l1=l1+1
do k=1,n
f1(k,l1)=f0(k,j)
enddo
endif
enddo
write(*,*)'Grids without terrain:'
write(*,*)'mg1=',l1
c-----Remove annual cycle.
if(km.eq.1)then
ny=n/nd
do i=1,l1
call initial(nd,ny,n,f1(1,i),f2(1,i))
end do
do i=1,l1
do k=1,n
f1(k,i)=f2(k,i)
end do
end do
end if
c-----Remove the grids whose variance equal to zero.
l2=0
do i=1,l1
call meanvar(n,f1(1,i),ax,g(i),vx)
if(g(i).gt.std)then
l2=l2+1
do k=1,n
f(l2,k)=f1(k,i)
enddo
endif
end do
write(*,*)'Grids without terrain and constant value:'
write(*,*)'mg2=',l2
c stop
c-----Call the subroutine.
write(*,*)'!!!!'
call eof(l2,n,mnl,f,ks,er,gvt,ecof)
write(*,*)'EOF ok and transform to the original form in the next!'
c-----Add the grids whose variance equal to zero.
l3=0
do i=1,mg2
if(g(i).gt.std)then
l3=l3+1
do k=1,mnl
h1(i,k)=gvt(l3,k)
end do
else
do k=1,mnl
h1(i,k)=undef
end do
endif
enddo
c-----Add the terrain or missing value.
l4=0
do i=1,m
if(f0(1,i).ne.undef)then
l4=l4+1
do k=1,mnl
egvt(i,k)=h1(l4,k)
end do
else
do k=1,mnl
egvt(i,k)=undef
end do
endif
enddo
c-----output the result.
c-----output eigenvactors matrix of EOF.
open(11,file='egvt1.dat',status='unknown'
&,form='unformatted',access='direct',recl=mx*my)
do j=1,1
write(11,rec=j)(egvt(i,j),i=1,m/2)
end do
write(*,*)' OK: output eigenvactors matrix of EOF!'
close(11)
open(1111,file='egvt2.dat',status='unknown'
&,form='unformatted',access='direct',recl=mx*my)
do j=1,1
write(1111,rec=j)(egvt(i,j),i=m/2+1,m/2)
end do
write(*,*)' OK: output eigenvactors matrix of EOF!'
close(1111)
下面是EOF的子程序,都没有什么改动,应该是上面这部分出错的,希望大家能帮忙找一下~
|
|