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

本人为南京信息工程大学大气科学系学生，在我大四上时，不幸选修课选了短期气候预测实习课程。其烦人程度超乎我的想象，一边在准备考研和找工作，一边还要花心思完成实习内容，当时根本喘不过气来。

为让学弟学妹们不重蹈覆辙，在此共享出我所用的程序，给你们一些帮助。

不过需要明确一点，这次实习对于个人编程水平的提高帮助很大，学有余力的同学，应仅把本文当做参考，理解基础上使用，而非不劳而获的资本。 以下为正文：

具体实习要求参见课本，在此不赘述

实习一：大气环流状况的表征

program EX1 real 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:\\1\\hgt500.grd',form='binary') !补充正确路径 open(4,file='d:\\1\\ave7.grd',form='binary') !补充数据输出路径 open(5,file='d:\\1\\dev.grd',form='binary') open(6,file='d:\\1\\latdev.grd',form='binary') do it=1,65 do imo=1,12 do j=1,73 do i=1,144 read(2)a(i,j,imo,it) enddo;enddo;enddo;enddo

cccccccccccccc 请完成以下的程序 !月时间平均 ave7 imo=7 do i=1,144 do j=1,73 do it=1,65 asum(i,j)=asum(i,j)+a(i,j,imo,it) enddo ave7(i,j)=asum(i,j)/65.0 enddo

enddo

!7月距平 deviation do i=1,144 do j=1,73 do it=1,65

dev(i,j,it)=a(i,j,imo,it)-ave7(i,j) enddo enddo enddo

!纬圈平均 latitude average do it=1,65 do imo=1,12 do j=1,73 do i=1,144

latsum(j,imo,it)=latsum(j,imo,it)+a(i,j,imo,it) enddo

latave(j,imo,it)=latsum(j,imo,it)/144.0 enddo enddo enddo

!纬向偏差 latitude deviation do it=1,65 do imo=1,12 do j=1,73 do i=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

实习二：大气环流分型

PROGRAM EOF

C THIS PROGRAM USES EOF FOR ANALYSING TIME SERIES C OF METEOROLOGICAL FIELD

C M:LENTH OF TIME SERIES !!!!!!!!!! m:时间序列长度 C N:NUMBER OF GRID-POINTS !!!!!!!!!! n:格点数

C KS=-1:SELF; KS=0:DEPATURE; KS=1:STANDERDLIZED DEPATURE C KV:NUMBER OF EIGENVALUES WILL BE OUTPUT

C KVT:NUMBER OF EIGENVECTORS AND TIME SERIES WILL BE OUTPUT C MNH=MIN(M,N)

C EGVT=EIGENVACTORS, ECOF=TIME COEFFICIENTS FOR EGVT. C ER(KV,1)=LAMDA,LAMDA EIGENVALUE C ER(KV,2)=ACCUMULATE LAMDA

C ER(KV,3)=THE SUM OF COMPONENTS VECTORS PROJECTED ONTO c 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:\\2\\hgt500.grd',form='binary',status='old') open(20,file='d:\\2\\egvt7.grd',form='binary') open(30,file='d:\\2\\t7.grd',form='binary') open(16,file='d:\\2\\eof7.txt')

cccccccccccccccccc读数据 do it=1,61 do k=1,12 do j=1,73 do i=1,144 read(10)hh(i,j,k,it) enddo;enddo;enddo;enddo write(*,*)'read data ok' !裁剪区域 k=1 do it=1,61 do j=1,21 do i=1,41 h(i,j,it)=hh(i+16,j+44,k,it) enddo enddo enddo write(*,*)'data narrowed' !二维空间场变一维数组，注意按照grads的XY顺序， !因为最后EIGENVACTORS文件里面直接按照该格式存的,不再经过这一步变化 do it=1,M ii=1 do j=1,21 do i=1,41 F(ii,it)=h(i,j,it) ii=ii+1 enddo enddo enddo

CCCCCCCCCCCCCCCCINPUT 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 7'

ccccccccccccc存储数据 do j=1,m do i=1,kvt write(30)ecof(j,i) enddo;enddo

do it=1,kvt do j=1,n write(20)egvt(j,it) enddo;enddo write(*,*)'ok 8'

cccccccccccc

END

ccccccccccccccccccccccccc子程序

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 EIGENVALUES C 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 ER

C ER(KV,1) FOR SEQUENCE OF EIGENVALUE FROM BIG TO SMALL C ER(KV,2) FOR EIGENVALUE FROM BIG TO SMALL

C 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 EIGENVECTORS C 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 CONTINUE

C 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 CONTINUE

