2014年气象统计实习报告

1、实习报告书课程名称：气象统计方法课程实践姓名：班级：级气科班*实习一求500hPa高度场气候场、距平场和均方差场实习时间：第9周周三1、2节资料介绍有一500hPa高度场资料，文件名h500.dat,范围：60150昌040N.时段：1982.11985.1狭48个月。水平分辨率：2.5*2.5，格点数：37*17。要求编fortran程序，求500hPa高度场的(1) 气候场；(2) 距平场；(3) 均方差场。并能用Grads做出图形，实习报告中气候场、距平场、均方差场任意给出两张图，图注要活楚，即要注明是哪个时间的图形，并做简单分析。注：h500.Fori合出了如何用fortran读取AS

2、C码资料h500.dat.FORTRANprogramsx1implicitnoneintegernx,ny,mo,yrparameter(nx=37,ny=17,mo=12,yr=4)realvar(nx,ny,mo,yr)realat(nx,ny,mo),xd(nx,ny,mo,yr),sx(nx,ny,mo)integeri,j,m,t,it,iy,irecopen(5,file='d:studyformshixioneh500.dat')doiy=1,4dom=1,12read(5,1000)read(5,3000)(var(i,j,m,iy),i=1,nx),j=1,

3、ny)enddoenddoclose(5)!计算气候场atdot=1,12doj=1,nydoi=1,nxat(i,j,t)=0doit=1,4at(i,j,t)=at(i,j,t)+var(i,j,t,it)enddoat(i,j,t)=at(i,j,t)/4enddoenddoenddo!求距平场xddot=1,12doj=1,nydoi=1,nxxd(i,j,t,1)=0doit=1,4xd(i,j,t,it)=var(i,j,t,it)-at(i,j,t)enddoenddoenddoenddo!求均方差场sxdot=1,12doj=1,nydoi=1,nxsx(i,j,t)=0doi

4、t=1,4sx(i,j,t)=sx(i,j,t)+(var(i,j,t,it)-at(i,j,t)*2enddosx(i,j,t)=sqrt(sx(i,j,t)/4)enddoenddoenddo!写入气候场open(10,file='d:studyformshixioneat.grd',form='unformatted',access='direct',recl=nx*ny)irec=0dot=1,12irec=irec+1write(10,rec=irec)(at(i,j,t),i=1,nx),j=1,ny)enddoclose(10)!写

5、入距平场open(11,file='d:studyformshixionexd.grd',form='unformatted',access,'direct',recl=nx*ny)irec=0doit=1,4dot=1,12irec=irec+1write(11,rec=irec)(xd(i,j,t,it),i=1,nx),j=1,ny)enddoenddoclose(11)!写入均方差场open(12,file='d:studyformshixionesx.grd',form='unformatted',acc

6、ess,'direct',recl=nx*ny)irec=0dot=1,12irec=irec+1write(12,rec=irec)(sx(i,j,t),i=1,nx),j=1,ny)enddoclose(12)1000format(2i7)2000format(37F6.2)3000format(37f8.1)4000format(37f7.2)endprogramsx1运行结果：WtJat.9rd|h500gh500Plh500writesx.gmfJsxxd.gmfIxdat.gmf-1atH500.DATJh500di1wJsx.grd_IxdVxd.grdGrad或

7、;件气候场'reinit''enableprintd:studyformshixioneat.gmf''opend:studyformshixioneat.ctl''setgridoff'setgradsoff'setlat040''setlon60150''setlev500'mon=1while(mon<=12)'sett'mon'''dh''drawtitle1982year'mon'month

8、9;'print''c'mon=mon+1endwhile'disableprint'J距平场reinit'enableprintd:studyformshixionesx.gmf''opend:studyformshixionesx.ctl''setgridoff'setgradsoff'setlat040''setlon60150''setlev500'year=1982while(year<=1984)mon=1while(mon<=

9、12)'sett'mon'''dh''drawtitle500hPa'year'year'mon'monthanomaly'print''c'mon=mon+1endwhileyear=year+1endwhile'disableprint'SOOhFa1982year1monthonomaly均万差'reinit''enableprintd:studyformshixionesx.gmf''opend:studyfor

10、mshixionesx.ctl''setgridoff'setgradsoff'setlat040''setlon60150''setlev500'mon=1while(mon<=12)'sett'mon'''dh'drawtitle500hPa1982year'mon'monthMean-squareDeviation'print''c'mon=mon+1endwhile'disableprint'Vat

