短期气候预测实习程序总结.doc

本人为南京信息工程大学大气科学系学生，在我大四上时，不幸选修课选了短期气候预测实习课程。其烦人程度超乎我的想象，一边在准备考研和找工作，一边还要花心思完成实习内容，当时根本喘不过气来。为让学弟学妹们不重蹈覆辙，在此共享出我所用的程序，给你们一些帮助。不过需要明确一点，这次实习对于个人编程水平的提高帮助很大，学有余力的同学，应仅把本文当做参考，理解基础上使用，而非不劳而获的资本。以下为正文：具体实习要求参见课本，在此不赘述实习一：大气环流状况的表征 program EX1real a(144,73,12,65),ave1(144,73),ave7(144,73),asum(144,73)real dev(144,73,65),latave(73,12,65),latsum(73,12,65)real latdev(144,73,12,65)open(2,file=d:1hgt500.grd,form=binary) !补充正确路径open(4,file=d:1ave7.grd,form=binary) !补充数据输出路径open(5,file=d:1dev.grd,form=binary)open(6,file=d:1latdev.grd,form=binary) do it=1,65do imo=1,12do j=1,73do i=1,144read(2)a(i,j,imo,it)enddo;enddo;enddo;enddocccccccccccccc 请完成以下的程序!月时间平均 ave7imo=7do i=1,144do j=1,73do it=1,65asum(i,j)=asum(i,j)+a(i,j,imo,it)enddoave7(i,j)=asum(i,j)/65.0enddoenddo!7月距平 deviationdo i=1,144do j=1,73do it=1,65dev(i,j,it)=a(i,j,imo,it)-ave7(i,j)enddoenddoenddo!纬圈平均 latitude averagedo it=1,65do imo=1,12do j=1,73do i=1,144latsum(j,imo,it)=latsum(j,imo,it)+a(i,j,imo,it)enddolatave(j,imo,it)=latsum(j,imo,it)/144.0enddoenddoenddo!纬向偏差 latitude deviationdo it=1,65do imo=1,12do j=1,73do i=1,144latdev(i,j,imo,it)=a(i,j,imo,it)-latave(j,imo,it)enddoenddoenddoenddo!写数据write(4) (ave7(i,j),i=1,144),j=1,73)write(5) (dev(i,j,it),i=1,144),j=1,73),it=1,65)write(6) (latdev(i,j,imo,it),i=1,144),j=1,73),imo=1,12),it=1,65)close(4)close(5)close(6)End实习二：大气环流分型 PROGRAM EOFC THIS PROGRAM USES EOF FOR ANALYSING TIME SERIESC OF METEOROLOGICAL FIELDC M:LENTH OF TIME SERIES ! m:时间序列长度C N:NUMBER OF GRID-POINTS ! n:格点数C 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 EGVT=EIGENVACTORS, ECOF=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=61,N=41*21,MNH=61,KS=1,KV=8,KVT=8,pi=3.1415926)C DIMENSION F(N,M),A(MNH,MNH),S(MNH,MNH),ER(MNH,4), * DF(N),V(MNH),AVF(N),EGVT(N,KVT),ECOF(M,KVT)dimension hh(144,73,12,61),h(41,21,61) open(10,file=d:2hgt500.grd,form=binary,status=old)open(20,file=d:2egvt7.grd,form=binary)open(30,file=d:2t7.grd,form=binary)open(16,file=d:2eof7.txt)cccccccccccccccccc读数据do it=1,61do k=1,12do j=1,73do i=1,144read(10)hh(i,j,k,it)enddo;enddo;enddo;enddo write(*,*)read data ok!裁剪区域k=1do it=1,61do j=1,21do i=1,41h(i,j,it)=hh(i+16,j+44,k,it)enddoenddoenddowrite(*,*)data narrowed!二维空间场变一维数组，注意按照grads的XY顺序，!因为最后EIGENVACTORS文件里面直接按照该格式存的,不再经过这一步变化do it=1,Mii=1do j=1,21do i=1,41F(ii,it)=h(i,j,it)ii=ii+1enddoenddoenddoCCCCCCCCCCCCCCCCINPUT DATA CCCCCCCCCCCCCCCCCCC CALL TRANSF(N,M,F,AVF,DF,KS)write(*,*)ok program 1 CALL FORMA(N,M,MNH,F,A)write(*,*)ok program 2 CALL JCB(MNH,A,S,0.00001)write(*,*)ok program 3 CALL ARRANG(KV,MNH,A,ER,S)write(*,*)ok program 4 CALL TCOEFF(KVT,KV,N,M,MNH,S,F,V,ER)write(*,*)ok program 5 CALL OUTER(KV,ER,MNH)write(*,*)ok program 6 CALL OUTVT(KVT,N,M,MNH,S,F,EGVT,ECOF)write(*,*)ok program 7ccccccccccccc存储数据do j=1,mdo i=1,kvtwrite(30)ecof(j,i)enddo;enddo do it=1,kvtdo j=1,nwrite(20)egvt(j,it)enddo;enddo write(*,*)ok 8cccccccccccc ENDccccccccccccccccccccccccc子程序 SUBROUTINE TRANSF(N,M,F,AVF,DF,KS)C THIS SUBROUTINE PROVIDES INITIAL F BY KS DIMENSION F(N,M),AVF(N),DF(N) DO 5 I=1,N AVF(I)=0.0 5 DF(I)=0.0 IF(KS) 30,10,10 10 DO 14 I=1,N DO 12 J=1,M 12 AVF(I)=AVF(I)+F(I,J) AVF(I)=AVF(I)/M DO 14 J=1,M F(I,J)=F(I,J)-AVF(I) 14 CONTINUE IF(KS.EQ.0) THEN RETURN ELSE DO 24 I=1,N 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 F(I,J)=F(I,J)/DF(I) 24 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 standard 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(*,250) TR 250 FORMAT(/5X,TOTAL SQUARE ERROR=,F20.5) RETURN END SUBROUTINE TCOEFF(KVT,KV,N,M,MNH,S,F,V,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) 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 380 IS=1,KVT DO 380 I=1,N 380 F(IS,J)=F(IS,J)+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 S(J,JS)=S(J,JS)*SQRT(ER(JS,1) 420 CONTINUE DO 430 I=1,N F(I,JS)=F(I,JS)/SQRT(ER(JS,1) 430 CONTINUE ENDIF 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(/10X,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 STANDARD EIGENVECTORSC AND ITS TIME-COEFFICENT SERIES DIMENSION F(N,M),S(MNH,MNH),EGVT(N,KVT),ECOF(M,KVT) WRITE(16,560) 560 FORMAT(10X,STANDARD EIGENVECTORS) WRITE(16,570) (IS,IS=1,KVT) 570 FORMAT(3X,10i7) DO 550 I=1,N IF(M.GE.N) THEN WRITE(16,580) I,(S(I,JS),JS=1,KVT) 580 FORMAT(1X,I3,10F7.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,I5,10F7.3) DO 12 JS=1,KVT EGVT(I,JS)=F(I,JS) 12 CONTINUE ENDIF 550 CONTINUEC WRITE(16,590) I,(F(I,JS),JS=1,KVT)! WRITE(20)(F(I,JS),i=1,n),JS=1,KVT) WRITE(16,720) 720 FORMAT(/) WRITE(16,610) 610 FORMAT(10X,TIME-COEFFICENT SERIES OF S. E.) WRITE(16,620) (IS,IS=1,KVT) 620 FORMAT(3X,5i12) DO 600 J=1,M IF(M.GE.N) THEN WRITE(16,630) J,(f(is,j),is=1,kvt) 630 FORMAT(1X,I3,5F12.3) 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,10F12.3) DO 14 IS=1,KVT ECOF(J,IS)=S(J,IS) 14 CONTINUE ENDIF 600 CONTINUEC WRITE(30)(S(J,IS),j=1,m),IS=1,KVT) RETURN END实习三：大气遥相关program ex3parameter nt=60real t(160,nt),eu(63),h(144,73,63),hsum(144,73),have(144,73),euave,rh(144,73),r(144,73),r2(160)real a(144,73,12,63),ave1(144,73),ave7(144,73),rup(144,73),rh_2(144,73),reu_2,eu2,reu,tave(160),tsum(160)real rh2(160),rh2_2(160),rup2(160)real lat(160),lon(160),timcharacter*8 id(160)integer nlev,nflag!读数据open(3,file=D:3t1601.txt)open(2,file=d:3hgt500.grd,form=binary) open(4,file=D:3rheu.grd,form=binary) !h和eu的相关open(5,file=d:3eu.grd,form=binary) !eu指数open(6,file=D:3rteu.grd,status=replace,form=binary) !t和eu的相关open(7,file=D:3lat_lon.txt)read(3,*)(t(i,j),i=1,160),j=1,nt)!读160站温度do i=1,160do j=1,ntt(i,j)=t(i,j)/10.0enddoenddodo i=1,160read(7,*)lat(i),lon(i)!读经纬度enddo do it=1,63do imo=1,12do j=1,73do i=1,144read(2)a(i,j,imo,it)!读高度场enddo;enddo;enddo;enddo!计算EU指数do it=1,63eu(it)=-0.25*a(9,59,1,it)+0.5*a(31,59,1,it)-0.25*a(59,53,1,it)enddo!计算EU指数与高度场的相关系数!h-高度场 hsum-高度场和 have-高度场平均值 eu-EU指数,eusum,euave类似!1,提取1月份高度场do it=1,63do j=1,73do i=1,144h(i,j,it)=a(i,j,1,it)enddoenddoenddo!2,计算高度场和EU指数的平均值do j=1,73do i=1,144do it=1,63hsum(i,j)=hsum(i,j)+h(i,j,it)enddohave(i,j)=hsum(i,j)/63enddoenddodo it=1,63eusum=eusum+eu(it)enddoeuave=eusum/63!3,计算相关系数各部：分子、分母、分母(对照相关系数公式)!rup-分子 rh-分母h reu-分母eu r-相关系数do j=1,73do i=1,144reu_2=0do it=1,63rup(i,j)=rup(i,j)+(eu(it)-euave)*(h(i,j,it)-have(i,j)rh_2(i,j)=rh_2(i,j)+(h(i,j,it)-have(i,j)*2reu_2=reu_2+(eu(it)-euave)*2enddorh(i,j)=sqrt(rh_2(i,j)reu=sqrt(reu_2)enddoenddoprint*,reudo j=1,73do i=1,144r(i,j)=rup(i,j)/(rh(i,j)*reu)enddoenddo!计算EU指数和气温的相关系数!1,计算温度场的平均值do i=1,160do it=1,nttsum(i)=tsum(i)+t(i,it)enddotave(i)=tsum(i)/ntenddo!print*, (tave(i),i=1,160)!2,计算相关系数各部：分子、分母、分母(对照相关系数公式)!rup2-分子 rh2-分母h reu-分母eu r2-相关系数eusum=0do it=1,nteusum=eusum+eu(it+3)!之所以加3，是因为在分析资料和观测资料起始年份差3年enddoeuave=eusum/ntreu_2=0do it=1,ntreu_2=reu_2+(eu(it+3)-euave)*2enddoreu=sqrt(reu_2)do i=1,160do it=1,ntrup2(i)=rup2(i)+(eu(it+3)-euave)*(t(i,it)-tave(i)rh2_2(i)=rh2_2(i)+(t(i,it)-tave(i)*2enddorh2(i)=sqrt(rh2_2(i)enddodo i=1,160r2(i)=rup2(i)/(rh2(i)*reu)enddo!计算完毕，写数据!写站点数据do j=1,160id(j)=char(j)tim=0.0nlev=1nflag=1write(6)id(j),lat(j),lon(j),tim,nlev,nflag,r2(j)enddotim=0.0nlev=0nflag=1write(6)id(j-1),lat(j-1),lon(j-1),tim,nlev,nflagprint*, (r2(j),j=1,160)write(5)(eu(it),it=1,nt)!写EU指数write(4)(r(i,j),i=1,144),j=1,73)!写h和EU的相关系数close(2)close(3)close(4)close(5)close(6)close(7)end实习四：预测因子的选择合成分析方法（1） program EX4real a6(160,63),a7(160,62),a8(160,62),lat(160),a(160,60), & lon(160),rap(160,60),ddi(60,3),rsum(160),r(160,60),rave(160)real num,rapave(160,3),h(144,73,12,65),hw(144,73,60)real hwavereal hwa60(144,73,60)real hwa(144,73,3),hwat(144,73,60,3)real rapt(160,60,3),rap1(160,22),rap2(160,19),rap3(160,19)real t(144,73,3)character*8 id(160)open(3,file=D:4r1606.txt) !请修改路径 open(4,file=D:4r1607.txt)open(5,file=D:4r1608.txt)open(6,file=D:4lat_lon.txt)open(7,file=D:4ddi)open(8,file=D:4rapave.grd,form=binary)open(9,file=D:4hgt500.grd,form=binary)open(10,file=D:4hwa.grd,form=binary)open(11,file=D:4t.grd,form=binary)ccccccccccccccc 读数据（经纬度、160站降水、雨型） read(3,*)(a6(i,j),i=1,160),j=1,63) read(4,*)(a7(i,j),i=1,160),j=1,62) read(5,*)(a8(i,j),i=1,160),j=1,62)do i=1,160read(6,*)lat(i),lon(i)enddodo it=1,60read(7,*)(ddi(it,j),j=1,3)enddo!print*,(ddi(it,j),j=1,3),it=1,60)ccccccccccccccc 编程求合成!读高度场do it=1,65do imo=1,12do j=1,73do i=1,144read(9)h(i,j,imo,it)!读高度场enddo;enddo;enddo;enddo!30年降水平均值 do i=1,160 rsum(i)=0.0do it=21,50rsum(i)=rsum(i)+a6(i,it)+a7(i,it)+a8(i,it)enddorave(i)=rsum(i)/30.0enddoccccccccccc 请补充!求夏季降水(3个月相加)do i=1,160do it=1,60a(i,it)=a6(i,it)+a7(i,it)+a8(i,it)enddoenddo!计算百分率do i=1,160do it=1,60rap(i,it)=(a(i,it)-rave(i)/rave(i)enddoenddo!挑雨型,并求各个雨型的百分率平均值do i=1,160do j=1,3num=0.0do it=1,60num=num+ddi(it,j)rapave(i,j)=rapave(i,j)+rap(i,it)*ddi(it,j)enddorapave(i,j)=rapave(i,j)/numenddoenddo! 3类雨型各自前期的冬季高度场距平!前期，该夏季之前的那个冬季!求冬季平均高度场hw(height of winter),时间往后移3年，1948+3=1951do it=1,60do imo=1,12do j=1,73do i=1,144hw(i,j,it)=(h(i,j,1,it+3)+h(i,j,2,it+3)+h(i,j,12,it+3-1)/3.0enddo;enddo;enddo;enddo!求3类雨型前期冬季距平hwa1;hwa2;hwa3 height of winter abnormal!先求60年的距平hwa60do j=1,73do i=1,144hwa=0call cal_ave(hw(i,j,:),60,hwave)do it=1,60hwa60(i,j,it)=hw(i,j,it)-hwaveenddoenddoenddo!再挑相应雨型年份的距平do j=1,73do i=1,144do k=1,3num=0.0do it=1,60num=num+ddi(it,k)hwa(i,j,k)=hwa(i,j,k)+hwa60(i,j,it)*ddi(it,k)enddohwa(i,j,k)=hwa(i,j,k)/numenddoenddoenddo!对各个雨型的距平进行t检验!t检验之前的准备：制作样本数组 hwat:Height of Winter Abnormal for Testdo j=1,73do i=1,144do k=1,3do it=1,60hwat(i,j,it,k)=hwa60(i,j,it)*ddi(it,k)enddo;enddo;enddo;enddodo j=1,73do i=1,144do k=1,3call forward_push(hwat(i,j,:,k),60)enddo;enddo;enddo!开始T检验t=0do j=1,73do i=1,144do k=1,3if(k=1) kk=22else kk=19call t_test(hwat(i,j,:,k),hwa60(i,j,:),kk,60,t(i,j,k)enddoenddoenddo!写前期冬季距平数据write(10)(hwa(i,j,k),i=1,144),j=1,73),k=1,3)write(10)(t(i,j,k),i=1,144),j=1,73),k=1,3)close(10)ccccccccccccccccccc写站点数据!写3类雨型百分率合成值rapave rain abnormal percentige averagedo j=1,160id(j)=char(j)tim=0.0nlev=1nflag=1write(8)id(j),lat(j),lon(j),tim,nlev,nflag,(rapave(j,i),i=1,3)enddotim=0.0nlev=0nflag=1write(8)id(j-1),lat(j-1),lon(j-1),tim,nlev,nflag!写t检验结果end!以下为子程序subroutine t_test(sam,set,nsam,nset,t)!t检验子程序!sam样本数组 set总体数组 nsam样本量 nset总体量!t检验值real sam(nsam),set(nset),t,samave,setave,samsinteger nsam,nsetsetave=0samave=0sams=0call cal_ave(sam,nsam,samave)call cal_ave(set,nset,setave)call cal_s(sam,nsam,sams)t=(samave-setave)/sams)*sqrt(real(nsam)end subroutinesubroutine cal_ave(a,na,ave)!求平均值子程序!a数组 na数组量 ave平均值real a(na),ave,asuminteger naave=0asum=0do i=1,naasum=asum+a(i)enddoave=asum/real(na)endsubroutinesubroutine cal_s(a,na,s)!计算标准差子程序!a数组 na数组量 s标准差real a(na),s,aveinteger naave=0s=0call cal_ave(a,na,ave)do i=1,nas=s+(a(i)-ave)*2enddos=s/real(na)s=sqrt(s)endsubroutinesubroutine forward_push(a,na)!数组前缩子程序!效果，将数组中0值后移，其他值前移integer nareal a(na),b(na)b=0j=1do i=1,naif(a(i)/=0) thenb(j)=a(i)j=j+1endifenddoa=0do i=1,naa(i)=b(i)enddoendsubroutine实习五：预测因子的选择合成分析方法（2） program EX5real a(180,89,12,64),sst12(180,89,60),ssta(180,89,60)real sstave,ddi(60,3),sstaa(180,89,3),numreal sstat(180,89,60,3),t(180,89,3)real sas(180,89,3)integer k1,k2,kk1,kk2open(3,file=D:5sst.grb,form=binary) !修改open(7,file=D:5ddi)open(8,file=D:5sst12.grb,form=binary)ccccccccccccccc 读数据（经纬度、160站降水、雨型） do it=1,64do imo=1,12do j=1,89do i=1,180read(3,end=10)a(i,j,imo,it)enddo;enddo;enddo;enddo10 continuedo it=1,60read(7,*)(ddi(it,j),j=1,3)enddoccccccccccccccc 编程求合成!计算前期12月距平合成图!先提取60年所有前期12月距平 sst12:年份往后移do j=1,89do i=1,180do it=1,60sst12(i,j,it)=a(i,j,12,it+4-1