短期气候实习报告六之欧阳歌谷创编.docx
《短期气候实习报告六之欧阳歌谷创编.docx》由会员分享,可在线阅读,更多相关《短期气候实习报告六之欧阳歌谷创编.docx(10页珍藏版)》请在冰豆网上搜索。
短期气候实习报告六之欧阳歌谷创编
南京信息工程大学短期气候实验(实习)报告
欧阳歌谷(2021.02.01)
实验名称夏季区域降水的定量预测日期2017.6.6得分指导教师
系大气科学专业大气科学年级2014级班次2姓名车楚玉学号20141301043
一、目的要求:
目的:
掌握短期气候预测中物理统计预测的基本步骤。
要求:
能运用提供的资料和方法子程序,编写或补充完成程序当中的部分片断,了解区域降水的预测方法及其建立过程,输出实验要求的相应结果,并就方法对区域降水的拟合及试验预测效果进行分析。
二、资料和方法:
资料:
1、前期1月的Nino3.4指数(来自CPC);
2、西太平洋副高脊线、西太平洋副高西伸脊点、亚洲极涡面积、南方涛动指数(来自中国气象局整编的74个环流指数);
3、夏季华北区域10站的降水量,距平百分率。
方法:
回归分析(mregrssion.for)是用来寻找若干变量之间统计关系的一种方法,利用所找到的统计关系对某一变量作出未来时刻的估计,称为回归预报值。
三、实习步骤:
①说明所用资料方法;
②计算方法的简单介绍;
③输出反映回归效果的参数及回归系数,并就相关参数分析回归效果;
④预测量与回归方程计算的估计值和观测值的历年曲线变化图(1952~2001年),并附简单的说明;
⑤输出独立预测试验的观测与预测值。
四、结果(图表并解释说明):
程序如下:
PROGRAMMAIN
INTEGER,PARAMETER:
:
N=50
INTEGER,PARAMETER:
:
K=5
REAL,DIMENSION(K,N):
:
X
REAL,DIMENSION(N):
:
Y
REAL,DIMENSION(K+1):
:
A
REAL,DIMENSION(K+1,K+1):
:
B
REAL,DIMENSION(K):
:
V
REALQ,S,R,U,ind(6,60),year(50),y1(50),pre(7)
integeri,j
Cx(5,50)y(50)a(6)b(6,6),v(5)
COPENTHEINPUTDATAFILE
OPEN(10,FILE='D:
\duanqi\shixi6\shixi.txt')
CREADTHEDATAandgivedatatoXandY
doi=1,50
read(10,*)year(i),x(1:
5,i),y(i)
enddo
MM=K+1
callDYHG(X,Y,K,MM,N,A,Q,S,R,V,U,B,DYY)
write(*,88)A
(1)
88format(/1x,'b0=',f19.5)
do89j=2,MM
89write(*,100)j-1,A(j)
100format(1x,'b',i2,'=',f9.5)
open(11,FILE='D:
\duanqi\shixi6\guji.txt')
open(12,FILE='D:
\duanqi\shixi6\guji.grd',form='binary')
open(13,FILE='D:
\duanqi\shixi6\prediction.txt')
doj=1,50
y1(j)=a
(1)+a
(2)*x(1,j)+a(3)*x(2,j)+a(4)*x(3,j)+a(5)*x(4,j)+a(6)
*x(5,j)
enddo
pre
(1)=a
(1)+a
(2)*26.50+a(3)*13.00+a(4)*100.00+a(5)*196.00+a(6)*
2.00
pre
(2)=a
(1)+a
(2)*27.76+a(3)*12.00+a(4)*90.0+a(5)*197.00+a(6)*(-1.00)
pre(3)=a
(1)+a
(2)*26.74+a(3)*12.0+a(4)*120.00+a(5)*199.00+a(6)*
(-11.00)
pre(4)=a
(1)+a
(2)*27.10+a(3)*13.00+a(4)*110.00+a(5)*220.00+a(6)*3.00
pre(5)=a
(1)+a
(2)*25.64+a(3)*14.00+a(4)*105.00+a(5)*215.00+a(6)*13.00
pre(6)=a
(1)+a
(2)*27.26+a(3)*15.00+a(4)*90.00+a(5)*185.00+a(6)*
(-8.00)
pre(7)=a
(1)+a
(2)*24.71+a(3)*17.00+a(4)*125.00+a(5)*222.00+a(6)*13.00
print*,pre
write(11,*)(y(i),y1(i),i=1,50)
write(12)(y(i),y1(i),i=1,50)
write(13,*)(pre(i),i=1,7)
cccccccccccccccccccccccccccccccccccccccccccccccccccccc
write(*,20)Q,S,R
20format(1x,'Q=',f13.6,3x,'S=',f13.6,3x,'R=',f13.6)
write(*,22)U,DYY
22format(1x,'U=',f13.6,3x,'DYY=',f13.6)
write(*,30)(i,V(i),i=1,K)
30format(1x,'V(',i2,')=',f13.6)
write(*,40)U
40format(1x,'U=',f13.6)
open(6,file='table')!
outputdata
write(6,180)
180format(/2x,'regressioncoefficients:
')
write(6,88)A
(1)
do189j=2,MM
189write(6,100)j-1,A(j)
write(6,200)
200format(/1x,'GenericAnalysisofVarianceTablefortheMultiple
LinearRegression')
write(6,202)
202format(/1x,'-----------------------------------------------------
---------------')
write(6,204)
204format(/3x,'SourcedfSSMS')
write(6,202)
write(6,206)N-1,DYY
206format(/1x,'Totaln-1=',i2,'SST=',f13.4)
u2=U/real(K)
write(6,208)K,U,U2
208format(/1x,'RegressionK=',i2,'SSR=',f13.4,'MSR=SSR/K='
f13.4)
q2=q/real(n-k-1)
write(6,209)n-k-1,q,q2
209format(/1x,'Residualn-k-1=',i2,'SSE=',f13.4,'MSE=SSE/(n-k-1)
=',f13.4)
f=(U/real(K))/(Q/real(N-K-1))
write(6,220)f
220format(/1x,'F=MSR/MSE=',f13.4)
write(6,202)
close(6)
stop
end
subroutineDYHG(x,y,m,mm,n,a,q,s,r,v,u,b,dyy)
dimensionx(m,n),y(n),a(mm),b(mm,mm),v(m)
b(1,1)=n
do20j=2,mm
b(1,j)=0.0
do10i=1,n
10b(1,j)=b(1,j)+x(j-1,i)
b(j,1)=b(1,j)
20continue
do50i=2,mm
do40j=i,mm
b(i,j)=0.0
do30k=1,n
30b(i,j)=b(i,j)+x(i-1,k)*x(j-1,k)
b(j,i)=b(i,j)
40continue
50continue
a
(1)=0.0
do60i=1,n
60a
(1)=a
(1)+y(i)
do80i=2,mm
a(i)=0.0
do70j=1,n
70a(I)=a(i)+x(i-1,j)*y(j)
80continue
callCHOLESKY(b,mm,1,a,l)
yy=0.0
do90i=1,n
90yy=yy+y(i)/n
q=0.0
dyy=0.0
u=0.0
cccccccccccccccccccccccccccccccccc
do110i=1,n
p=a
(1)
do100j=1,m
100p=p+a(j+1)*x(j,i)
q=q+(y(i)-p)*(y(i)-p)
dyy=dyy+(y(i)-yy)*(y(i)-yy)
u=u+(yy-p)*(yy-p)
110continue
cccccccccccccccccccccccccccccccccc
s=sqrt(q/n)
r=sqrt(1.0-q/dyy)
do150j=1,m
p=0.0
do140i=1,n
pp=a
(1)
do130k=1,m
if(k.ne.j)pp=pp+a(k+1)*x(k,i)
130continue
p=p+(y(i)-pp)*(y(i)-pp)
140continue
v(j)=sqrt(1.0-q/p)
150continue
return
end
subroutineCHOLESKY(a,n,m,d,l)!
PerformtheCHOLESKYDecomposition
dimensiona(n,n),d(n,m)
l=1
if(a(1,1)+1.0.eq.1.0)then
l=0
write(*,30)
return
endif
a(1,1)=sqrt(a(1,1))
do10j=2,n
10a(1,j)=a(1,j)/a(1,1)
do100i=2,n
do20j=2,i
20a(i,i)=a(i,i)-a(j-1,i)*a(j-1,i)
if(a(i,i)+1.0.eq.1.0)then
l=0
write(*,30)
return
endif
30format(1x,'fail')
a(i,i)=sqrt(a(i,i))
if(i.ne.n)then
do50j=I+1,n
do40k=2,i
40a(i,j)=a(i,j)-a(k-1,i)*a(k-1,j)
50a(i,j)=a(i,j)/a(I,i)
endif
100continue
do130j=1,m
d(1,j)=d(1,j)/a(1,1)
do120i=2,n
do110k=2,i
110d(i,j)=d(i,j)-a(k-1,i)*d(k-1,j)
d(i,j)=d(i,j)/a(i,i)
120continue
130continue
do160j=1,m
d(n,j)=d(n,j)/a(n,n)
do150k=n,2,-1
do140i=k,n
140d(k-1,j)=d(k-1,j)-a(k-1,i)*d(i,j)
d(k-1,j)=d(k-1,j)/a(k-1,k-1)
150continue
160continue
return
end
表一1952-2001年回归方程计算的估计值和观测值的历年曲线变化图
2002200320042005200620072008
观测值-29.75-17.450.64-10.63-22.96-19.403.30
预测值30.8053339.2322031.9930128.4085024.2936114.32182-2.514375
表二2002-2008年预测试验的观测与预测值
五、结果讨论:
分析表一,白色实线为1952-2001年的观测真实值,绿色实线为1952-2001年的回归方程计算估计值,由表可以看出,某些年份,预测值与观测值的趋势大致相同,而在某些年份趋势相反,如1959年为极大值,但是估计值却为极小值,且各个年份的两个值的差别也很大,可以看出,通过回归方程预测的值不是很理想。
通过表二可以看出,预报观测值与回归方程的预测值差别十分显著,进一步表明了多元线形回归方程的预测值并不理想。