- 积分
- 223
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2011-9-25
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
求助~~好奇怪 我计算指数与场的相关 如果用原始场没有问题 可是把场进行谐波分析后,再用如下程序运行发现好多格点方差变为零啦~~好不理解,然后相关系数就变成NA啦~~~求大神解答~~~!!!!
parameter (nx=180,ny=89,m=54)
real sst(nx,ny,m),in(m),rxy(nx,ny),t(nx,ny)
do k=1,m
read(1)((sst(i,j,k),i=1,nx),j=1,ny)
enddo
open(2,file='inslow.grd',form='binary')
read(2)(in(i),i=1,m)
enddo
call xiangguan(sst,in,rxy,t)
open(3,file='slow-index.grd',form='binary')
open(4,file='slow-t.grd',form='binary')
write(3)((rxy(i,j),i=1,nx),j=1,ny)
write(4)((t(i,j),i=1,nx),j=1,ny)
end
subroutine xiangguan(a,b,r,t)
parameter (nx=180,ny=89,m=54,undef=32767)
real a(nx,ny,m),b(m),sa(nx,ny),sb,sab(nx,ny),r(nx,ny),avea(nx,ny),aveb,t(nx,ny)
do k=1,m
aveb=aveb+b(k)/m
enddo
do j=1,ny
do i=1,nx
do k=1,m
if(a(i,j,k)/=undef)then
avea(i,j)=avea(i,j)+a(i,j,k)/real(m)
else
avea(i,j)=undef
endif
enddo
enddo
enddo
do k=1,m
do j=1,ny
do i=1,nx
if(a(i,j,1)/=undef)then
sa(i,j)=sa(i,j)+(a(i,j,k)-avea(i,j))**2/real(m)
sab(i,j)=sab(i,j)+(a(i,j,k)-avea(i,j))*(b(k)-aveb)/real(m)
else
sa(i,j)=undef
sab(i,j)=undef
endif
enddo
enddo
sb=sb+(b(k)-aveb)**2/real(m)
enddo
open(100,file='test.txt')
write(100,*)((sa(i,j),i=1,nx),j=1,ny)
do j=1,ny
do i=1,nx
if(sa(i,j)/=undef.and.sa(i,j)/=0)then
r(i,j)=sab(i,j)/sqrt(sa(i,j))/sqrt(sb)
else
r(i,j)=undef
endif
!if(abs(r(i,j))>=1)then
!write(*,*)r(i,j)
!endif
enddo
enddo
do j=1,ny
do i=1,nx
if(r(i,j)/=undef)then
!if(abs(r(i,j))<=1)then
t(i,j)=r(i,j)/sqrt(1-r(i,j)**2)*sqrt(real(m-2))
!endif
else
t(i,j)=undef
endif
enddo
enddo
end
|
|