2023年气象统计实习报告

实习报告书课程名称：气象统计方法课程实践姓名：学号：班级：级气科班\\*实习一求500hPa高度场气候场、距平场和均方差场实习时间：第9周周三1、2节资料介绍有一500hPa高度场资料，文件名h500.dat，范围：60~150E，0~40N.时段：1982.1~1985.12共48个月。水平分辨率：2.5*2.5，格点数：37*17。2．要求编fortran程序，求500hPa高度场的（1）气候场；（2）距平场；（3）均方差场。并能用Grads做出图形,实习报告中气候场、距平场、均方差场任意给出两张图，图注要清楚，即要注明是哪个时间的图形，并做简单分析。注：h500.For给出了如何用fortran读取ASCII码资料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:\study\form\shixione\h500.dat')doiy=1,4dom=1,12read(5,1000)read(5,3000)((var(i,j,m,iy),i=1,nx),j=1,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)=0doit=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:\study\form\shixione\at.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)!写入距平场open(11,file='d:\study\form\shixione\xd.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:\study\form\shixione\sx.grd',form='unformatted',access='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运行结果：Grads文件气候场'reinit''enableprintd:\study\form\shixione\at.gmf''opend:\study\form\shixione\at.ctl''setgridoff''setgradsoff''setlat040''setlon60150''setlev500'mon=1while(mon<=12)'sett'mon'''dh''drawtitle1982year'mon'month''print''c'mon=mon+1endwhile'disableprint';距平场'reinit''enableprintd:\study\form\shixione\sx.gmf''opend:\study\form\shixione\sx.ctl''setgridoff''setgradsoff''setlat040''setlon60150''setlev500'year=1982while(year<=1984)mon=1while(mon<=12)'sett'mon'''dh''drawtitle500hPa'year'year'mon'monthanomaly''print''c'答案：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,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),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-1989年全年平均气温与冬季平均气温的简单相关系数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.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,/,'滞后时间长度k=',I2)",rxy_tw(k),kk=0END运行成果：*实习四求给定数据的一元线性回归方程实习时间：第12周周三1、2节利用下表数据，以环流指标为预报因子，气温为预报量，计算气温和环流指标之间的一元线性回归方程，并对回归方程进行检验。年份气温T环流指标19510.93219521.22519532.22019542.4261955-0.52719562.5241957-1.128195802419596.21519602.71619613.2241962-1.13019632.52219641.23019651.82419660.63319672.42619682.52019691.2321970-0.835答案：F=20.18>Fα=4.41，回归方程显著programsx4implicitnoneintegerireal::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,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=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,回归方程不显著End运行成果：*实习七计算给定数据的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=0integeri,jopen(2,file='d:\study\form\shiyanseven\ma.dat')open(3,file='d:\study\form\shiyanseven\maresult.dat')open(4,file='d:\study\form\shiyanseven\accresult.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)enddodoi=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图形显示：*实习八对给定的海温数据进行EOF分析实习时间：第16周周三1、2节给出海表温度距平数据资料sstpx.grd，以及相应的数据描述文件sstpx.ctl，对其进行EOF分析，资料的时空范围可以根据sstpx.ctl获知。数据在文件夹中单独给出，距平或者标准化距平处理后再进行EOF。Zhunsst.for给出了如何读取资料，Ssteof.for为对距平或者标准化距平处理后的资料进行EOF分析。要求：实习报告中给出第一特征向量及其时间系数，并分析其时空特征。programmainparameter(mo=12,my=43,nx=18,ny=12,mt=516)dimensionavf(mo,nx,ny),df(mo,nx,ny)dimensionsst(nx,ny,mt),sst3(nx,ny,mt)open(1,file='d:\study\form\shixieight\sstpx.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)+sst(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,it)=(sst(i,j,it)-avf(k,i,j))/df(k,i,j)endifenddoenddoenddoenddo!outputfileopen(2,file='d:\study\form\shixieight\standard.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运行程序后：得到zhunsst.grd文件，即图中standard下面为老师给出程序，为运行方便将部分简略。PROGRAMEOFPWPARAMETER(M=516,N=216,MNH=216,KS=-1,KV=10,KVT=2)PARAMETER(ff=-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,KVT),tCF(M,KVT),data(Nx,ny),nf(N)open(11,file='d:\study\form\shixieight\zhunsst.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)CALLTest1(n,m,ff,f,nf)write(*,*)'ok2'CALLTRANSF(N,M,F,nf,AVF,DF,KS)write(*,*)'ok3'CALLFORMA(N,M,MNH,F,A)write(*,*)'ok4'CALLJCB(MNH,A,S,0.00001)write(*,*)'ok5'CALLARRANG(KV,MNH,A,ER,S)write(*,*)'ok6'CALLTCOEFF(KVT,KV,N,M,MNH,S,F,V,evf,tcf,ER)write(*,*)'ok7'calltest3(N,ff,nf,evf,kvt)write(*,*)'ok8'open(21,file='d:\study\form\shixieight\evf.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:\study\form\shixieight\tcf.grd',form='unformatted',access='direct',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:\study\form\shixieight\dats.dat')write(21,106)(er(iiii,3),iiii=1,kv)close(21)STOPENDSUBROUTINETest1(n1,m,ff,f,nf)DIMENSIONF(N1,M)DIMENSIONnF(N1)doi=1,n1nf(i)=0.0enddodoi=1,n1dok=1,mif(f(i,k).eq.ff)thenf(i,k)=0.0nf(i)=nf(i)+1endifenddoenddoreturnendSUBROUTINETRANSF(n1,m,f,nf,avf,df,ks)DIMENSIONF(N1,M),AVF(N1)DIMENSIONDF(N1)DIMENSIONnF(N1)if(ks.eq.-1)thengoto30endifdoi=1,n1avf(i)=0.0enddoif(ks.eq.0)thengoto5endifdoi=1,n1df(i)=0.0enddo5continueDO141I=1,N1if(nf(i).ne.0)goto141do12j=1,m12AVF(I)=AVF(I)+F(I,J)AVF(I)=AVF(I)/MDO14J=1,M14F(I,J)=F(I,J)-AVF(I)141CONTINUEIF(KS.EQ.0)THENRETURNELSEDO241I=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)241CONTINUEENDIF30CONTINUERETURNENDSUBROUTINEFORMA(N,M,MNH,F,A)DIMENSIONF(N,M),A(MNH,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)54CONTINUERETURNENDSUBROUTINEJCB(N,A,S,EPS)DIMENSIONA(N,N),S(N,N)DO30I=1,NDO30J=1,IIF(I-J)20,10,2010S(I,J)=1.GOTO3020S(I,J)=0.S(J,I)=0.30CONTINUEG=0.DO40I=2,NI1=I-1DO40J=1,I140G=G+2.*A(I,J)*A(I,J)S1=SQRT(G)S2=EPS/FLOAT(N)*S1S3=S1L=050S3=S3/FLOAT(N)60DO130IQ=2,NIQ1=IQ-1DO130IP=1,IQ1IF(ABS(A(IP,IQ)).LT.S3)GOTO130L=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)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(IQ,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=0GOTO60150IF(S3.GT.S2)GOTO50RETURNENDSUBROUTINEARRANG(KV,MNH,A,ER,S)DIMENSIONA(MNH,MNH),ER(mnh,4),S(MNH,MNH)TR=0.0DO200I=1,MNHTR=TR+A(I,I)200ER(I,1)=A(I,I)MNH1=MNH-1DO210K1=MNH1,1,-1DO210K2=K1,MNH1IF(ER(K2,1).LT.ER(K2+1,1))THENC=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)/TR230CONTINUEWRITE(6,250)TR250FORMAT(/5X,'TOTALSQUAREERROR=',F20.5)RETURNENDSUBROUTINETCOEFF(KVT,KV,N,M,MNH,S,F,V,evf,tcf,ER)DIMENSIONS(MNH,MNH),F(N,M),V(MNH),ER(mnh,4),evf(n,kvt),tcf(m,kvt)DO360J=1,KVTC=0.DO350I=1,MNH350C=C+S(I,J)*S(I,J)C=SQRT(C)DO160I=1,MNHs(I,J)=S(I,J)/C160evf(I,J)=S(I,J)/C360CONTINUEIF(N.LE.M)THENDO390J=1,MDO370I=1,NV(I)=F(I,J)F(I,J)=0.370CONTINUEdo371is=1,kvttcf(j,is)=0.371continueDO380IS=1,KVTDO380I=1,Nf(is,j)=F(IS,J)+V(I)*S(I,IS)380tcf(j,is)=tcf(J,is)+V(I)*S(I,IS)390CONTINUEELSEDO410I=1,NDO400J=1,MV(J)=F(I,J)F(I,J)=0.400CONTINUEDO410JS=1,KVTDO410J=1,Mf(I,JS)=F(I,JS)+V(J)*S(J,JS)410CONTINUEDO430JS=1,KVTDO420J=1,Mtcf(J,JS)=S(J,JS)*SQRT(ER(JS,1))420CONTINUEDO430I=1,Nevf(I,JS)=F(I,JS)/SQRT(ER(JS,1))430CONTINUEt=0.0do3650i=1,m3650t=t+tcf(i,1)*tcf(i,2)write(*,*)tt=0.0do3651i=1,n3651t=t+evf(i,1)*evf(i,2)write(*,*)tENDIFRETURNENDSUBROUTINEtest3(N1,ff,nf,evf,kvt)dimensionnf(n1),evf(n1,kvt)doi=1,n1if(nf(i).ne.0)thendok=1,kvtevf(i,k)=ffenddoendifenddoreturnendSUBROUTINEOUTER(KV,ER,mnh)DIMENSIONER(mnh,4)WRITE(16,510)510FORMAT(/30X,'EIGENVALUEANDANALYSISERROR',/)WRITE(16,520)520FORMAT(10X,1HH,8X,5HLAMDA,10X,6HSLAMDA,11X,2HPH,12X,3HSPH)WRITE(16,530)(IS,(ER(IS,J),J=1,4),IS=1,KV)530FORMAT(1X,I10,4F15.5)WRITE(16,540)540FORMAT(//)RETURNENDSUBROUTINEOUTVT(KVT,N,M,MNH,S,F,EGVT,ECOF)DIMENSIONF(N,M),S(MNH,MNH),EGVT(N,KVT),ECOF(M,KVT)WRITE(16,560)560FORMAT(30X,'STANDARDEIGENVECTORS',/)WRITE(16,570)(IS,IS=1,KVT)570FORMAT(3X,80I7)DO550I=1,NIF(M.GE.N)THENWRITE(16,580)I,(S(I,JS),JS=1,KVT)580FORMAT(1X,I3,80F7.3,/)DO11JS=1,KVTEGVT(I,JS)=S(I,JS)11CONTINUEELSEWRITE(16,590)I,(F(I,JS),JS=1,KVT)590FORMAT(1X,I3,80F7.3)DO12JS=1,KVTEGVT(I,JS)=F(I,JS)12CONTINUEENDIF550CONTINUEWRITE(16,720)720FORMAT(//)WRITE(16,610)610FORMAT(30X,'TIME-COEFFICENTSERIESOFS.E.'/)WRITE(16,620)(IS,IS=1,KVT)620FORMAT(3X,80I7)DO600J=1,MIF(M.GE.N)THENWRITE(16,630)J,(F(IS,J),IS=1,KVT)630FORMAT(1X,I3,80F7.1)DO13IS=1,KVTECOF(J,IS)=F(IS,J)13CONTINUEELSEWRITE(16,640)J,(S(J,IS),IS=1,KVT)640FORMAT(1X,I3,80F7.1)DO14IS=1,KVTECOF(J,IS)=S(J,IS)14CO