C WRITE(30)((S(J,IS),j=1,m),IS=1,KVT) RETURN END

实习三：大气遥相关

program ex3 parameter nt=60

real 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),tsu

m(160)

real rh2(160),rh2_2(160),rup2(160) real lat(160),lon(160),tim character*8 id(160) integer nlev,nflag !读数据

open(3,file='D:\\3\\t1601.txt')

open(2,file='d:\\3\\hgt500.grd',form='binary')

open(4,file='D:\\3\\rheu.grd',form='binary') !h和eu的相关 open(5,file='d:\\3\\eu.grd',form='binary') !eu指数

open(6,file='D:\\3\\rteu.grd',status='replace',form='binary') !t和eu的相关 open(7,file='D:\\3\\lat_lon.txt') read(3,*)((t(i,j),i=1,160),j=1,nt)!读160站温度 do i=1,160 do j=1,nt t(i,j)=t(i,j)/10.0 enddo enddo do i=1,160 read(7,*)lat(i),lon(i)!读经纬度 enddo

do it=1,63 do imo=1,12 do j=1,73 do i=1,144 read(2)a(i,j,imo,it)!读高度场 enddo;enddo;enddo;enddo !计算EU指数 do it=1,63

eu(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,63 do j=1,73 do i=1,144

h(i,j,it)=a(i,j,1,it) enddo

enddo enddo

!2,计算高度场和EU指数的平均值 do j=1,73 do i=1,144 do it=1,63

hsum(i,j)=hsum(i,j)+h(i,j,it) enddo

have(i,j)=hsum(i,j)/63 enddo enddo

do it=1,63

eusum=eusum+eu(it) enddo

euave=eusum/63

!3,计算相关系数各部：分子、分母、分母(对照相关系数公式) !rup--分子 rh--分母h reu--分母eu r--相关系数 do j=1,73 do i=1,144 reu_2=0 do it=1,63

rup(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))**2 reu_2=reu_2+(eu(it)-euave)**2 enddo

rh(i,j)=sqrt(rh_2(i,j)) reu=sqrt(reu_2) enddo enddo

print*,reu

do j=1,73 do i=1,144

r(i,j)=rup(i,j)/(rh(i,j)*reu) enddo enddo

!!!!计算EU指数和气温的相关系数 !1,计算温度场的平均值

do i=1,160

do it=1,nt

tsum(i)=tsum(i)+t(i,it) enddo

tave(i)=tsum(i)/nt enddo

!print*, (tave(i),i=1,160)

!2,计算相关系数各部：分子、分母、分母(对照相关系数公式) !rup2--分子 rh2--分母h reu--分母eu r2--相关系数 eusum=0 do it=1,nt

eusum=eusum+eu(it+3)!之所以加3，是因为在分析资料和观测资料起始年份差3年 enddo

euave=eusum/nt

reu_2=0 do it=1,nt

reu_2=reu_2+(eu(it+3)-euave)**2 enddo

reu=sqrt(reu_2)

do i=1,160 do it=1,nt

rup2(i)=rup2(i)+(eu(it+3)-euave)*(t(i,it)-tave(i)) rh2_2(i)=rh2_2(i)+(t(i,it)-tave(i))**2 enddo

rh2(i)=sqrt(rh2_2(i)) enddo

do i=1,160

r2(i)=rup2(i)/(rh2(i)*reu) enddo

!计算完毕，写数据 !写站点数据 do j=1,160 id(j)=char(j)

tim=0.0 nlev=1 nflag=1 write(6)id(j),lat(j),lon(j),tim,nlev,nflag,r2(j) enddo tim=0.0 nlev=0 nflag=1 write(6)id(j-1),lat(j-1),lon(j-1),tim,nlev,nflag

