短期气候实习报告六

1、南京信息工程大学短期气候实验（实习）报告实验名称 夏季区域降水的定量预测 日期2017.6.6得分指导教师 系大气科学专业大气科学年级2014级班次2姓名车楚玉学号20141301043一、目的要求：目的：掌握短期气候预测中物理统计预测的基本步骤。要求：能运用提供的资料和方法子程序，编写或补充完成程序当中的部分片断，了解区域降水的预测方法及其建立过程，输出实验要求的相应结果，并就方法对区域降水的拟合及试验预测效果进行分析。二、资料和方法：资料：1、前期1月的nino3.4指数（来自cpc）;2、西太平洋副高脊线、西太平洋副高西伸脊点、亚洲极涡面积、南方涛动指数（来自中国气象局整编的 74个环流

2、指数）；3、夏季华北区域10站的降水量，距平百分率。方法：回归分析（mregrssion.for）是用来寻找若干变量之间统计关系的一种方法，利用所找 到的统计关系对某一变量作出未来时刻的估计，称为回归预报值。三、实习步骤：说明所用资料方法； 计算方法的简单介绍； 输出反映回归效果的参数及回归系数，并就相关参数分析回归效果；预测量与回归方程计算的估计值和观测值的历年曲线变化图（19522001年），并附简单的说明；输出独立预测试验的观测与预测值。四、结果（图表并解释说明）程序如下：programainintegerparameter=50integerparameter=5 realdimens

3、io(k,n):x realdimensionj):y realdimensionk+1):a realdimensionk+1,k+1):b realdimensionk):vrealq,s,r,u,ind(6,60),year(50),y1(50),pre(7) integer i,jc x(5,50) y(50) a(6) b(6,6),v(5)c open the input data fileopen0,file= 'd:duanqishixi6shixi.txt')c read the data and give data to x and y do i=1,50

4、|read (10,pre(1)=a(1)+a(2)*26.50+a(3)*13.00+a(4)*100.00+a(5)*196.00+a(6)* 180 format (/2x, 'regression coefficients:' )write (6,88) a(1)do 189 j=2,mm189 write (6,100) j-1,a(j)write (6,200)200 format (/1x, 'generic analysis of variance table for the multiple) year(i),x(1:5,i),y(i)enddomm=

5、k+1call dyhg(x,y,k,mm,n,a,q,s,r,v,u,b,dyy) write (*,88) a(1)88 format (/1x, 'b 0=' ,f19.5) do 89 j=2,mm89write (*,100) j-1,a(j)100format (1x, 'b' ,i2, '=' ,f9.5)open(11,file= 'd:duanqishixi6guji.txt' ) open(12,file= 'd:duanqishixi6guji.grd' ,form= 'binary&

6、#39;) open(13,file= 'd:duanqishixi6prediction.txt')do j=1,50y1(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)enddopre(4)=a(1)+a(2)do 140 i=k,n140d(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)150continue160continuereturnend27.10+a(3)*13.00+a(4)*110.00+a(5

7、)*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 *,prewrite (11,*)(y(i),y1(i),i=1,50)write (12)(y(i),y1(i),i=1,50)wri

8、te (13,*)(pre(i),i=1,7)ccccccccccccccccccccccccccccccccccccccccccccccccccccccwrite (*,20)q,s,r20 format (1x, 'q=' ,f13.6,3x,'s=' ,f13.6,3x,'r=' ,f13.6)write (*,22)u,dyy22 format (1x, 'u=' ,f13.6,3x, 'dyy=' ,f13.6) write (*,30)(i,v(i),i=1,k)30 format (1x, '

9、v(' ,i2, ')=' ,f13.6) write (*,40)u40 format (1x, 'u=' ,f13.6)open(6,file= 'table' )! output datawrite (6,180)f=(u/ real (k)/(q/ real (n-k-1) write (6,220) f220 format (/1x, 'f=msr/mse=',f13.4)write (6,202) close (6) stop endsubroutine dyhg(x,y,m,mm,n,a,q,s,r,v,u,

10、b,dyy) dimension x(m,n),y(n),a(mm),b(mm,mm),v(m) b(1,1)=ndo 20 j=2,mm b(1,j)=0.0 do 10 i=1,n 10b(1,j)=b(1,j)+x(j-1,i)b(j,1)=b(1,j)20 continue do 50 i=2,mm do 40 j=i,mm b(i,j)=0.0do 100 j=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)110 continuecccccccc

11、cccccccccccccccccccccccccc s= sqrt (q/n) r= sqrt (1.0-q/dyy) | do 150 j=1,m |p=0.0do 140 i=1,n pp=a(1)do 130 k=1,mj if (k.ne.j)pp=pp+a(k+1)*x(k,i)130 continuep=p+(y(i)-pp)*(y(i)-pp)140continuev(j)= sqrt (1.0-q/p)endif30 format (1x, 'fail') a(i,i sqrt (a(i,i) if (i.ne.n) then do 50 j=i+1,n do

12、 40 k=2,i40a(i,j)=a(i,j)-a(k-1,i)*a(k-1,j)50a(i,j)=a(i,j)/a(i,i)endif100 continue do 130 j=1,m. d(1,j尸d(1,j)/a(1,1) do 120 i=2,n do 110 k=2,i110d(i,j)=d(i,j)-a(k-1,i)*d(k-1,j)d(i,j)=d(i,j)/a(i,i)120continue130continuedo 160 j=1,m d(n,j)=d(n,j)/a(n,n)do 150 k=n,2,-1表一 1952-2001年回归方程计算的估计值和观测值的历年曲线变化图

13、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年为极大值，但是估计值却为极小值，且各个年份的两个值的差别也很大，可以看出，通过回归方

14、程预测的值不是很理想。通过表二可以看出，预报观测值与回归方程的预测值差别十分显著，进一步表明了多元线形回归方程的预测值并不理想。id* 2.00pre(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)* linear regression' )write (6,202)202 format (/1x,'*')write (6,204)204 for

15、mat (/3x, 'source df ss ms')write (6,202)write (6,206) n-1,dyy206 format (/1x, 'total n-1=',i2, ' sst=' ,f13.4)u2=u/ real (k)write (6,208) k,u,u2208 format (/1x, 'regression k=',i2, ' ssr=' ,f13.4, ' msr=ssr/k='*,f13.4)q2=q/ real (n-k-1)write (6,209)n-

16、k-1,q,q2209 format (/1x, 'residual n-k-1=',i2, ' sse=' ,f13.4, ' mse=sse/(n-k-1)*=' ,f13.4)do 30 k=1,n30b(i,j尸b(i,j)+x(i-1,k)*x(j-1,k)b(j,i)=b(i,j)40continue50 continue a(1)=0.0do 60 i=1,n60a(1)=a(1)+y(i)do 80 i=2,mm a(i)=0.0do 70 j=1,n70a(i)=a(i)+x(i-1,j)*y(j)80continuecall cholesky(b,mm,1,a,l) yy=0.0do 90 i=1,n90 yy=yy+y(i)/n q=0.0 dyy=0.0 u=0.0cccccccccccccccccccccccccccccccccc do 110 i=1,n p=a(1)150continuereturnendsubroutine cholesky(a,n,m,d,l) ! perform the chole