2014年气象统计实习报告

实习报告书课程名称：气象统计方法课程实践姓名：学号： 班级： 级气科 班*实习一 求500hPa高度场气候场、距平场和均方差场 实习时间：第9周周三1、2节1. 资料介绍有一500hPa高度场资料，文件名h500.dat，范围：60150E，040N.时段：1982.11985.12共48个月。水平分辨率：2.5*2.5，格点数：37*17。2要求编fortran程序，求500hPa高度场的（1）气候场；（2）距平场；（3）均方差场。并能用Grads做出图形,实习报告中气候场、距平场、均方差场任意给出两张图，图注要清楚，即要注明是哪个时间的图形，并做简单分析。注：h500.For给出了如何用fortran读取ASCII码资料h500.dat.n FORTRANprogram sx1implicit noneinteger nx,ny,mo,yrparameter(nx=37,ny=17,mo=12,yr=4)real var(nx,ny,mo,yr)real at(nx,ny,mo),xd(nx,ny,mo,yr),sx(nx,ny,mo)integer i,j,m,t,it,iy,irecopen(5,file=d:studyformshixioneh500.dat) do iy=1,4 do m=1,12 read(5,1000) read(5,3000)(var(i,j,m,iy),i=1,nx),j=1,ny)enddoenddo close(5)!计算气候场atdo t=1,12 do j=1,ny do i=1,nx at(i,j,t)=0 do it=1,4 at(i,j,t)=at(i,j,t)+var(i,j,t,it) enddo at(i,j,t)=at(i,j,t)/4 enddo enddoenddo!求距平场xddo t=1,12 do j=1,ny do i=1,nx xd(i,j,t,1)=0 do it=1,4 xd(i,j,t,it)=var(i,j,t,it)-at(i,j,t) enddo enddo enddoenddo!求均方差场sxdo t=1,12 do j=1,ny do i=1,nx sx(i,j,t)=0 do it=1,4 sx(i,j,t)=sx(i,j,t)+(var(i,j,t,it)-at(i,j,t)*2 enddo sx(i,j,t)=sqrt(sx(i,j,t)/4) enddo enddoenddo!写入气候场open(10,file=d:studyformshixioneat.grd,form=unformatted,access=direct,recl=nx*ny)irec=0do t=1,12irec=irec+1write(10,rec=irec)(at(i,j,t),i=1,nx),j=1,ny)enddoclose(10)!写入距平场open(11,file=d:studyformshixionexd.grd,form=unformatted,access=direct,recl=nx*ny)irec=0do it=1,4 do t=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,access=direct,recl=nx*ny)irec=0do t=1,12irec=irec+1write(12,rec=irec)(sx(i,j,t),i=1,nx),j=1,ny)enddoclose(12)1000 format(2i7)2000 format(37F6.2)3000 format(37f8.1)4000 format(37f7.2)end program sx1n 运行结果：n Grads文件气候场reinitenable print d:studyformshixioneat.gmfopen d:studyformshixioneat.ctlset grid offset grads offset lat 0 40set lon 60 150set lev 500mon=1while(mon=12)set t mond hdraw title 1982yearmonmonthprintcmon=mon+1endwhiledisable print;距平场reinitenable print d:studyformshixionesx.gmfopen d:studyformshixionesx.ctlset grid offset grads offset lat 0 40set lon 60 150set lev 500year=1982while(year=1984)mon=1while(mon=12)set t mond hdraw title 500hPa yearyearmonmonth anomalyprintcmon=mon+1endwhileyear=year+1endwhiledisable print;均方差reinitenable print d:studyformshixionesx.gmfopen d:studyformshixionesx.ctlset grid offset grads offset lat 0 40set lon 60 150set lev 500mon=1while(monmax_y) THEN max_y=rxy_ty(ty) k=ty END IFEND DOPRINT (全年平均气温绝对值最大自相关系数rxy_ty=,f7.4,/,滞后时间长度 k=,I2),rxy_ty(k),kk=0DO tw=1,N/2 DO i=1,N-tw rtw(i)=(w(i)-avr_w)/sw)*(w(i+tw)-avr_w)/sw) rxy_tw(tw)=rxy_tw(tw)+rtw(i) END DO rxy_tw(tw)=rxy_tw(tw)/(N-tw) rxy_tw(tw)=ABS(rxy_tw(tw) IF(rxy_tw(tw)max_w) THEN max_w=rxy_tw(tw) k=tw END IFEND DOPRINT (冬季平均气温绝对值最大自相关系数rxy_tw=,f7.4,/,滞后时间长度 k=,I2),rxy_tw(k),kk=0ENDn 运行成果：*实习四 求给定数据的一元线性回归方程实习时间：第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.18F=4.41，回归方程显著program sx4implicit noneinteger ireal: 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/do i=1,20f1=f1+x(i)*y(i)end doavx=sum(x)/20 !求均值avy=sum(y)/20do i=1,20f2=f2+x(i)*2end dob=(20*f1-sum(x)*sum(y)/(20*f2-sum(x)*2) !求bb0=avy-b*avxprint *,b=,bprint *,b0=,b0print*,y=,b0,b,xdo i=1,20 sx=sx+(x(i)-avx)*2sy=sy+(y(i)-avy)*2sxy=sxy+(x(i)-avx)*(y(i)-avy)end do sx=sqrt(sx/20) sy=sqrt(sy/20)sxy=sxy/20rxy=sx*b/syf=rxy*2*18/(1-rxy*2)if(f4.41) print *,F=,f,F=4.41，回归方程显著if(f=4.41) print *,F=,f,F=4.41,回归方程不显著Endn 运行成果：*实习七 计算给定数据的11年滑动平均和累积距平 实习时间：第15周周三1、2节利用数据ma.dat，编写11点滑动平均的程序，ma.for给出了阅读资料的fortran程序。数据在文件夹中单独给出。要求：实习报告中附出程序，并给出原数据和滑动后数据的图形（1张图）和累积距平数据图形（1张图）program sx7 implicit noneinteger n,ih,nyearparameter(n=85,ih=11,nyear=1922)real:x(n),x1(n-ih+1)=0,x2(n)=0,averx=0integer i,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)do j=1,n-ih+1 !moving average do i=1,ih x1(j)=x1(j)+x(i+j-1) end do x1(j)=x1(j)/ihend dodo i=1,n-ih+1write(3,*) x1(i)end dodo i=1,n !average of x averx=averx+x(i)end doaverx=averx/ndo i=1,n !accelarate do j=1,i x2(i)=x2(i)+(x(j)-averx) end doend do do i=1,n write(4,*) x2(i)end do Endn 图形显示：*实习八 对给定的海温数据进行EOF分析 实习时间：第16周周三1、2节给出海表温度距平数据资料sstpx.grd，以及相应的数据描述文件sstpx.ctl，对其进行EOF分析，资料的时空范围可以根据sstpx.ctl获知。 数据在文件夹中单独给出，距平或者标准化距平处理后再进行EOF。Zhunsst.for给出了如何读取资料，Ssteof.for为对距平或者标准化距平处理后的资料进行EOF分析。要求：实习报告中给出第一特征向量及其时间系数，并分析其时空特征。 program main parameter(mo=12,my=43,nx=18,ny=12,mt=516) dimension avf(mo,nx,ny),df(mo,nx,ny) dimension sst(nx,ny,mt),sst3(nx,ny,mt) open(1,file=d:studyformshixieightsstpx.grd,form=unformatted,access=direct,recl=nx*ny) irec=1 do it=1,mt read(1,rec=irec)(sst(i,j,it),i=1,nx),j=1,ny) irec=irec+1 end do! average do j=1,ny do i=1,nx do k=1,mo do it=k,mt,12 avf(k,i,j)=avf(k,i,j)+sst(i,j,it) end do avf(k,i,j)=avf(k,i,j)/my end do end do end do! variance do j=1,ny do i=1,nx do k=1,mo do it=k,mt,12 df(k,i,j)=df(k,i,j)+(sst(i,j,it)-avf(k,i,j)*2 end do df(k,i,j)=sqrt(df(k,i,j)/my) end do end do end do ! standardizing do j=1,ny do i=1,nx do k=1,mo do it=k,mt,12 if(sst(i,j,it)=-999.0)then sst3(i,j,it)=-999.0 else sst3(i,j,it)=(sst(i,j,it)-avf(k,i,j)/df(k,i,j) end if end do end do end doend do! output file open(2,file=d:studyformshixieightstandard.grd,form=unformatted,access=direct,recl=nx*ny) irec=1 do it=1,mt write(2,rec=irec)(sst3(i,j,it),i=1,nx),j=1,ny) irec=irec+1 end do close(2) close(1)endn 运行程序后：得到zhunsst.grd文件，即图中standard下面为老师给出程序，为运行方便将部分简略。 PROGRAM EOFPW PARAMETER(M=516,N=216,MNH=216,KS=-1,KV=10,KVT=2) PARAMETER(ff=-999.0,nx=18,ny=12) DIMENSION F(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:studyformshixieightzhunsst.grd,form=unformatted,access=direct,recl=nx*ny) do 132 it=1,m read(11,rec=it)(data(i,j),i=1,nx),j=1,ny) do 132 jj=1,ny do 132 ii=1,nx kkkk=nx*(jj-1)+ii f(kkkk,it)=data(ii,jj)132 continue close(11) CALL Test1(n,m,ff,f,nf)write(*,*)ok2 CALL TRANSF(N,M,F,nf,AVF,DF,KS) write(*,*)ok3 CALL FORMA(N,M,MNH,F,A)write(*,*)ok4 CALL JCB(MNH,A,S,0.00001)write(*,*)ok5 CALL ARRANG(KV,MNH,A,ER,S)write(*,*)ok6 CALL TCOEFF(KVT,KV,N,M,MNH,S,F,V,evf,tcf,ER)write(*,*)ok7 call test3(N,ff,nf,evf,kvt)write(*,*)ok8 open(21,file=d:studyformshixieightevf.grd,form=unformatted,access=direct,recl=nx*ny) irec=0 do 668 kk=1,kvt irec=irec+1668 write(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,recl=kvt) irec=0 do 345 it=1,m irec=irec+1345 write(21,rec=irec) (tcf(it,iik),iik=1,kvt) close(21)106 format(10f8.4) open(21,file=d:studyformshixieightdats.dat) write(21,106)(er(iiii,3),iiii=1,kv) close(21) STOP END SUBROUTINE Test1(n1,m,ff,f,nf) DIMENSION F(N1,M) DIMENSION nF(N1) do i=1,n1 nf(i)=0.0 enddo do i=1,n1 do k=1,m if(f(i,k).eq.ff)then f(i,k)=0.0 nf(i)=nf(i)+1 endif enddo enddo return end SUBROUTINE TRANSF(n1,m,f,nf,avf,df,ks) DIMENSION F(N1,M),AVF(N1) DIMENSION DF(N1) DIMENSION nF(N1) if(ks.eq.-1)then goto 30 endif do i=1,n1 avf(i)=0.0 enddo if(ks.eq.0)then goto 5 endif do i=1,n1 df(i)=0.0 enddo5 continue DO 141 I=1,N1 if(nf(i).ne.0) goto 141 do 12 j=1,m 12 AVF(I)=AVF(I)+F(I,J) AVF(I)=AVF(I)/M DO 14 J=1,M 14 F(I,J)=F(I,J)-AVF(I) 141 CONTINUE IF(KS.EQ.0) THEN RETURN ELSE DO 241 I=1,N1 if(nf(i).ne.0) goto 241 DO 22 J=1,M 22 DF(I)=DF(I)+F(I,J)*F(I,J) DF(I)=SQRT(DF(I)/M) DO 24 J=1,M 24 F(I,J)=F(I,J)/DF(I) 241 CONTINUE ENDIF 30 CONTINUE RETURN END SUBROUTINE FORMA(N,M,MNH,F,A) DIMENSION F(N,M),A(MNH,MNH) IF(M-N) 40,50,50 40 DO 44 I=1,MNH DO 44 J=I,MNH A(I,J)=0.0 DO 42 IS=1,N 42 A(I,J)=A(I,J)+F(IS,I)*F(IS,J) A(J,I)=A(I,J) 44 CONTINUE RETURN 50 DO 54 I=1,MNH DO 54 J=I,MNH A(I,J)=0.0 DO 52 JS=1,M 52 A(I,J)=A(I,J)+F(I,JS)*F(J,JS) A(J,I)=A(I,J) 54 CONTINUE RETURN END SUBROUTINE JCB(N,A,S,EPS) DIMENSION A(N,N),S(N,N) DO 30 I=1,N DO 30 J=1,I IF(I-J) 20,10,20 10 S(I,J)=1. GO TO 30 20 S(I,J)=0. S(J,I)=0. 30 CONTINUE G=0. DO 40 I=2,N I1=I-1 DO 40 J=1,I1 40 G=G+2.*A(I,J)*A(I,J) S1=SQRT(G) S2=EPS/FLOAT(N)*S1 S3=S1 L=0 50 S3=S3/FLOAT(N) 60 DO 130 IQ=2,N IQ1=IQ-1 DO 130 IP=1,IQ1 IF(ABS(A(IP,IQ).LT.S3) GOTO 130 L=1 V1=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) DO 110 I=1,N G=A(I,IP)*CT-A(I,IQ)*ST A(I,IQ)=A(I,IP)*ST+A(I,IQ)*CT A(I,IP)=G G=S(I,IP)*CT-S(I,IQ)*ST S(I,IQ)=S(I,IP)*ST+S(I,IQ)*CT 110 S(I,IP)=G DO 120 I=1,N A(IP,I)=A(I,IP) 120 A(IQ,I)=A(I,IQ) G=2.*V2*ST*CT A(IP,IP)=V1*CT*CT+V3*ST*ST-G A(IQ,IQ)=V1*ST*ST+V3*CT*CT+G A(IP,IQ)=(V1-V3)*ST*CT+V2*(CT*CT-ST*ST) A(IQ,IP)=A(IP,IQ) 130 CONTINUE IF(L-1) 150,140,150 140 L=0 GO TO 60 150 IF(S3.GT.S2) GOTO 50 RETURN END SUBROUTINE ARRANG(KV,MNH,A,ER,S) DIMENSION A(MNH,MNH),ER(mnh,4),S(MNH,MNH) TR=0.0 DO 200 I=1,MNH TR=TR+A(I,I) 200 ER(I,1)=A(I,I) MNH1=MNH-1 DO 210 K1=MNH1,1,-1 DO 210 K2=K1,MNH1 IF(ER(K2,1).LT.ER(K2+1,1) THEN C=ER(K2+1,1) ER(K2+1,1)=ER(K2,1) ER(K2,1)=C DO 205 I=1,MNH C=S(I,K2+1) S(I,K2+1)=S(I,K2) S(I,K2)=C 205 CONTINUE ENDIF 210 CONTINUE ER(1,2)=ER(1,1) DO 220 I=2,KV ER(I,2)=ER(I-1,2)+ER(I,1) 220 CONTINUE DO 230 I=1,KV ER(I,3)=ER(I,1)/TR ER(I,4)=ER(I,2)/TR 230 CONTINUE WRITE(6,250) TR 250 FORMAT(/5X,TOTAL SQUARE ERROR=,F20.5) RETURN END SUBROUTINE TCOEFF(KVT,KV,N,M,MNH,S,F,V,evf,tcf,ER) DIMENSION S(MNH,MNH),F(N,M),V(MNH),ER(mnh,4),evf(n,kvt),tcf(m,kvt) DO 360 J=1,KVT C=0. DO 350 I=1,MNH 350 C=C+S(I,J)*S(I,J) C=SQRT(C) DO 160 I=1,MNH s(I,J)=S(I,J)/C 160 evf(I,J)=S(I,J)/C 360 CONTINUE IF(N.LE.M) THEN DO 390 J=1,M DO 370 I=1,N V(I)=F(I,J) F(I,J)=0. 370 CONTINUE do 371 is=1,kvttcf(j,is)=0.371continue DO 380 IS=1,KVT DO 380 I=1,N f(is,j)=F(IS,J)+V(I)*S(I,IS) 380 tcf(j,is)=tcf(J,is)+V(I)*S(I,IS) 390 CONTINUE ELSE DO 410 I=1,N DO 400 J=1,M V(J)=F(I,J) F(I,J)=0. 400 CONTINUE DO 410 JS=1,KVT DO 410 J=1,M f(I,JS)=F(I,JS)+V(J)*S(J,JS) 410 CONTINUE DO 430 JS=1,KVT DO 420 J=1,M tcf(J,JS)=S(J,JS)*SQRT(ER(JS,1) 420 CONTINUE DO 430 I=1,N evf(I,JS)=F(I,JS)/SQRT(ER(JS,1) 430 CONTINUE t=0.0 do 3650 i=1,m3650 t=t+tcf(i,1)*tcf(i,2) write(*,*)tt=0.0 do 3651 i=1,n3651 t=t+evf(i,1)*evf(i,2) write(*,*)t ENDIF RETURN END SUBROUTINE test3(N1,ff,nf,evf,kvt) dimension nf(n1),evf(n1,kvt) do i=1,n1 if(nf(i).ne.0)then do k=1,kvt evf(i,k)=ff enddo endif enddo return end SUBROUTINE OUTER(KV,ER,mnh) DIMENSION ER(mnh,4) WRITE(16,510) 510 FORMAT(/30X,EIGENVALUE AND ANALYSIS ERROR,/) WRITE(16,520) 520 FORMAT(10X,1HH,8X,5HLAMDA,10X,6HSLAMDA,11X,2HPH,12X,3HSPH) WRITE(16,530) (IS,(ER(IS,J),J=1,4),IS=1,KV) 530 FORMAT(1X,I10,4F15.5) WRITE(16,540) 540 FORMAT(/) RETURN END SUBROUTINE OUTVT(KVT,N,M,MNH,S,F,EGVT,ECOF) DIMENSION F(N,M),S(MNH,MNH),EGVT(N,KVT),ECOF(M,KVT) WRITE(16,560) 560 FORMAT(30X,STANDARD EIGENVECTORS,/) WRITE(16,570) (IS,IS=1,KVT) 570 FORMAT(3X,80I7) DO 550 I=1,N IF(M.GE.N) THEN WRITE(16,580) I,(S(I,JS),JS=1,KVT) 580 FORMAT(1X,I3,80F7.3,/) DO 11 JS=1,KVT EGVT(I,JS)=S(I,JS) 11 CONTINUE ELSE WRITE(16,590) I,(F(I,JS),JS=1,KVT) 590 FORMAT(1X,I3,80F7.3) DO 12 JS=1,KVT EGVT(I,JS)=F(I,JS) 12 CONTINUE ENDIF 550 CONTINUE WRITE(16,720) 720 FORMAT(/) WRITE(16,610) 610 FORMAT(30X,TIME-COEFFICENT SERIES OF S. E./) WRITE(16,620) (IS,IS=1,KVT) 620 FORMAT(3X,80I7) DO 600 J=1,M IF(M.GE.N) THEN WRITE(16,630) J,(F(IS,J),IS=1,KVT) 630 FORMAT(1X,I3,80F7.1) DO 13 IS=1,KVT ECOF(J,IS)=F(IS,J) 13 CONTINUE ELSE WRITE(16,640) J,(S(J,IS),IS=1,KVT) 640 FORMAT(1X,I3,80F7.1) DO 14 IS=1,KVT ECOF(J,IS)=S(J,IS) 14 CONTINUE ENDIF 600 CONTINUE RETURN ENDn 运行结果n 编写CTL文件Evf.ctldset d:studyformshixieightevf.grdti