print*, (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 EX4 real 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 hwave real 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:\\4\\r1606.txt') !!!请修改路径 open(4,file='D:\\4\\r1607.txt') open(5,file='D:\\4\\r1608.txt') open(6,file='D:\\4\\lat_lon.txt')

open(7,file='D:\\4\\ddi') open(8,file='D:\\4\\rapave.grd',form='binary') open(9,file='D:\\4\\hgt500.grd',form='binary') open(10,file='D:\\4\\hwa.grd',form='binary') open(11,file='D:\\4\\t.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,160 read(6,*)lat(i),lon(i) enddo do it=1,60 read(7,*)(ddi(it,j),j=1,3) enddo !print*,((ddi(it,j),j=1,3),it=1,60)

ccccccccccccccc 编程求合成 !读高度场 do it=1,65 do imo=1,12 do j=1,73 do i=1,144 read(9)h(i,j,imo,it)!读高度场 enddo;enddo;enddo;enddo !30年降水平均值 do i=1,160 rsum(i)=0.0 do it=21,50 rsum(i)=rsum(i)+a6(i,it)+a7(i,it)+a8(i,it) enddo rave(i)=rsum(i)/30.0 enddo

ccccccccccc 请补充 !求夏季降水(3个月相加) do i=1,160 do it=1,60 a(i,it)=a6(i,it)+a7(i,it)+a8(i,it) enddo

enddo

!计算百分率 do i=1,160 do it=1,60

rap(i,it)=(a(i,it)-rave(i))/rave(i) enddo enddo

!挑雨型,并求各个雨型的百分率平均值 do i=1,160 do j=1,3 num=0.0 do it=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(height of winter),时间往后移3年，1948+3=1951 do it=1,60 do imo=1,12 do j=1,73 do i=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))/3.0 enddo;enddo;enddo;enddo

!求3类雨型前期冬季距平hwa1;hwa2;hwa3 height of winter abnormal !先求60年的距平hwa60 do j=1,73

do i=1,144

hwa=0

call cal_ave(hw(i,j,:),60,hwave) do it=1,60

hwa60(i,j,it)=hw(i,j,it)-hwave enddo enddo enddo

!再挑相应雨型年份的距平

do j=1,73 do i=1,144 do k=1,3 num=0.0 do it=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:Height of Winter Abnormal for Test do j=1,73 do i=1,144 do k=1,3 do it=1,60

hwat(i,j,it,k)=hwa60(i,j,it)*ddi(it,k) enddo;enddo;enddo;enddo

do j=1,73 do i=1,144 do k=1,3

call forward_push(hwat(i,j,:,k),60) enddo;enddo;enddo