11、.gmfIatH500.DAT-JhSOOcIiHxsx.grdIxdxd.grd运行结果：7jStt.grdQh500g|h500gjhSOOwritexd.gmt日xd上面图中只展示了部分而未全部添加*实习二计算给定数据资料的简单相关系数和自相关系数实习时间：第10周周三1、2节根据下表中年平均气温和冬季平均气温的等级数据进行下列计算：1）计算两个气温之间的简单相关系数。2）分别找出两个气温数据自相关系数绝对值最大的滞后时间长度。（滞后长度r最大取10）中国1970T989年年平蝴冬季平均气温钢1970-1979V40侦120290280360却。2803.001980-19893JO3.

12、001902.703.50203.101802.902.90轴冬季我气温197M9793.243.143.262.383322.712.843.942.751.831980-19892.802.812.633.203.603.401071.872.632.47要求：实习报告中附出简单相关系数或自相关系数程序答案：r=0.47年平均气温在滞后长度j=7、冬季序列在j=4最大。FORTRAN计算相关系数rPROGRAMEXAMIMPLICITNONEINTEGER,PARAMETER:N=20INTEGERi,j,k,ty,tw,tywREAL:avr_y=0,avr_w=0,sy=0,sw=0,

13、rxy=0,max_y=0,max_w=0,max_yw=0REALy(N),w(N)DATAy/3.4,3.3,3.2,2.9,3.4,2.8,3.6,3.0,2.8,3.0,3.1,3.0,2.9,2.7,3.5,3.2,3.1,2.8,2.9,2.9/DATAw/3.24,3.14,3.26,2.38,3.32,2.71,2.84,3.94,2.75,1.83,2.80,2.81,2.63,3.20,3.60,3.40,3.07,1.87,2.63,2.47/REALsyy(N),sww(N),r(N),rty(N),rtw(N),rtyw(N),rxy_ty(N),rxy_tw(N),

14、rxy_tyw(N)!求两数组平均值DOi=1,Navr_y=avr_y+y(i)avr_w=avr_w+w(i)ENDDOavr_y=avr_y/Navr_w=avr_w/N!简单相关系数DOj=1,Nsyy(j)=(y(j)-avr_y)*2sy=sy+syy(j)sww(j)=(w(j)-avr_w)*2sw=sw+sww(j)ENDDOsy=sqrt(sy/N)sw=sqrt(sw/N)DOj=1,Nr(j)=(y(j)-avr_y)/sy)*(w(j)-avr_w)/sw)rxy=rxy+r(j)ENDDOrxy=rxy/NPRINT"(/'1970-198外全年平

15、均气温与冬季平均气温的简单相关系数rxy=',f5.2)",rxyk=0!自相关系数DOty=1,N/2DOi=1,N-tyrty(i)=(y(i)-avr_y)/sy)*(y(i+ty)-avr_y)/sy)rxy_ty(ty)=rxy_ty(ty)+rty(i)ENDDOrxy_ty(ty)=rxy_ty(ty)/(N-ty)rxy_ty(ty)=ABS(rxy_ty(ty)IF(rxy_ty(ty)>max_y)THENmax_y=rxy_ty(ty)k=tyENDIFENDDOPRINT"侔年平均气温绝对值最大自相关系数rxy_ty=',f7.

16、4,/,'滞后时间长度k=',I2)",rxy_ty(k),kk=0DOtw=1,N/2DOi=1,N-twrtw(i)=(w(i)-avr_w)/sw)*(w(i+tw)-avr_w)/sw)rxy_tw(tw)=rxy_tw(tw)+rtw(i)ENDDOrxy_tw(tw)=rxy_tw(tw)/(N-tw)rxy_tw(tw)=ABS(rxy_tw(tw)IF(rxy_tw(tw)>max_w)THENmax_w=rxy_tw(tw)k=twENDIFENDDOPRINT"佟季平均气温绝对值最大自相关系数rxy_tw=',f7.4,/,

17、'滞后时间长度k=',I2)",rxy_tw(k),kk=0END运行成果：*实习四求给定数据的一元线性回归方程实习时间：第12周周三1、2节利用下表数据，以环流指标为预报因子，气温为预报量，计算气温和环流指标之间的一元线性回归方程，并对回归方程进行检验。环流指标19510.93219521.22519532.22019542.4261955-0.52719562.5241957-1.128195802419596.21519602.71619613.2241962-1.130答案：?7.5-0.23xF=20.18>Fa=4.41,回归方程显著programs

