登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
parameter(m=180,n=89,nyr=8,num=55,undef=32767)
real x(m,n,num),aa(m,n),x1(m,n,nyr)
real s(m,n),ax(m,n),sx(m,n),t(m,n)
integer na(nyr)
! data na/19,37,42,43,45,47,52,53,54,55/
data na/3,6,10,11,12,17,21,34/
open(1,file='ssts.grd',form='binary')
read(1)(((x(i,j,k),i=1,m),j=1,n),k=1,num)
do k=1,nyr
do j=1,n
do i=1,m
if(x(i,j,k)/=undef)then
x1(i,j,k)=x(i,j,na(k))
else
x1(i,j,k)=undef
endif
enddo
enddo
enddo
call jps(x1,aa,s,undef,nyr,m,n)
call jps(x,ax,sx,undef,num,m,n)
do i=1,m
do j=1,n
if(aa(i,j)/=undef)then
t(i,j)=(aa(i,j)-ax(i,j))/sqrt(s(i,j))*sqrt(real(nyr))
else
t(i,j)=undef
endif
enddo
enddo
open(2,file='qiang.grd',form='binary')
write(2)((aa(i,j),i=1,m),j=1,n)
write(2)((t(i,j),i=1,m),j=1,n)
end
subroutine jps(x1,aa,s,undef,num,m,n)
real p(m,n,num),s(m,n),ave(m,n),aa(m,n),x1(m,n,num)
do k=1,num
do j=1,n
do i=1,m
if(x1(i,j,k)/=undef)then
ave(i,j)=ave(i,j)+x1(i,j,k)/num
else
ave(i,j)=undef
endif
enddo
enddo
enddo
do k=1,num
do j=1,n
do i=1,m
if(x1(i,j,k)/=undef)then
p(i,j,k)=x1(i,j,k)-ave(i,j)
else
p(i,j,k)=undef
endif
enddo
enddo
enddo
do i=1,m
do j=1,n
do k=1,num
if(x1(i,j,k)/=undef)then
aa(i,j)=aa(i,j)+p(i,j,k)/num
s(i,j)=s(i,j)+p(i,j,k)**2/num
else
aa(i,j)=undef
s(i,j)=undef
endif
enddo
enddo
enddo
end
|