!开始T检验 t=0 do j=1,73 do i=1,144 do k=1,3 if(k==1) kk=22 else kk=19 call t_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类雨型百分率合成值rapave rain abnormal percentige average do j=1,160 id(j)=char(j) tim=0.0 nlev=1 nflag=1 write(8)id(j),lat(j),lon(j),tim,nlev,nflag,(rapave(j,i),i=1,3) enddo tim=0.0 nlev=0 nflag=1 write(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,sams integer nsam,nset setave=0 samave=0

sams=0

call 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 subroutine

subroutine cal_ave(a,na,ave) !求平均值子程序

!a数组 na数组量 ave平均值 real a(na),ave,asum integer na ave=0 asum=0 do i=1,na

asum=asum+a(i) enddo

ave=asum/real(na) endsubroutine

subroutine cal_s(a,na,s) !计算标准差子程序

!a数组 na数组量 s标准差 real a(na),s,ave integer na ave=0 s=0

call cal_ave(a,na,ave) do i=1,na

s=s+(a(i)-ave)**2 enddo

s=s/real(na) s=sqrt(s)

endsubroutine

subroutine forward_push(a,na) !数组前缩子程序

!效果，将数组中0值后移，其他值前移 integer na

real a(na),b(na) b=0

j=1

do i=1,na

if(a(i)/=0) then b(j)=a(i) j=j+1 endif enddo a=0

do i=1,na a(i)=b(i) enddo endsubroutine

实习五：预测因子的选择——合成分析方法（2）

program EX5 real a(180,89,12,64),sst12(180,89,60),ssta(180,89,60) real sstave,ddi(60,3),sstaa(180,89,3),num real sstat(180,89,60,3),t(180,89,3) real sas(180,89,3) integer k1,k2,kk1,kk2 open(3,file='D:\\5\\sst.grb',form='binary') !修改 open(7,file='D:\\5\\ddi') open(8,file='D:\\5\\sst12.grb',form='binary')

ccccccccccccccc 读数据（经纬度、160站降水、雨型）

do it=1,64 do imo=1,12 do j=1,89 do i=1,180 read(3,end=10)a(i,j,imo,it) enddo;enddo;enddo;enddo 10 continue do it=1,60 read(7,*)(ddi(it,j),j=1,3) enddo

ccccccccccccccc 编程求合成 !计算前期12月距平合成图 !先提取60年所有前期12月距平 sst12:年份往后移 do j=1,89 do i=1,180 do it=1,60 sst12(i,j,it)=a(i,j,12,it+4-1) enddo enddo enddo !算距平ssta:sst abnormal do j=1,89 do i=1,180 sstave=0 call cal_ave(sst12(i,j,:),60,sstave) do it=1,60 ssta(i,j,it)=sst12(i,j,it)-sstave enddo enddo enddo !挑年份,求合成sstaa: sst abnormal asemble do j=1,89 do i=1,180 if(ssta(i,j,1)==0)then sstaa(i,j,:)=32767 else do k=1,3 num=0.0 do it=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 :Sst Asemble Subtract !1和2类雨型 sas(:,:,1) k=1

do j=1,89 do i=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

do j=1,89 do i=1,180

sas(i,j,k)=sstaa(i,j,2)-sstaa(i,j,3) enddo enddo

!1和3类雨型 sas(:,:,3) k=3

do j=1,89 do i=1,180

sas(i,j,k)=sstaa(i,j,1)-sstaa(i,j,3) enddo enddo

!做t检验

!准备t检验数组sstat : sst Abnormal for Test do j=1,89 do i=1,180 do k=1,3 do it=1,60

sstat(i,j,it,k)=ssta(i,j,it)*ddi(it,k) enddo;enddo;enddo;enddo

do j=1,89 do i=1,180 do k=1,3

call forward_push(sstat(i,j,:,k),60) enddo;enddo;enddo

!开始t检验

do j=1,89 do i=1,180 do k=1,3 if(k==1) then k1=1 k2=2 kk1=22 kk2=19 elseif(k==2)then k1=2 k2=3 kk1=19 kk2=19 else k1=1 k2=3 kk1=22 kk2=19 endif call t_test2(sstat(i,j,:,k1),kk1,sstat(i,j,:,k2),kk2,t(i,j,k)) enddo enddo;enddo

! print*,t(:,:,1)

ccccccccccccccccccc写站点数据 do k=1,3 do j=1,89 do i=1,180 write(8)sstaa(i,j,k) enddo;enddo;enddo do k=1,3 do j=1,89 do i=1,180 write(8)t(i,j,k) enddo;enddo;enddo do k=1,3 do j=1,89 do i=1,180 write(8)sas(i,j,k) enddo;enddo;enddo

close(8) end

!以下为子程序

subroutine t_test2(sam1,nsam1,sam2,nsam2,t) !t检验子程序（样本之间）

!sam1,sam2两个样本数组 nsam1,nsam2样本数组容量 !t--t统计量

integer nsam1,nsam2

real sam1(nsam1),sam2(nsam2),t real ave1,ave2,s1,s2,ma,mb t=0 ave1=0 ave2=0 s1=0 s2=0 ma=1 mb=1

call cal_ave(sam1,nsam1,ave1) call cal_ave(sam2,nsam2,ave2) call cal_s(sam1,nsam1,s1) call cal_s(sam2,nsam2,s2)

!缺测值的处理：先判断，再赋值 if(s1==0.and.s2==0)then t=32767 else

ma=sqrt(((nsam1-1)*(s1**2)+(nsam2-1)*(s2**2))/(nsam1+nsam2-2)) mb=sqrt((1.0/real(nsam1))+(1.0/real(nsam2))) t=(ave1-ave2)/(ma*mb) endif end subroutine

subroutine cal_ave(a,na,ave) !求平均值子程序

!a数组 na数组量 ave平均值 real a(na),ave,asum integer na ave=0 asum=0 do i=1,na

asum=asum+a(i) enddo

ave=asum/real(na) endsubroutine

subroutine cal_s(a,na,s) !计算标准差子程序

!a数组 na数组量 s标准差 real a(na),s,ave integer na ave=0 s=0

call cal_ave(a,na,ave) do i=1,na

s=s+(a(i)-ave)**2 enddo

s=s/real(na) s=sqrt(s)

endsubroutine

subroutine forward_push(a,na) !数组前缩子程序

!效果，将数组中0值后移，其他值前移 integer na

real a(na),b(na) b=0 j=1

do i=1,na

if(a(i)/=0) then b(j)=a(i) j=j+1 endif enddo a=0

do i=1,na a(i)=b(i) enddo endsubroutine

ave=asum/real(na) endsubroutine

subroutine cal_s(a,na,s) !计算标准差子程序

!a数组 na数组量 s标准差 real a(na),s,ave integer na ave=0 s=0

call cal_ave(a,na,ave) do i=1,na

s=s+(a(i)-ave)**2 enddo

s=s/real(na) s=sqrt(s)

endsubroutine

subroutine forward_push(a,na) !数组前缩子程序

!效果，将数组中0值后移，其他值前移 integer na

real a(na),b(na) b=0 j=1

do i=1,na

if(a(i)/=0) then b(j)=a(i) j=j+1 endif enddo a=0

do i=1,na a(i)=b(i) enddo endsubroutine

本文来源：https://www.bwwdw.com/article/l44r.html