气象统计方法实习BD

实习一：气候场、距平场、均方差场编程如下:parameter(ii=37,jj=17,mon=12,year=4) real var(ii,jj,mon,year),ave(ii,jj,mon),jp(ii,jj,mon,year) real s(ii,jj,mon) integer i,j,iy,m open(5,file=d:ex1h500.dat) open(6,file=d:ex1ave.grd,form=binary)open(7,file=d:ex1jp.grd,form=binary)open(8,file=d:ex1s.grd,form=binary)open(12,file=d:ex1outall.grd,form=binaryopen(9,file=d:ex1ave.txt)open(10,file=d:ex1jp.txt) open(11,file=d:ex1s.txt) ！读数据 DO iy=1,4 do m=1,12!ccc read h500 read(5,1000) read(5,2000) (var(i,j,m,iy),i=1,ii),j=1,jj) enddoenddo！计算气候场 do j=1,jj do i=1,ii do m=1,12 ave(i,j,m)=var(i,j,m,1)+var(i,j,m,2)+var(i,j,m,3)+var(i,j,m,4) ave(i,j,m)=ave(i,j,m)/4.0 enddo enddo enddo！计算距平场 do iy=1,4 do m=1,12 do j=1,jj do i=1,ii jp(i,j,m,iy)=var(i,j,m,iy)-ave(i,j,m) enddo enddo enddo enddo！计算均方差场 do j=1,jj do i=1,ii do m=1,12 s(i,j,m)=jp(i,j,m,1)*jp(i,j,m,1)+jp(i,j,m,2)*jp(i,j,m,2)+jp(i,j /,m,3)*jp(i,j,m,3)+jp(i,j,m,4)*jp(i,j,m,4) s(i,j,m)=s(i,j,m)/4.0 s(i,j,m)=sqrt(s(i,j,m) enddo enddo enddo do iy=1,4 do m=1,12 write(6)(ave(i,j,m),i=1,ii),j=1,jj) write(7)(jp(i,j,m,iy),i=1,ii),j=1,jj) write(8)(s(i,j,m),i=1,ii),j=1,jj) write(9,2000)(ave(i,j,m),i=1,ii),j=1,jj) write(10,2000)(jp(i,j,m,iy),i=1,ii),j=1,jj) write(11,2000)(s(i,j,m),i=1,ii),j=1,jj) write(12)(ave(i,j,m),i=1,ii),j=1,jj) write(12)(jp(i,j,m,iy),i=1,ii),j=1,jj) write(12)(s(i,j,m),i=1,ii),j=1,jj) enddoenddo1000 format(2i7)2000 format(37f8.1) close(5) close(6) close(7) close(8) close(9) close(10) close(11) close(12) end给ave配的ctl文件：dset d:ex1ave.grdundef -9.99E+33title NCEP/NCAR REANALYSIS PROJECTxdef 37 linear 60.000 2.500ydef 17 linear 0.000 2.500zdef 1 levels 500tdef 12 linear JAN1982 12movars 1ave 1 99 H500endvars给ave配的gs文件：reinitopen d:ex1ave.ctlenable print d:ex1ave.gmfmon=1while(mon=12)set t mond avedraw title qihouchang of mon printcmon=mon+1endwhiledisable print;气候场图：一月份高度的气候场呈现南高北低的状态，陆地上的高度场比较稀疏，而在西太平洋上高度场比较密集。八月份高度的气候场呈现东高西低的状态，在我国东北部以北以及印度东北部出现低压中心，而在赤道西太平洋地区出现高压中心。35N以北高度分布很密集，而35N以南比较稀疏。给jp配的ctl文件：dset d:ex1jp.grdundef -9.99E+33title NCEP/NCAR REANALYSIS PROJECTxdef 37 linear 60.000 2.500ydef 17 linear 0.000 2.500zdef 1 levels 500tdef 48 linear JAN1982 1movars 1jp 1 99 H500endvars给jp配的gs文件：reinitopen d:ex1jp.ctlenable print d:ex1jp.gmfyear=1982while(year=1985)mon=1while(mon=12)set t mond jpdraw title jupingchang of year.monprintcmon=mon+1endwhileyear=year+1endwhiledisable print;距平场图：1983年6月距平场在日本地区出现低压中心，在我国南部出现高压中心,在亚洲西北部也有高压中心。赤道至25N间以及25N-40N，60E-100E间基本都是正距平，而在25N-40N，100E-150E间基本都是负距平。1984年7月距平场在亚洲大陆西部、日本地区、赤道西太平洋地区形成低压中心，太平洋西北部形成高压中心。给s配的ctl文件：dset d:ex1s.grdundef -9.99E+33title NCEP/NCAR REANALYSIS PROJECTxdef 37 linear 60.000 2.500ydef 17 linear 0.000 2.500zdef 1 levels 500tdef 12 linear JAN1982 12movars 1s 1 99 H500endvars给s配的gs文件：reinitopen d:ex1s.ctlenable print d:ex1s.gmfmon=1while(monabs(zxgxs1(t1)t1=i if(abs(zxgxs2(i)abs(zxgxs2(t2)t2=i if(abs(lhxgxs(i)abs(lhxgxs(t3)t3=ienddoprint *,中国1970-1989年年平均气温自相关系数绝对值最大的滞后时间长度为：,t1,zxgxs1(t1)print *,中国1970-1989年冬季平均气温自相关系数绝对值最大的滞后时间长度为：,t2,zxgxs2(t2)print *,中国1970-1989年平均气温和冬季平均气温之间落后相关系数绝对值最大的滞后时间长度为：,t3,lhxgxs(t3)end程序运行如下：中国1970-1989年年平均气温自相关系数绝对值最大的滞后时间长度是7，自相关系数是-0.3724，说明在这些年的数据中，滞后7年的序列与原序列的相关系数绝对值最大，呈反相关。中国1970-1989年冬季平均气温自相关系数绝对值最大的滞后时间长度是4，自相关系数是-0.3678，说明在这些年的数据中，滞后4年的序列与原序列的相关系数绝对值最大，呈反相关。中国1970-1989年平均气温和冬季平均气温之间落后相关系数绝对值最大的滞后时间长度是3，自相关系数是-0.4066，说明在这些年的数据中，滞后3年的冬季平均气温序列与原平均气温序列的相关系数绝对值最大，呈反相关。实习三：一元线性回归Fortran程序如下：program ex3integer,parameter:n=17 integer i,jreal:suma1=0,suma2=0,avea1,avea2,s1=0,s2=0,b0,b,x=14.5,yinteger a1(n)real a2(n),a3(n)data a1/0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16/data a2/30.0,29.1,28.4,28.1,28.0,27.7,27.5,27.2,27.0,26.8,26.5,26.3,26.1,25.7,25.3,24.8,24.0/data a3/0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/open(5,file=d:3ys.txt)open(15,file=d:3ys.grd,form=binary)open(6,file=d:3hs.txt)open(16,file=d:3hs.grd,form=binary)!求平均值do i=1,n suma1=suma1+a1(i) suma2=suma2+a2(i)enddoavea1=suma1/navea2=suma2/n!求b,b0值do i=1,n s1=s1+a1(i)*a2(i) s2=s2+a1(i)*a1(i)enddob=(s1-suma1*suma2/n)/(s2-suma1*suma1/n)b0=avea2-avea1*by=b0+b*xprint*,y=,b0,+,b,xprint*print*,y!求回归数据do i=1,na3(i)=b0+b*a1(i)enddowrite(5,1000)(a1(i),i=1,n)write(5,2000)(a2(i),i=1,n)write(15)(a1(i),i=1,n)write(15)(a2(i),i=1,n)write(6,1000)(a1(i),i=1,n)write(6,2000)(a3(i),i=1,n)write(16)(a1(i),i=1,n)write(16)(a3(i),i=1,n)1000 format(i7)2000 format(f8.1)end程序运行：所给数据做出的线性回归曲线的斜率是-0.3012，截距是29.3804，说明所给数据y随着x递减。实习四：滑动平均Fortran程序如下：PROGRAM MAinteger iinteger,parameter:n=85,ih=11,nyear=1922integer:yr(n)=0 real X(n),X1(n)real:s(75)=0C*C* N: SAMPLE SIZE OF THE TIME SERIES *C* IH: MOVING LENGTH *C * NYEAR: FIRST YEAR OF THE SERIES *C* X(N): OROGINAL TIME SERIES *C * X1(N-IH+1): MOVED SERIES *C*OPEN(2,FILE=d:4MA.DAT)open(12,file=d:4ma.grd,form=binary) OPEN(3,FILE=d:4hd.DAT)open(13,file=d:4hd.grd,form=binary) !年份 do i=1,nyr(i)=1922-1+ienddo!读入数据READ(2,*)(X(I),I=1,N) !计算滑动平均do i=1,11 s(1)=s(1)+x(i)enddox1(6)=s(1)/ihdo i=7,80 s(i-5)=s(i-6)+x(i+ih-6)-x(i-6) x1(i)=s(i-5)/ihenddo!写入数据 write(3,(1x,i5,1X,f5.1,1X,f5.1)(yr(i),x(i),x1(i),i=1,n) write(12)(x(i),i=1,n)write(13)(x1(i),i=1,n) close(2)close(3)close(12)close(13)end原数据和滑动后数据的图形：由图可知，所给数据在1922-1955年之间呈波动下降趋势，在1955-1968年呈波动上升趋势，上升幅度较大，而1968-2006年之间大致在同一水平上波动，没有升降趋势。实习五：eofFortran程序：C$large PROGRAM EOFPWC THIS PROGRAM USES EOF FOR ANALYSING TIME SERIESC OF METEOROLOGICAL FIELDC M:LENTH OF TIME SERIESC N:NUMBER OF GRID-POINTSC KS=-1:SELF; KS=0:DEPATURE; KS=1:STANDERDLIZED DEPATUREC KV:NUMBER OF EIGENVALUES WILL BE OUTPUTC KVT:NUMBER OF EIGENVECTORS AND TIME SERIES WILL BE OUTPUTC MNH=MIN(M,N)C Evf=EIGENVACTORS, tcF=TIME COEFFICIENTS FOR EGVT.C ER(KV,1)=LAMDA,LAMDA EIGENVALUEC ER(KV,2)=ACCUMULATE LAMDAC ER(KV,3)=THE SUM OF COMPONENTS VECTORS PROJECTED ONTOc EIGENVACTOR.C ER(KV,4)=ACCUMULATE ER(KV,3)C PARAMETER(M=516,N=216,MNH=216,KS=-1,KV=10,KVT=2) PARAMETER(ff=-999.0,nx=18,ny=12)C 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)CCCCCCCCCCCCCCCCINPUT DATA open(11,file=d:5sstpx.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)CCCCCCCCCCCCCCCCINPUT DATA CCCCCCCCCCCCCCCCCCCccccccccccccccccccccccccccccccccccccc 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:5evf.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:5tcf.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)c CALL OUTER(KV,ER,mnh) open(21,file=d:5dats.dat) write(21,106)(er(iiii,3),iiii=1,kv) close(21)c CALL OUTVT(KVT,N,M,MNH,S,F,Evf,tcF) 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)C THIS SUBROUTINE PROVIDES INITIAL F BY KS and kv. 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 enddocccccccccccccccccc5 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)C THIS SUBROUTINE FORMS A BY F 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)C THIS SUBROUTINE COMPUTS EIGENVALUES AND EIGENVECTORS OF A 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)C THIS SUBROUTINE PROVIDES A SERIES OF EIGENVALUESC FROM MAX TO MIN 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)C THIS SUBROUTINE PROVIDES STANDARD EIGENVECTORS (M.GE.N,SAVED IN S;C M.LT.N,SAVED IN F) AND ITS TIME COEFFICENTS SERIES (M.GE.N,C SAVED IN F; M.LT.N,SAVED IN S) 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 CONTINUEccccccccccccccccccccccccccccccccccccccccccc t=0.0c do 365 i=1,mnhc365 t=t+s(i,1)*s(i,2)c write(*,*)tcccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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 CONTINUEccccccccccccccccccccccccccccccccccccccccccccccccccccc 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)c this subroutine sent undefine value ff to evf 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)C THIS SUBROUTINE PRINTS ARRAY ERC ER(KV,1) FOR SEQUENCE OF EIGENVALUE FROM BIG TO SMALLC ER(KV,2) FOR EIGENVALUE FROM BIG TO SMALLC ER(KV,3) FOR SMALL LO=(LAMDA/TOTAL VARIANCE)C ER(KV,4) FOR BIG LO=SUM OF SMALL LO) 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)C THIS SUBROUTINE PRINTS STANDA