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

短期气候预测实习程序总结短期气候预测实习程序总结短期气候预测实习程序总结xxx公司短期气候预测实习程序总结文件编号：文件日期：修订次数：第1.0次更改批准审核制定方案设计，管理制度本人为南京信息工程大学大气科学系学生，在我大四上时，不幸选修课选了短期气候预测实习课程。其烦人程度超乎我的想象，一边在准备考研和找工作，一边还要花心思完成实习内容，当时根本喘不过气来。为让学弟学妹们不重蹈覆辙，在此共享出我所用的程序，给你们一些帮助。不过需要明确一点，这次实习对于个人编程水平的提高帮助很大，学有余力的同学，应仅把本文当做参考，理解基础上使用，而非不劳而获的资本。以下为正文：具体实习要求参见课本，在此不赘述实习一：大气环流状况的表征programEX1 reala(144,73,12,65),ave1(144,73),ave7(144,73),asum(144,73) realdev(144,73,65),latave(73,12,65),latsum(73,12,65) reallatdev(144,73,12,65) open(2,file='d:\1\',form='binary')!补充正确路径 open(4,file='d:\1\',form='binary')!补充数据输出路径 open(5,file='d:\1\',form='binary') open(6,file='d:\1\',form='binary')doit=1,65 doimo=1,12 doj=1,73 doi=1,144 read(2)a(i,j,imo,it) enddo;enddo;enddo;enddocccccccccccccc请完成以下的程序 !月时间平均ave7 imo=7 doi=1,144 doj=1,73 doit=1,65 asum(i,j)=asum(i,j)+a(i,j,imo,it) enddo ave7(i,j)=asum(i,j)/ enddo enddo !7月距平deviation doi=1,144 doj=1,73 doit=1,65 dev(i,j,it)=a(i,j,imo,it)-ave7(i,j) enddo enddo enddo !纬圈平均latitudeaverage doit=1,65 doimo=1,12 doj=1,73 doi=1,144 latsum(j,imo,it)=latsum(j,imo,it)+a(i,j,imo,it) enddo latave(j,imo,it)=latsum(j,imo,it)/ enddo enddo enddo !纬向偏差latitudedeviation doit=1,65 doimo=1,12 doj=1,73 doi=1,144 latdev(i,j,imo,it)=a(i,j,imo,it)-latave(j,imo,it) enddo enddo enddo enddo !写数据 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实习二：大气环流分型PROGRAMEOFCTHISPROGRAMUSESEOFFORANALYSINGTIMESERIESCOFMETEOROLOGICALFIELDCM:LENTHOFTIMESERIES!!!!!!!!!!m:时间序列长度CN:NUMBEROFGRID-POINTS!!!!!!!!!!n:格点数CKS=-1:SELF;KS=0:DEPATURE;KS=1:STANDERDLIZEDDEPATURECKV:NUMBEROFEIGENVALUESWILLBEOUTPUTCKVT:NUMBEROFEIGENVECTORSANDTIMESERIESWILLBEOUTPUTCMNH=MIN(M,N)CEGVT=EIGENVACTORS,ECOF=TIMECOEFFICIENTSFOREGVT.CER(KV,1)=LAMDA,LAMDAEIGENVALUECER(KV,2)=ACCUMULATELAMDACER(KV,3)=THESUMOFCOMPONENTSVECTORSPROJECTEDONTOcEIGENVACTOR.CER(KV,4)=ACCUMULATEER(KV,3)CPARAMETER(M=61,N=41*21,MNH=61,KS=1,KV=8,KVT=8,pi=CDIMENSIONF(N,M),A(MNH,MNH),S(MNH,MNH),ER(MNH,4),*DF(N),V(MNH),AVF(N),EGVT(N,KVT),ECOF(M,KVT) dimensionhh(144,73,12,61),h(41,21,61)open(10,file='d:\2\',form='binary',status='old') open(20,file='d:\2\',form='binary') open(30,file='d:\2\',form='binary') open(16,file='d:\2\')cccccccccccccccccc读数据 doit=1,61 dok=1,12 doj=1,73 doi=1,144 read(10)hh(i,j,k,it) enddo;enddo;enddo;enddowrite(*,*)'readdataok' !裁剪区域 k=1 doit=1,61 doj=1,21 doi=1,41 h(i,j,it)=hh(i+16,j+44,k,it) enddo enddo enddo write(*,*)'datanarrowed' !二维空间场变一维数组，注意按照grads的XY顺序， !因为最后EIGENVACTORS文件里面直接按照该格式存的,不再经过这一步变化 doit=1,M ii=1 doj=1,21 doi=1,41 F(ii,it)=h(i,j,it) ii=ii+1 enddo enddo enddoCCCCCCCCCCCCCCCCINPUTDATACCCCCCCCCCCCCCCCCCCCALLTRANSF(N,M,F,AVF,DF,KS) write(*,*)'okprogram1'CALLFORMA(N,M,MNH,F,A) write(*,*)'okprogram2'CALLJCB(MNH,A,S, write(*,*)'okprogram3'CALLARRANG(KV,MNH,A,ER,S) write(*,*)'okprogram4'CALLTCOEFF(KVT,KV,N,M,MNH,S,F,V,ER) write(*,*)'okprogram5'CALLOUTER(KV,ER,MNH) write(*,*)'okprogram6'CALLOUTVT(KVT,N,M,MNH,S,F,EGVT,ECOF) write(*,*)'okprogram7'ccccccccccccc存储数据 doj=1,m doi=1,kvt write(30)ecof(j,i) enddo;enddodoit=1,kvt doj=1,n write(20)egvt(j,it) enddo;enddo write(*,*)'ok8'ccccccccccccENDccccccccccccccccccccccccc子程序SUBROUTINETRANSF(N,M,F,AVF,DF,KS)CTHISSUBROUTINEPROVIDESINITIALFBYKSDIMENSIONF(N,M),AVF(N),DF(N)DO5I=1,NAVF(I)=5DF(I)=IF(KS)30,10,1010DO14I=1,NDO12J=1,M12AVF(I)=AVF(I)+F(I,J)AVF(I)=AVF(I)/MDO14J=1,MF(I,J)=F(I,J)-AVF(I)14CONTINUEIFTHENRETURNELSEDO24I=1,NDO22J=1,M22DF(I)=DF(I)+F(I,J)*F(I,J)DF(I)=SQRT(DF(I)/M)DO24J=1,MF(I,J)=F(I,J)/DF(I)24CONTINUEENDIF30CONTINUERETURNENDSUBROUTINEFORMA(N,M,MNH,F,A)CTHISSUBROUTINEFORMSABYFDIMENSIONF(N,M),A(MNH,MNH)IF(M-N)40,50,5040DO44I=1,MNHDO44J=I,MNHA(I,J)=DO42IS=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)=DO52JS=1,M52A(I,J)=A(I,J)+F(I,JS)*F(J,JS)A(J,I)=A(I,J)54CONTINUERETURNENDSUBROUTINEJCB(N,A,S,EPS)CTHISSUBROUTINECOMPUTSEIGENVALUESANDstandardEIGENVECTORSOFADIMENSIONA(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)).GOTO130L=1V1=A(IP,IP)V2=A(IP,IQ)V3=A(IQ,IQ)U=*(V1-V3)IFG=1.IF(ABS(U).G=-SIGN(1.,U)*V2/SQRT(V2*V2+U*U)ST=G/SQRT(2.*(1.+SQRT*G)))CT=SQRT*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=0GOTO60150IFGOTO50RETURNENDSUBROUTINEARRANG(KV,MNH,A,ER,S)CTHISSUBROUTINEPROVIDESASERIESOFEIGENVALUESCFROMMAXTOMINDIMENSIONA(MNH,MNH),ER(MNH,4),S(MNH,MNH)TR=DO200I=1,MNHTR=TR+A(I,I)200ER(I,1)=A(I,I)MNH1=MNH-1DO210K1=MNH1,1,-1DO210K2=K1,MNH1IF(ER(K2,1).(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(*,250)TR250FORMAT(/5X,'TOTALSQUAREERROR=',RETURNENDSUBROUTINETCOEFF(KVT,KV,N,M,MNH,S,F,V,ER)CTHISSUBROUTINEPROVIDESSTANDARDEIGENVECTORSINS;CINF)ANDITSTIMECOEFFICENTSSERIESSAVEDINF;INS)DIMENSIONS(MNH,MNH),F(N,M),V(MNH),ER(MNH,4) IFTHENDO390J=1,MDO370I=1,NV(I)=F(I,J)F(I,J)=0.370CONTINUEDO380IS=1,KVTDO380I=1,N380F(IS,J)=F(IS,J)+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,MS(J,JS)=S(J,JS)*SQRT(ER(JS,1))420CONTINUEDO430I=1,NF(I,JS)=F(I,JS)/SQRT(ER(JS,1))430CONTINUEENDIFRETURNENDSUBROUTINEOUTER(KV,ER,MNH)CTHISSUBROUTINEPRINTSARRAYERCER(KV,1)FORSEQUENCEOFEIGENVALUEFROMBIGTOSMALLCER(KV,2)FOREIGENVALUEFROMBIGTOSMALLCER(KV,3)FORSMALLLO=(LAMDA/TOTALVARIANCE)CER(KV,4)FORBIGLO=SUMOFSMALLLO)DIMENSIONER(MNH,4)WRITE(16,510)510FORMAT(/10X,'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,WRITE(16,540)540FORMATE.')WRITE(16,620)(IS,IS=1,KVT)620FORMAT(3X,5i12)DO600J=1,MIFTHENWRITE(16,630)J,(f(is,j),is=1,kvt)630FORMAT(1X,I3,DO13IS=1,KVTECOF(J,IS)=F(IS,J)13CONTINUEELSEWRITE(16,640)J,(S(J,IS),IS=1,KVT)640FORMAT(1X,I3,DO14IS=1,KVTECOF(J,IS)=S(J,IS)14CONTINUEENDIF600CONTINUECWRITE(30)((S(J,IS),j=1,m),IS=1,KVT)RETURNEND实习三：大气遥相关programex3parameternt=60realt(160,nt),eu(63),h(144,73,63),hsum(144,73),have(144,73),euave,rh(144,73),r(144,73),r2(160)reala(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)realrh2(160),rh2_2(160),rup2(160)reallat(160),lon(160),timcharacter*8id(160)integernlev,nflag!读数据open(3,file='D:\3\')open(2,file='d:\3\',form='binary')open(4,file='D:\3\',form='binary')!h和eu的相关open(5,file='d:\3\',form='binary')!eu指数open(6,file='D:\3\',status='replace',form='binary')!t和eu的相关open(7,file='D:\3\') read(3,*)((t(i,j),i=1,160),j=1,nt)!读160站温度 doi=1,160 doj=1,nt t(i,j)=t(i,j)/ enddo enddo doi=1,160 read(7,*)lat(i),lon(i)!读经纬度 enddodoit=1,63 doimo=1,12 doj=1,73 doi=1,144 read(2)a(i,j,imo,it)!读高度场 enddo;enddo;enddo;enddo!计算EU指数doit=1,63eu(it)=*a(9,59,1,it)+*a(31,59,1,it)*a(59,53,1,it)enddo!!!!计算EU指数与高度场的相关系数!h--高度场hsum--高度场和have--高度场平均值eu--EU指数,eusum,euave类似!1,提取1月份高度场doit=1,63doj=1,73doi=1,144h(i,j,it)=a(i,j,1,it)enddoenddoenddo!2,计算高度场和EU指数的平均值doj=1,73doi=1,144doit=1,63hsum(i,j)=hsum(i,j)+h(i,j,it)enddohave(i,j)=hsum(i,j)/63enddoenddodoit=1,63eusum=eusum+eu(it)enddoeuave=eusum/63!3,计算相关系数各部：分子、分母、分母(对照相关系数公式)!rup--分子rh--分母hreu--分母eur--相关系数doj=1,73doi=1,144reu_2=0doit=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*,reudoj=1,73doi=1,144r(i,j)=rup(i,j)/(rh(i,j)*reu)enddoenddo!!!!计算EU指数和气温的相关系数!1,计算温度场的平均值doi=1,160doit=1,nttsum(i)=tsum(i)+t(i,it)enddotave(i)=tsum(i)/ntenddo!print*,(tave(i),i=1,160)!2,计算相关系数各部：分子、分母、分母(对照相关系数公式)!rup2--分子rh2--分母hreu--分母eur2--相关系数eusum=0doit=1,nteusum=eusum+eu(it+3)!之所以加3，是因为在分析资料和观测资料起始年份差3年enddoeuave=eusum/ntreu_2=0doit=1,ntreu_2=reu_2+(eu(it+3)-euave)**2enddoreu=sqrt(reu_2)doi=1,160doit=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))enddodoi=1,160r2(i)=rup2(i)/(rh2(i)*reu)enddo!计算完毕，写数据!写站点数据 doj=1,160 id(j)=char(j) tim= nlev=1 nflag=1 write(6)id(j),lat(j),lon(j),tim,nlev,nflag,r2(j) enddo tim= nlev=0 nflag=1 write(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）programEX4 reala6(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) realnum,rapave(160,3),h(144,73,12,65),hw(144,73,60) realhwave realhwa60(144,73,60) realhwa(144,73,3),hwat(144,73,60,3) realrapt(160,60,3),rap1(160,22),rap2(160,19),rap3(160,19) realt(144,73,3) character*8id(160) open(3,file='D:\4\')!!!请修改路径open(4,file='D:\4\') open(5,file='D:\4\') open(6,file='D:\4\') open(7,file='D:\4\ddi') open(8,file='D:\4\',form='binary') open(9,file='D:\4\',form='binary') open(10,file='D:\4\',form='binary') open(11,file='D:\4\',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) doi=1,160 read(6,*)lat(i),lon(i) enddo doit=1,60 read(7,*)(ddi(it,j),j=1,3) enddo !print*,((ddi(it,j),j=1,3),it=1,60) ccccccccccccccc编程求合成 !读高度场 doit=1,65 doimo=1,12 doj=1,73 doi=1,144 read(9)h(i,j,imo,it)!读高度场 enddo;enddo;enddo;enddo !30年降水平均值doi=1,160rsum(i)= doit=21,50 rsum(i)=rsum(i)+a6(i,it)+a7(i,it)+a8(i,it) enddo rave(i)=rsum(i)/ enddoccccccccccc请补充 !求夏季降水(3个月相加) doi=1,160 doit=1,60 a(i,it)=a6(i,it)+a7(i,it)+a8(i,it) enddo enddo !计算百分率 doi=1,160 doit=1,60 rap(i,it)=(a(i,it)-rave(i))/rave(i) enddo enddo !挑雨型,并求各个雨型的百分率平均值 doi=1,160 doj=1,3 num= doit=1,60 num=num+ddi(it,j) rapave(i,j)=rapave(i,j)+rap(i,it)*ddi(it,j) enddo rapave(i,j)=rapave(i,j)/num enddo enddo !3类雨型各自前期的冬季高度场距平 !前期，该夏季之前的那个冬季 !求冬季平均高度场hw(heightofwinter),时间往后移3年，1948+3=1951 doit=1,60 doimo=1,12 doj=1,73 doi=1,144 hw(i,j,it)=(h(i,j,1,it+3)+h(i,j,2,it+3)+h(i,j,12,it+3-1))/ enddo;enddo;enddo;enddo !求3类雨型前期冬季距平hwa1;hwa2;hwa3heightofwinterabnormal !先求60年的距平hwa60 doj=1,73 doi=1,144 hwa=0 callcal_ave(hw(i,j,:),60,hwave) doit=1,60 hwa60(i,j,it)=hw(i,j,it)-hwave enddo enddo enddo !再挑相应雨型年份的距平 doj=1,73 doi=1,144 dok=1,3 num= doit=1,60 num=num+ddi(it,k) hwa(i,j,k)=hwa(i,j,k)+hwa60(i,j,it)*ddi(it,k) enddo hwa(i,j,k)=hwa(i,j,k)/num enddo enddo enddo !对各个雨型的距平进行t检验 !t检验之前的准备：制作样本数组hwat:HeightofWinterAbnormalforTest doj=1,73 doi=1,144 dok=1,3 doit=1,60 hwat(i,j,it,k)=hwa60(i,j,it)*ddi(it,k) enddo;enddo;enddo;enddo doj=1,73 doi=1,144 dok=1,3 callforward_push(hwat(i,j,:,k),60) enddo;enddo;enddo !开始T检验 t=0 doj=1,73 doi=1,144 dok=1,3 if(k==1)kk=22 elsekk=19 callt_test(hwat(i,j,:,k),hwa60(i,j,:),kk,60,t(i,j,k)) enddo enddo enddo !写前期冬季距平数据 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类雨型百分率合成值rapaverainabnormalpercentigeaverage doj=1,160 id(j)=char(j) tim= nlev=1 nflag=1 write(8)id(j),lat(j),lon(j),tim,nlev,nflag,(rapave(j,i),i=1,3) enddo tim= nlev=0 nflag=1 write(8)id(j-1),lat(j-1),lon(j-1),tim,nlev,nflag !写t检验结果 end!以下为子程序 subroutinet_test(sam,set,nsam,nset,t) !t检验子程序 !sam样本数组set总体数组nsam样本量nset总体量 !t检验值 realsam(nsam),set(nset),t,samave,setave,sams integernsam,nset setave=0 samave=0 sams=0 callcal_ave(sam,nsam,samave) callcal_ave(set,nset,setave) callcal_s(sam,nsam,sams) t=((samave-setave)/sams)*sqrt(real(nsam)) endsubroutine subroutinecal_ave(a,na,ave) !求平均值子程序 !a数组na数组量ave平均值 reala(na),ave,asum integerna ave=0 asum=0 doi=1,na asum=asum+a(i) enddo ave=asum/real(na) endsubroutine subroutinecal_s(a,na,s) !计算标准差子程序 !a数组na数组量s标准差 reala(na),s,ave integerna ave=0 s=0 callcal_ave(a,na,ave) doi=1,na s=s+(a(i)-ave)**2 enddo s=s/real(na) s=sqrt(s) endsubroutine subroutineforward_push(a,na) !数组前缩子程序 !效果，将数组中0值后移，其他值前移 integerna reala(na),b(na) b=0 j=1 doi=1,na if(a(i)/=0)then b(j)=a(i) j=j+1 endif enddo a=0 doi=1,na a(i)=b(i) enddo endsubroutine实习五：预测因子的选择——合成分析方法（2）programEX5 reala(180,89,12,64),sst12(180,89,60),ssta(180,89,60) realsstave,ddi(60,3),sstaa(180,89,3),num realsstat(180,89,60,3),t(180,89,3) realsas(180,89,3) integerk1,k2,kk1,kk2 open(3,file='D:\5\',form='binary')!修改 open(7,file='D:\5\ddi') open(8,file='D:\5\',form='binary')ccccccccccccccc读数据（经纬度、160站降水、雨型）doit=1,64 doimo=1,12 doj=1,89 doi=1,180 read(3,end=10)a(i,j,imo,it) enddo;enddo;enddo;enddo10continue doit=1,60 read(7,*)(ddi(it,j),j=1,3) enddoccccccccccccccc编程求合成 !计算前期12月距平合成图 !先提取60年所有前期12月距平sst12:年份往后移 doj=1,89 doi=1,180 doit=1,60 sst12(i,j,it)=a(i,j,12,it+4-1) enddo enddo enddo !算距平ssta:sstabnormal doj=1,89 doi=1,180 sstave=0 callcal_ave(sst12(i,j,:),60,sstave) doit=1,60 ssta(i,j,it)=sst12(i,j,it)-sstave enddo enddo enddo !挑年份,求合成sstaa:sstabnormalasemble doj=1,89 doi=1,180 if(ssta(i,j,1)==0)then sstaa(i,j,:)=32767 else dok=1,3 num= doit=1,60 num=num+ddi(it,k) sstaa(i,j,k)=sstaa(i,j,k)+ssta(i,j,it)*ddi(it,k) enddo sstaa(i,j,k)=sstaa(i,j,k)/num enddo endif enddo enddo !做雨型合成差值sas:SstAsembleSubtract !1和2类雨型sas(:,:,1) k=1 doj=1,89 doi=1,180 sas(i,j,k)=sstaa(i,j,1)-sstaa(i,j,2) if(sas(i,j,k)==0)sas(i,j,k)=32767 enddo enddo !2和3类雨型sas(:,:,2) k=2 doj=1,89 doi=1,180 sas(i,j,k)=sstaa(i,j,2)-sstaa(i,j,3) enddo enddo !1和3类雨型sas(:,:,3) k=3 doj=1,89 doi=1,180 sas(i,j,k)=sstaa(i,j,1)-sstaa(i,j,3) enddo enddo !做t检验 !准备t检验数组sstat:sst