18、x4implicitnone19632.52219641.23019651.82419660.63319672.42619682.52019691.2321970-0.835integerireal:x(20),y(20),f1=0,f2=0,b,avx,avy,b0,sx=0,sy=0,sxy=0,f,rxydata(x(i),i=1,20)/32,25,20,26,27,24,28,24,15,16,24,30,22,30,24,33,26,20,32,35/data(y(i),i=1,20)/0.9,1.2,2.2,2.4,-0.5,2.5,-1.1,0,6.2,2.7,3.2,-1.1

19、,2.5,1.2,1.8,0.6,2.4,2.5,1.2,-0.8/doi=1,20f1=f1+x(i)*y(i)enddoavx=sum(x)/20!求均值|avy=sum(y)/20doi=1,20f2=f2+x(i)*2enddob=(20*f1-sum(x)*sum(y)/(20*f2-sum(x)*2)求？bb0=avy-b*avxprint*,'b=',bprint*,'b0=',b0print*,'y=',b0,b,'x'doi=1,20sx=sx+(x(i)-avx)*2sy=sy+(y(i)-avy)*2sxy=

20、sxy+(x(i)-avx)*(y(i)-avy)enddosx=sqrt(sx/20)sy=sqrt(sy/20)sxy=sxy/20rxy=sx*b/syf=rxy*2*18/(1-rxy*2)if(f>4.41)print*,'F=',f,'>F也=4.41,回归方程显著if(f<=4.41)print*,'F=',f,'<F也=4.41,0归方程不显著End运行成果：C:UsersadminDocumentsVisualStudio2010Ptb=-0.2342303b0=7.5055%y=7.5095H6-0.2

21、342993xF=20.4052>FQ=L41,回归方程显著请按任意键壅续.*实习七计算给定数据的11年滑动平均和累积距平实习时间：第15周周三1、2节利用数据ma.dat编写11点滑动平均的程序，ma.for给出了阅读资料的fortran程序。数据在文件夹中单独给出。要求：实习报告中附出程序，并给出原数据和滑动后数据的图形(1张图)和累积距平数据图形(1张图)programsx7implicitnoneintegern,ih,nyearparameter(n=85,ih=11,nyear=1922)real:x(n),x1(n-ih+1)=0,x2(n)=0,averx=0intege

22、ri,jopen(2,file='d:studyformshiyansevenma.dat')open(3,file='d:studyformshiyansevenmaresult.dat')open(4,file='d:studyformshiyansevenaccresult.dat')read(2,*)(x(i),i=1,n)doj=1,n-ih+1!movingaveragedoi=1,ihx1(j)=x1(j)+x(i+j-1)enddox1(j)=x1(j)/ihenddodoi=1,n-ih+1write(3,*)x1(i)endd

23、odoi=1,n!averageofxaverx=averx+x(i)enddoaverx=averx/ndoi=1,n!accelaratedoj=1,ix2(i)=x2(i)+(x(j)-averx)enddoenddodoi=1,nwrite(4,*)x2(i)enddoEnd图形显示：原数据和滑动后数据4JDGOOOQ站OODOO30000002.50000020000001.500000mooooo0500000O.OGOOOC原数据和累积距平数B.DOOOOO4.0000002.0000000,000000-2.000000-4,000000-6.000000-8.000000*实

24、习八对给定的海温数据进行EO吩析实习时间：第16周周三1、2节给出海表温度距平数据资料sstpx.grq以及相应的数据描述文件sstpx.ct|对其进行EO奇析，资料的时空范围可以根据sstpx.ct获知。数据在文件夹中单独给出，距平或者标准化距平处理后再进行EOFZhunsst.fo给出了如何读取资料，Ssteof.fo为对距平或者标准化距平处理后的资料进行EO吩析。要求：实习报告中给出第一特征向量及其时间系数，并分析其时空特征。programmainparameter(mo=12,my=43,nx=18,ny=12,mt=516)dimensionavf(mo,nx,ny),df(mo,n

25、x,ny)dimensionsst(nx,ny,mt),sst3(nx,ny,mt)open(1,file='d:studyformshixieightsstpx.grd',form='unformatted',access='direct',recl=nx*ny)irec=1doit=1,mtread(1,rec=irec)(sst(i,j,it),i=1,nx),j=1,ny)irec=irec+1enddo!averagedoj=1,nydoi=1,nxdok=1,modoit=k,mt,12avf(k,i,j)=avf(k,i,j)+ss

