- 积分
- 5736
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2015-4-19
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
我需要计算空间线性趋势分布,因为数据存在缺测值,我在Fortran处理的过程中遇到一些问题,希望各位大神可以帮我解决怎么处理缺测的数值,感激不尽 ~另外,帮我看看程序是否也有问题啊
!x=bt+a
!四季发生频次线性趋势(春夏秋冬季频次累加)
program station
parameter (nx=128,ny=72,n=51,nk=4,nn=624,undef=-999*3)
real h(nx,ny,nn),t(n),a1(nx,ny,n,nk),b1(nx,ny,nk),x(n),i,j,it,mon
s=0;s1=0;s2=0;s3=0;s4=0;a=0;b=0;num=0
open(1,file='F:\data\new\china-monthly-tx90p-all.grid.grd',form='binary')
open(2,file='F:\data\new\china-monthly-tx90p-all-b.grd',form='binary')
do j=1,ny
do i=1,nx
do mon=1,nn
read(1) (h(i,j,mon))
enddo;enddo;enddo
!4季
do i=1,nx
do j=1,ny
do it=1,n
do k=1,nk
a1(i,j,it,k)=h(i,j,(it-1)*12+(k-1)*3+3)+h(i,j,(it-1)*12+(k-1)*3+4)+h(i,j,(it-1)*12+(k-1)*3+5)
enddo;enddo;enddo;enddo
do i=1,nx
do j=1,ny
do k=1,nk
num=0
do it=1,n
if(a1(i,j,it,k).ne.undef)then !这里没有很好的处理缺测数据,单纯的没用undef的频次来计算,忽略了其中有缺测和没缺测月份的累加
num=num+1
x(num)=a1(i,j,it,k)
t(num)=num
else
b1(i,j,k)=-999
endif
enddo
call SLR(x,t,num,a,b)
b1(i,j,k)=b
enddo
enddo
enddo
WRITE(2) b1
end
SUBROUTINE SLR(x,t,num,a,b)
real x(num),t(num),s1,s2,s3,s4,a,b
s1=0;s2=0;s3=0;s4=0
do k=1,num
s1=s1+t(k)*x(k)
s2=s2+t(k)
s3=s3+x(k)
s4=s4+t(k)**2
enddo
b=(s1*num-s2*s3)/(s4*num-s2**2)
a=(s3-b*s2)/num
end subroutine
|
|