短期气候实习报告六之欧阳地创编.docx
《短期气候实习报告六之欧阳地创编.docx》由会员分享，可在线阅读，更多相关《短期气候实习报告六之欧阳地创编.docx（10页珍藏版）》请在冰豆网上搜索。
短期气候实习报告六之欧阳地创编
南京信息工程大学短期气候实验（实习）报告
时间：
2021.03.04
创作：
欧阳地
实验名称夏季区域降水的定量预测日期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年为极大值，但是估计值却为极小值，且各个年份的两个值的差别也很大，可以看出，通过回归方程预测的值不是很理想。
通过表二可以看出，预报观测值与回归方程的预测值差别十分显著，进一步表明了多元线形回归方程的预测值并不理想。
时间：
2021.03.04
创作：
欧阳地
升级会员 