26、t(i,j,it)enddoavf(k,i,j)=avf(k,i,j)/myenddoenddoenddo!variancedoj=1,nydoi=1,nxdok=1,modoit=k,mt,12df(k,i,j)=df(k,i,j)+(sst(i,j,it)-avf(k,i,j)*2enddodf(k,i,j)=sqrt(df(k,i,j)/my)enddoenddoenddo!standardizingdoj=1,nydoi=1,nxdok=1,modoit=k,mt,12if(sst(i,j,it)=-999.0)thensst3(i,j,it)=-999.0elsesst3(i,j,i

27、t)=(sst(i,j,it)-avf(k,i,j)/df(k,i,j)endifenddoenddoenddoenddo!outputfileopen(2,file='d:studyformshixieightstandard.grd',form='unformatted',access='direct',recl=nx*ny)irec=1doit=1,mtwrite(2,rec=irec)(sst3(i,j,it),i=1,nx),j=1,ny)irec=irec+1enddoclose(2)close(1)end运行程序后：得到zhunss

28、t.gr叔件，即图中standardjdyform<hmeight名称"Q)sstEof功SSTPXisstpx.grdJJJstandard,grd-BlZhunssti下面为老师给出程序，为运行方便将部分简略。PROGRAMOFPWPARAMETER=516,N=216,MNH=216,KS=-1,KV=10,KVT=2)PARAMET昨-999.0,nx=18,ny=12)DIMENSIONF(N,M),A(MNH,MNH),S(MNH,MNH),ER(mnh,4),DF(N),V(MNH),AVF(N),evf(N,KVTINXfn(M,KVT),open(11,fil

29、e='d:studyformshixieightzhunsst.grd',form='unformatted',access='direct',recl=nx*ny)do132it=1,mread(11,rec=it)(data(i,j),i=1,nx),j=1,ny)do132jj=1,nydo132ii=1,nxkkkk=nx*(jj-1)+iif(kkkk,it)=data(ii,jj)132continueclose(11)CALTest1(n,m,ff,f,nf)write(*,*)'ok2'CALTRANSF(N,M,

30、F,nf,AVF,DF,KS)write(*,*)'ok3'CALFORMA(N,M,MNH,F,A)write(*,*)'ok4'CALJCB(MNH,A,S,0.00001)write(*,*)'ok5'CALIARRANG(KV,MNH,A,ER,S)write(*,*)'ok6'CALTCOEFF(KVT,KV,N,M,MNH,S,F,V,evf,tcf,ER)write(*,*)'ok7'calltest3(N,ff,nf,evf,kvt)write(*,*)'ok8'open(21,fi

31、le='d:studyformshixieightevf.grd',form='unformatted',access='direct',recl=nx*ny)irec=0do668kk=1,kvtirec=irec+1668write(21,rec=irec)(evf(nx*(j-1)+i,kk),i=1,nx),j=1,ny)close(21)open(21,file='d:studyformshixieighttcf.grd',form='unformatted',access='direct'

32、;,recl=kvt)irec=0do345it=1,mirec=irec+1345write(21,rec=irec)(tcf(it,iik),iik=1,kvt)close(21)106format(10f8.4)open(21,file='d:studyformshixieightdats.dat')write(21,106)(er(iiii,3),iiii=1,kv)close(21)STOPENDSUBROUTINEst1(n1,m,ff,f,nf)DIMENSIONN1,M)DIMENSIOF(NI)doi=1,n1nf(i)=0.0enddodoi=1,n1dok

