短期气候实习报告六

南京信息工程大学短期气候实验（实习）报告实验名称夏季区域降水的定量预测日期2017.6.6得分指导教师系大气科学专业大气科学年级2014级班次2姓名车楚玉学号20141301043一、目的要求：目的：掌握短期气候预测中物理统计预测的基本步骤。要求：能运用提供的资料和方法子程序，编写或补充完成程序当中的部分片断，了解区域降水的预测方法及其建立过程，输出实验要求的相应结果，并就方法对区域降水的拟合及试验预测效果进行分析。二、资料和方法：资料：1、前期1月的Nino3.4指数（来自CPC）；2、西太平洋副高脊线、西太平洋副高西伸脊点、亚洲极涡面积、南方涛动指数（来自中国气象局整编的74个环流指数）；3、夏季华北区域10站的降水量，距平百分率。方法：回归分析(mregrssion.for)是用来寻找若干变量之间统计关系的一种方法，利用所找到的统计关系对某一变量作出未来时刻的估计，称为回归预报值。三、实习步骤：①说明所用资料方法；②计算方法的简单介绍；③输出反映回归效果的参数及回归系数，并就相关参数分析回归效果；④预测量与回归方程计算的估计值和观测值的历年曲线变化图（1952~2001年），并附简单的说明；⑤输出独立预测试验的观测与预测值。四、结果（图表并解释说明）：程序如下：PROGRAMMAININTEGER,PARAMETER::N=50INTEGER,PARAMETER::K=5REAL,DIMENSION(K,N)::XREAL,DIMENSION(N)::YREAL,DIMENSION(K+1)::AREAL,DIMENSION(K+1,K+1)::BREAL,DIMENSION(K)::VREALQ,S,R,U,ind(6,60),year(50),y1(50),pre(7)integeri,jCx(5,50)y(50)a(6)b(6,6),v(5)COPENTHEINPUTDATAFILEOPEN(10,FILE='D:\duanqi\shixi6\shixi.txt')CREADTHEDATAandgivedatatoXandYdoi=1,50read(10,*)year(i),x(1:5,i),y(i)enddoMM=K+1callDYHG(X,Y,K,MM,N,A,Q,S,R,V,U,B,DYY)write(*,88)A(1)88format(/1x,'b0=',f19.5)do89j=2,MM89write(*,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,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(1)=a(1)+a(2)*26.50+a(3)*13.00+a(4)*100.00+a(5)*196.00+a(6)**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)pre(4)=a(1)+a(2)*27.10+a(3)*13.00+a(4)*110.00+a(5)*220.00+a(6)*3.00pre(5)=a(1)+a(2)*25.64+a(3)*14.00+a(4)*105.00+a(5)*215.00+a(6)*13.00pre(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.00print*,prewrite(11,*)(y(i),y1(i),i=1,50)write(12)(y(i),y1(i),i=1,50)write(13,*)(pre(i),i=1,7)ccccccccccccccccccccccccccccccccccccccccccccccccccccccwrite(*,20)Q,S,R20format(1x,'Q=',f13.6,3x,'S=',f13.6,3x,'R=',f13.6)write(*,22)U,DYY22format(1x,'U=',f13.6,3x,'DYY=',f13.6)write(*,30)(i,V(i),i=1,K)30format(1x,'V(',i2,')=',f13.6)write(*,40)U40format(1x,'U=',f13.6)open(6,file='table')!outputdatawrite(6,180)180format(/2x,'regressioncoefficients:')write(6,88)A(1)do189j=2,MM189write(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,DYY206format(/1x,'Totaln-1=',i2,'SST=',f13.4)u2=U/real(K)write(6,208)K,U,U2208format(/1x,'RegressionK=',i2,'SSR=',f13.4,'MSR=SSR/K='*,f13.4)q2=q/real(n-k-1)write(6,209)n-k-1,q,q2209format(/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)f220format(/1x,'F=MSR/MSE=',f13.4)write(6,202)close(6)stopendsubroutineDYHG(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)=ndo20j=2,mmb(1,j)=0.0do10i=1,n10b(1,j)=b(1,j)+x(j-1,i)b(j,1)=b(1,j)20continuedo50i=2,mmdo40j=i,mmb(i,j)=0.0do30k=1,n30b(i,j)=b(i,j)+x(i-1,k)*x(j-1,k)b(j,i)=b(i,j)40continue50continuea(1)=0.0do60i=1,n60a(1)=a(1)+y(i)do80i=2,mma(i)=0.0do70j=1,n70a(I)=a(i)+x(i-1,j)*y(j)80continuecallCHOLESKY(b,mm,1,a,l)yy=0.0do90i=1,n90yy=yy+y(i)/nq=0.0dyy=0.0u=0.0ccccccccccccccccccccccccccccccccccdo110i=1,np=a(1)do100j=1,m100p=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)110continueccccccccccccccccccccccccccccccccccs=sqrt(q/n)r=sqrt(1.0-q/dyy)do150j=1,mp=0.0do140i=1,npp=a(1)do130k=1,mif(k.ne.j)pp=pp+a(k+1)*x(k,i)130continuep=p+(y(i)-pp)*(y(i)-pp)140continuev(j)=sqrt(1.0-q/p)150continuereturnendsubroutineCHOLESKY(a,n,m,d,l)!PerformtheCHOLESKYDecompositiondimensiona(n,n),d(n,m)l=1if(a(1,1)+1.0.eq.1.0)thenl=0write(*,30)returnendifa(1,1)=sqrt(a(1,1))do10j=2,n10a(1,j)=a(1,j)/a(1,1)do100i=2,ndo20j=2,i20a(i,i)=a(i,i)-a(j-1,i)*a(j-1,i)if(a(i,i)+1.0.eq.1.0)thenl=0write(*,30)returnendif30format(1x,'fail')a(i,i)=sqrt(a(i,i))if(i.ne.n)thendo50j=I+1,ndo40k=2,i40a(i,j)=a(i,j)-a(k-1,i)*a(k-1,j)50a(i,j)=a(i,j)/a(I,i)endif100continuedo130j=1,md(1,j)=d(1,j)/a(1,1)do120i=2,ndo110k=2,i110d(i,j)=d(i,j)-a(k-1,i)*d(k-1,j)d(i,j)=d(i,j)/a(i,i)120continue130continuedo160j=1,md(n,j)=d(n,j)/a(n,n)do150k=n,2,-1do140i=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)150continue160continuereturnend表一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年为极大值，但是估计值却为极小值，且各个年份的两个值的差别也很大，可以看出，通过回归方程预测的值不是很理想。通过表二可以看出，预报观测值与回归方程的预测值差别十分显著，进一步表明了多元线形回归方程的预测值并不理想。