线性趋势

hfh-grads

我需要计算空间线性趋势分布，因为数据存在缺测值，我在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

言深深

初步看下来，程序的可读性还是比较ok的。
1、季节计算的似乎不对，主要是冬季的计算，你看看是不是将当年的12,1,2加起来了？应该是前一年12月，当年1,2月；
2、如果数据某一个月有缺测，那么在计算季节平均就错了，因此判断是否缺测那一步就判断不出来了。

hfh-grads

言深深 发表于 2015-5-5 18:30
初步看下来，程序的可读性还是比较ok的。
1、季节计算的似乎不对，主要是冬季的计算，你看看是不是将当年 ...

谢谢啊~我是把当年的12月份以及次年的1,2月份定义为冬季的

言深深

hfh-grads 发表于 2015-5-5 19:23
谢谢啊~我是把当年的12月份以及次年的1,2月份定义为冬季的

好吧···
但是问题2，似乎仍然存在。请核实

呼啦呼啦

言深深 发表于 2015-5-6 18:16
好吧···
但是问题2，似乎仍然存在。请核实

深深你好，关于这个问题我想问一下，如果我不考虑缺测和季节的话，想做线性趋势的空间分布，我自己照着上面的这个程序编了一个，能不能帮我看看有没有什么问题，非常感谢！我的资料是一个逐月的资料，32年，384个月
!x=bt+a
program station
parameter (nx=105,ny=14,mt=384)
real tmsu(nx,ny,mt),t(mt),x(mt),b1(nx,ny)
integer i,j,it,k

open(1,file='D:\by\evalution_MSU_AMSU_BT\evalution\MSU_AMSU\AMSU\NORM\Trend_Space_distribution\1_DATA\TLS_JP_lev100.grd',form='binary')
do it=1,mt
read(1) ((tmsu(i,j,it),i=1,nx),j=1,ny)
end do


do i=1,nx
do j=1,ny
num=0
do it=1,mt
num=num+1
x(num)=tmsu(i,j,it)
t(num)=num
enddo
call SLR(x,t,num,b)
b1(i,j)=b
enddo
enddo

open(2,file='D:\by\evalution_MSU_AMSU_BT\evalution\MSU_AMSU\AMSU\NORM\Trend_Space_distribution\b_space_distribution\b_grid.grd',form='binary')
write(2)((b1(i,j),i=1,nx),j=1,ny)
END

SUBROUTINE SLR(x,t,num,b)
real x(num),t(num),s1,s2,s3,s4,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)
end subroutine

言深深

呼啦呼啦 发表于 2015-5-29 14:30
深深你好，关于这个问题我想问一下，如果我不考虑缺测和季节的话，想做线性趋势的空间分布，我自己照着上 ...

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)
你这个公式似乎不对。

呼啦呼啦

言深深 发表于 2015-5-29 16:19
do k=1,num
s1=s1+t(k)*x(k)
s2=s2+t(k)

恩，好像是有问题，后来我又修改了一下，能帮我再看一下吗？
!x=bt+a
program station
parameter (nx=105,ny=13,mt=384)
real tmsu(nx,ny,mt),yr(mt)
real s1(nx,ny),s2(nx,ny),s3(nx,ny),s4(nx,ny),b(nx,ny)
integer i,j,it,k

open(1,file='D:\by\evalution_MSU_AMSU_BT\evalution\MSU_AMSU\AMSU\NORM\annual_trend\Total_ocean\1Cut_DATA\NCEP2_1981_2012_lev100.grd',form='binary')
do it=1,mt
read(1) ((tmsu(i,j,it),i=1,nx),j=1,ny)
end do

do it=1,mt
yr(it)=it
end do

do i=1,nx
do j=1,ny
do it=1,mt
s1(i,j)=s1(i,j)+tmsu(i,j,it)*yr(it)
s2(i,j)=s2(i,j)+yr(it)
s3(i,j)=s3(i,j)+tmsu(i,j,it)
s4(i,j)=s4(i,j)+yr(it)*yr(it)
enddo
b(i,j)=(s1(i,j)-(s2(i,j)*s3(i,j))/384.0)/(s4(i,j)-(s2(i,j)**2)/384.0)
enddo
enddo

open(2,file='D:\by\evalution_MSU_AMSU_BT\evalution\MSU_AMSU\AMSU\NORM\Trend_Space_distribution\b_NCEP2_1981_2012_lev100.grd',form='binary')
write(2)((b(i,j)*10,i=1,nx),j=1,ny)
END

言深深

呼啦呼啦 发表于 2015-6-1 09:25
恩，好像是有问题，后来我又修改了一下，能帮我再看一下吗？
!x=bt+a
program station

计算b(i,j)仍然是错的，公式如下：
a=(NΣxy-ΣxΣy)/(NΣx^2-(Σx)^2)

应该写为：

b(i,j)=(s1(i,j)*mt-s2(i,j)*s3(i,j))/(s4(i,j)*mt-s2(i,j)**2)

呼啦呼啦

言深深 发表于 2015-6-1 09:48
计算b(i,j)仍然是错的，公式如下：
a=(NΣxy-ΣxΣy)/(NΣx^2-(Σx)^2)

可是魏凤英老师的书里面P37上写的公式是：
(Σxt-(ΣxΣt)/n)/(Σt^2-((Σt)^2)/n)
我做的就是一个时间序列数据的线性趋势，不知道是什么情况呀

言深深

呼啦呼啦 发表于 2015-6-1 10:04
可是魏凤英老师的书里面P37上写的公式是：
(Σxt-(ΣxΣt)/n)/(Σt^2-((Σt)^2)/n)
我做的就是一个时间 ...

额，是的哎。是我看错了。
你试试程序看看，这两个都是一样的。