33、=1,mif(f(i,k).eq.ff)thenf(i,k)=0.0nf(i)=nf(i)+1endifenddoenddoreturnendSUBROUTINRANSF(n1,m,f,nf,avf,df,ks)DIMENSIONN1,M),AVF(N1)DIMENSIONF(N1)DIMENSIOF(N1)if(ks.eq.-1thengoto30endifdoi=1,n1avf(i)=0.0enddoif(ks.eq.0)hengoto5endifdoi=1,n1df(i)=0.0enddo5continueDO141I=1,N1if(nf(i).ne.0)goto141do12j=1,m

34、12AVF(I)=AVF(I)+F(I,J)AVF(I)=AVF(I)/MDO14J=1,M14F(I,J)=F(I,J)-AVF(I)141CONTINUEIF(KS.EQ.0l)HENRETURNELSEDO241I=1,N1if(nf(i).ne.0)goto241DO22J=1,M22DF(I)=DF(I)+F(I,J)*F(I,J)DF(I)=SQRT(DF(I)/M)DO24J=1,M24F(I,J)=F(I,J)/DF(I)241CONTINUEENDIF30CONTINUERETURNENDSUBROUTINORMA(N,M,MNH,F,A)DIMENSIONN,M),A(MN

35、H,MNH)IF(M-N)40,50,5040DO44I=1,MNHDO44J=I,MNHA(I,J)=0.0DO42IS=1,N42A(I,J)=A(I,J)+F(IS,I)*F(IS,J)A(J,I)=A(I,J)44CONTINUERETURN50DO54I=1,MNHDO54J=I,MNHA(I,J)=0.0DO52JS=1,M52A(I,J)=A(I,J)+F(I,JS)*F(J,JS)A(J,I)=A(I,J)54CONTINUERETURNENDSUBROUTINEB(N,A,S,EPS)DIMENSIONN,N),S(N,N)DO30I=1,NDO30J=1,I10203040

36、5060IF(I-J)20,10,20S(I,J)=1.GOTCB0S(I,J)=0.S(J,I)=0.CONTINUEG=0.DO40I=2,NI1=I-1DO40J=1,I1G=G+2.*A(I,J)*A(I,J)S1=SQRT(G)S2=EPS/FLOAT(N)*S1S3=S1L=0S3=S3/FLOAT(N)DO130IQ=2,NIQ1=IQ-1DO130IP=1,IQ1IF(ABS(A(IP,IQ).LT.GOTQ30L=1V1=A(IP,IP)V2=A(IP,IQ)V3=A(IQ,IQ)U=0.5*(V1-V3)IF(U.EQ.0.0)G=1.IF(ABS(U).GE.1E-10)

37、G=-SIGN(1.,U)*V2/SQRT(V2*V2+U*U)ST=G/SQRT(2.*(1.+SQRT(1.-G*G)CT=SQRT(1.-ST*ST)DO110I=1,NG=A(I,IP)*CT-A(I,IQ)*STA(I,IQ)=A(I,IP)*ST+A(I,IQ)*CTA(I,IP)=GG=S(I,IP)*CT-S(I,IQ)*STS(I,IQ)=S(I,IP)*ST+S(I,IQ)*CT110S(I,IP)=GDO120I=1,NA(IP,I)=A(I,IP)120A(IQ,I)=A(I,IQ)G=2.*V2*ST*CTA(IP,IP)=V1*CT*CT+V3*ST*ST-GA(I

38、Q,IQ)=V1*ST*ST+V3*CT*CT+GA(IP,IQ)=(V1-V3)*ST*CT+V2*(CT*CT-ST*ST)A(IQ,IP)=A(IP,IQ)130CONTINUEIF(L-1)150,140,150140L=0GOT(60150IF(S3.GT.S2)OTO50RETURNENDSUBROUTIARRANG(KV,MNH,A,ER,S)DIMENSIONMNH,MNH),ER(mnh,4),S(MNH,MNH)TR=0.0DO200I=1,MNHTR=TR+A(I,I)200ER(I,1)=A(I,I)MNH仁MNH-1DO210K1=MNH1,1,-1DO210K2=K

39、1,MNH1IF(ER(K2,1).LT.ER(K2+1THENC=ER(K2+1,1)ER(K2+1,1)=ER(K2,1)ER(K2,1)=CDO205I=1,MNHC=S(I,K2+1)S(I,K2+1)=S(I,K2)S(I,K2)=C205CONTINUEENDIF210CONTINUEER(1,2)=ER(1,1)DO220I=2,KVER(I,2)=ER(I-1,2)+ER(I,1)220CONTINUEDO230I=1,KVER(I,3)=ER(I,1)/TRER(I,4)=ER(I,2)/TR230CONTINUEWRIT(6,250)TR250FORMAT5X,'TOTALSQUAREERROR=',F20.5)RETURNENDSUBROUTINEOEFF(KVT,KV,N,M,MNH,S,F,V,evf,tcf,ER)DI