2014年氣象統(tǒng)計(jì)實(shí)習(xí)報(bào)告

1、實(shí)習(xí)報(bào)告書課程名稱：氣象統(tǒng)計(jì)方法課程實(shí)踐姓名：班級：級氣科班*實(shí)習(xí)一求500hPa高度場氣候場、距平場和均方差場實(shí)習(xí)時(shí)間：第9周周三1、2節(jié)資料介紹有一500hPa高度場資料，文件名h500.dat,范圍：60150昌040N.時(shí)段：1982.11985.1狹48個(gè)月。水平分辨率：2.5*2.5，格點(diǎn)數(shù)：37*17。要求編fortran程序，求500hPa高度場的(1) 氣候場；(2) 距平場；(3) 均方差場。并能用Grads做出圖形，實(shí)習(xí)報(bào)告中氣候場、距平場、均方差場任意給出兩張圖，圖注要活楚，即要注明是哪個(gè)時(shí)間的圖形，并做簡單分析。注：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)!計(jì)算氣候場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運(yùn)行結(jié)果：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運(yùn)行結(jié)果：7jStt.grdQh500g|h500gjhSOOwritexd.gmt日xd上面圖中只展示了部分而未全部添加*實(shí)習(xí)二計(jì)算給定數(shù)據(jù)資料的簡單相關(guān)系數(shù)和自相關(guān)系數(shù)實(shí)習(xí)時(shí)間：第10周周三1、2節(jié)根據(jù)下表中年平均氣溫和冬季平均氣溫的等級數(shù)據(jù)進(jìn)行下列計(jì)算：1）計(jì)算兩個(gè)氣溫之間的簡單相關(guān)系數(shù)。2）分別找出兩個(gè)氣溫?cái)?shù)據(jù)自相關(guān)系數(shù)絕對值最大的滯后時(shí)間長度。（滯后長度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要求：實(shí)習(xí)報(bào)告中附出簡單相關(guān)系數(shù)或自相關(guān)系數(shù)程序答案：r=0.47年平均氣溫在滯后長度j=7、冬季序列在j=4最大。FORTRAN計(jì)算相關(guān)系數(shù)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)!求兩數(shù)組平均值DOi=1,Navr_y=avr_y+y(i)avr_w=avr_w+w(i)ENDDOavr_y=avr_y/Navr_w=avr_w/N!簡單相關(guān)系數(shù)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、均氣溫與冬季平均氣溫的簡單相關(guān)系數(shù)rxy=',f5.2)",rxyk=0!自相關(guān)系數(shù)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"侔年平均氣溫絕對值最大自相關(guān)系數(shù)rxy_ty=',f7.

16、4,/,'滯后時(shí)間長度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"佟季平均氣溫絕對值最大自相關(guān)系數(shù)rxy_tw=',f7.4,/,

17、'滯后時(shí)間長度k=',I2)",rxy_tw(k),kk=0END運(yùn)行成果：*實(shí)習(xí)四求給定數(shù)據(jù)的一元線性回歸方程實(shí)習(xí)時(shí)間：第12周周三1、2節(jié)利用下表數(shù)據(jù)，以環(huán)流指標(biāo)為預(yù)報(bào)因子，氣溫為預(yù)報(bào)量，計(jì)算氣溫和環(huán)流指標(biāo)之間的一元線性回歸方程，并對回歸方程進(jìn)行檢驗(yàn)。環(huán)流指標(biāo)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運(yùn)行成果：C:UsersadminDocumentsVisualStudio2010Ptb=-0.2342303b0=7.5055%y=7.5095H6-0.2

21、342993xF=20.4052>FQ=L41,回歸方程顯著請按任意鍵壅續(xù).*實(shí)習(xí)七計(jì)算給定數(shù)據(jù)的11年滑動(dòng)平均和累積距平實(shí)習(xí)時(shí)間：第15周周三1、2節(jié)利用數(shù)據(jù)ma.dat編寫11點(diǎn)滑動(dòng)平均的程序，ma.for給出了閱讀資料的fortran程序。數(shù)據(jù)在文件夾中單獨(dú)給出。要求：實(shí)習(xí)報(bào)告中附出程序，并給出原數(shù)據(jù)和滑動(dòng)后數(shù)據(jù)的圖形(1張圖)和累積距平數(shù)據(jù)圖形(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圖形顯示：原數(shù)據(jù)和滑動(dòng)后數(shù)據(jù)4JDGOOOQ站OODOO30000002.50000020000001.500000mooooo0500000O.OGOOOC原數(shù)據(jù)和累積距平數(shù)B.DOOOOO4.0000002.0000000,000000-2.000000-4,000000-6.000000-8.000000*實(shí)

24、習(xí)八對給定的海溫?cái)?shù)據(jù)進(jìn)行EO吩析實(shí)習(xí)時(shí)間：第16周周三1、2節(jié)給出海表溫度距平數(shù)據(jù)資料sstpx.grq以及相應(yīng)的數(shù)據(jù)描述文件sstpx.ct|對其進(jìn)行EO奇析，資料的時(shí)空范圍可以根據(jù)sstpx.ct獲知。數(shù)據(jù)在文件夾中單獨(dú)給出，距平或者標(biāo)準(zhǔn)化距平處理后再進(jìn)行EOFZhunsst.fo給出了如何讀取資料，Ssteof.fo為對距平或者標(biāo)準(zhǔn)化距平處理后的資料進(jìn)行EO吩析。要求：實(shí)習(xí)報(bào)告中給出第一特征向量及其時(shí)間系數(shù)，并分析其時(shí)空特征。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運(yùn)行程序后：得到zhunss

28、t.gr叔件，即圖中standardjdyform<hmeight名稱"Q)sstEof功SSTPXisstpx.grdJJJstandard,grd-BlZhunssti下面為老師給出程序，為運(yùn)行方便將部分簡略。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