气象统计分析与预报方法 课程实验报告实验名称 实验二 经验正交函数分解系 别 大气科学 姓 名 学 号 班 级应气101实验地点 机房 实验日期 11月13日评 分指导老师肖国杰同组其他成员一、实验内容（含实验原理介绍）：实验所提供的资料为NCEP/NCAR 59年（1948年-2006年）逐年1～12月的850hPa 高度场资料，资料范围为（90N-90S ，0E -360E ），网格距为2.5*2.5，纬向格点数为144，经向格点数为73。

资料为NC 格式，资料从南到北、自西向东排列，每月为一个记录，按年逐月排放，注意读取方式以及记录长度。

对（0N -90N ，60E -120W ）850hPa 高度场进行经验正交展开（EOF.FOR ），输出分析主要参数指标；绘制环流型图和相应的时间系数序列图，并加以分析。

本实验运用EOF 方法：EOF （经验正交函数分解）是针对气象要素场进行的，其基本原理是把包含p 个空间点（变量）的场随时间变化进行分解。

设抽取样本容量为n 的资料.则场中任一空间点i 和任一时间点j 的距平观测值ij x 可看成由p 个空间函数ik v 和时间函数kj y (k=1,2,…,p)的线性组合,表示成11221pij ikkj i j i j ip pj k x vy v y v y v y ===+++∑EOF 功能是从一个气象场多次观测资料中识别出主要空间型及其时间演变规律。

EOF 展开就是将气象变量场分解为空间函数（V ）和时间函数（T ）两部分的乘积之和： X=VT 。

应用步骤：资料预处理(距平或标准化处理)计算协方差矩阵、用Jacobi 方法或迭代法计算协方差矩阵的特征值与特征向量、将特征值从大到小排列、计算特征向量的时间系数、计算每个特征向量的方差贡献、结果输出二、实验目的：经验正交函数分解（EOF）是统计天气分析中气象要素场最基础的研究模型，是必须理解和掌握的方法之一，是后续课程中许多气象要素场的计算结果的理解的基础理论，也是毕业设计和论文中的基本分析方法。

该方法用个数较少的几个空间分布模态来描述环流形势，而且基本涵盖环流场的信息，既能作为天气分析模型，其方法的延拓又能作为天气预报模型，在实际工作中也有极强的实用意义。

通过该实验，深刻理解气象要素场的统计模型的意义，掌握气象要素场分析的基本方法，为实际预报业务和科研工作打下一定的基础。

三、涉及实验的相关情况介绍（包含使用软件或实验设备等情况）：本实验要求运用Fortran6.5编译软件以及grads 计算机一台。

四、实验结果（含程序、数据记录及分析和实验总结等，可附页）：要求编写主程序，其中包括资料读入，范围截取，子程序调用。

注意：EOF的资料输入，时间场一维，空间场一维。

*********************（附程序，对关键部分标志出）**********************EOF程序C**********************************************************************C *C PROGRAM NOTES *C *C THIS PROGRAM USES EOF TO ANALYSIS TIME SERIES *C OF METEOROLOGICAL FIELD *C *C**********************************************************************C *C ******** Parameter Table ********* *C *C Mt===>LENTH OF TIME SERIES *C N ===>NUMBER OF GRID-POINTS ( or STATIONS ) *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 = Minimum(Mt,N) *C EGVT===>EIGENVECTORS, 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 *C**********************************************************************PARAMETER(N=73*37, MT=58, MNH=58)PARAMETER(KS=1, KV=10, KVT=10)REAL F(N,MT),AVF(N),DF(N),ER(MNH,4)REAL A(MNH,MNH),S(MNH,MNH),V(MNH)c**************************************************************************c INFN－输入数据文件名；OUTERA－输出特征值及方差贡献、累积方差贡献的文件名（文本）;c OUTTC1－输出时间系数文件（文本）；OUTTC2－输出时间系数文件（二进制）；c OUTTEVT－输出特征向量文件（二进制）；c**************************************************************************CHARACTER*50 INFN,OUTERA,OUTTC1,OUTTC2,OUTEVTDATA INFN/'hgt8501948-2005july.grd'/DATA OUTERA/'hgt_XT03ER3.DAT'/DATA OUTTC1/'hgt_XT03TC13.DAT'/DATA OUTTC2/'hgt_XT03TC23.DAT'/DATA OUTEVT/'hgt_XT03VT3.DAT'/C---------------- Read ORIGINAL DATA ----------------------------write(*,*)'Now is reading primative field !'OPEN (8,FILE=INFN,FORM='UNFORMATTED',ACCESS='DIRECT',RECL=N)DO IT=1,MTREAD (8,REC=IT)(F(IS,IT),IS=1,N)END DOpauseC************** START TO RUN EOF PROGRAM ******************************WRITE(*,*)write(*,*)' FIRST STEP'write(*,*)' forming the initial matrix (F) by using TRANSF !'CALL TRANSF(N,Mt,F,AVF,DF,KS)WRITE(*,*)write(*,*)' STEP 2'write(*,*)' achiving the covariance matrix by using the FORMA !' CALL FORMA(N,Mt,MNH,F,A)WRITE(*,*)write(*,*)' STEP 3 'write(*,*)' caculating the eigenvalue and eigenvectors 'WRITE(*,*)' by using Jacob method !'CALL JCB(MNH,A,S,0.001)WRITE(*,*)write(*,*)' STEP 4'write(*,*)' arrange the eigenvalue and eigenvectors' WRITE(*,*)' by using ARRANG !'CALL ARRANG(MNH,A,ER,S)WRITE(*,*)write(*,*)' STEP 5'write(*,*)' the caculation of standard eigenvectors'WRITE(*,*)' by using TCOEFF !'CALL TCOEFF(KVT,N,Mt,MNH,S,F,V,ER)write(*,*)write(*,*)' STEP 6'write(*,*)' outputing eigenvalue and accumulation using OUTER !' CALL OUTER(MNH,ER,OUTERA)WRITE(*,*)WRITE(*,*)' STEP 7'write(*,*)' outputing the time coefficient of the eigenvecters !' CALL OUTVT1(KVT,N,Mt,MNH,S,F,OUTTC1,OUTTC2)WRITE(*,*)WRITE(*,*)' STEP 8'write(*,*)' outputing the eigenvecters !'CALL OUTVT3(KVT,N,Mt,MNH,S,F,OUTEVT)ENDC *************** FINISH THE MAIN PROGRAM ***************************** C *C SUBROUTINE FUNCTION * C * 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) * C ********************************************************************* C -------- SAVING THE EIGENVALUE AND ERROR ---------------------------* SUBROUTINE OUTER(MNH,ER,OUTERA)DIMENSION ER(MNH,4)CHARACTER*50 OUTERAopen (30,file=OUTERA)WRITE(30,510)WRITE(30,520)WRITE(30,530) (IS,(ER(IS,J),J=1,4),IS=1,MNH)CLOSE(30)510 FORMAT(25X,'EIGENVALUE AND ANALYSIS ERROR')520 FORMAT(5X,'N',8X,'LAMDA',10X,'SLAMDA',11X,'PH',12X,'SPH')530 FORMAT(I6,2E15.6,2F15.5)RETURNENDC********************************************************************** C SUBROUTINE FUNCTION * C * C THIS SUBROUTINE PRINTS STANDARD EIGENVECTORS * C AND ITS TIME-COEFFICENT SERIES * C********************************************************************** C ------------- save time-coeffivcent seried of S.E. ------------SUBROUTINE OUTVT1(KVT,N,M,MNH,S,F,OUTTC1,OUTTC2)DIMENSION F(N,M),S(MNH,MNH)CHARACTER*50 OUTTC1,OUTTC2OPEN(31,file=OUTTC1)OPEN(32,file=OUTTC2,FORM='UNFORMATTED',ACCESS='DIRECT',RECL=KVT) WRITE(31,400)WRITE(31,200) (IS,IS=1,KVT)DO J=1,MIF(M.GE.N) THENWRITE(31,300) J,(F(IS,J),IS=1,KVT)WRITE(32,REC=J) (F(IS,J),IS=1,KVT)ELSEWRITE(31,300) J,(S(J,IS),IS=1,KVT)WRITE(32,REC=J) (s(J,IS),IS=1,KVT)ENDIFEND DOCLOSE(31)200 FORMAT(3X,10I15)300 FORMAT(I5,10E15.7)400 FORMAT(30X,'TIME-COEFFICENT SERIES OF S. E.')RETURNENDC --------- save standard eignvectors ------------------SUBROUTINE OUTVT3(KVT,N,M,MNH,S,F,OUTEVT)DIMENSION F(N,M),S(MNH,MNH)CHARACTER*50 OUTEVTOPEN(33,file=OUTEVT,FORM='UNFORMATTED',ACCESS='DIRECT',RECL=N)DO JS=1,KVTIF(M.GE.N) THENWRITE(33,REC=JS)(S(I,JS),I=1,N)ELSEWRITE(33,REC=JS)(F(I,JS),I=1,N)ENDIFEND DOCLOSE(33)RETURNENDC*********************************************************************** C SUBROUTINE FUNCTION * C * C THIS SUBROUTINE PROVIDES INITIAL F BY KS (optional parameter) * C ks=-1, 0, or 1 according to primative field * C*********************************************************************** SUBROUTINE TRANSF(N,M,F,AVF,DF,KS)REAL F(N,M),AVF(N),DF(N)IF (KS) 30,10,1010 DO I=1,NAVF(I)=0.0DF(I)=0.0END DODO I=1,NDO J=1,MAVF(I)=AVF(I)+F(I,J)END DOAVF(I)=AVF(I)/MDO J=1,MF(I,J)=F(I,J)-AVF(I)END DOEND DOIF (KS.EQ.1) THENDO I=1,NDO J=1,MDF(I)=DF(I)+F(I,J)*F(I,J)END DODF(I)=SQRT(DF(I)/M)DO J=1,MF(I,J)=F(I,J)/DF(I)END DOEND DOEND IF30 CONTINUERETURNENDC ----------------- FORMA -----------------------------SUBROUTINE FORMA(N,M,MNH,F,A)REAL F(N,M),A(MNH,MNH)IF (M-N) 40,50,5040 DO I=1,MNHDO J=1,IA(I,J)=0.0DO IS=1,NA(I,J)=A(I,J)+F(IS,I)*F(IS,J)END DOA(J,I)=A(I,J)END DOEND DORETURN50 DO I=1,MNHDO J=1,IA(I,J)=0.0DO JS=1,MA(I,J)=A(I,J)+F(I,JS)*F(J,JS)END DOA(J,I)=A(I,J)END DOEND DORETURNENDc*********************************************************************** C SUBROUTINE FUNCTION *C * C THIS SUBROUTINE COMPUTS EIGENVALUES AND EIGENVECTORS OF A * c*********************************************************************** SUBROUTINE JCB(N,A,S,EPS)DIMENSION A(N,N),S(N,N)DO I=1,NDO J=1,NIF (I.EQ.J) THENS(I,J)=1.0ELSES(I,J)=0.0END IFEND DOEND DOG=0.0DO I=2,NI1=I-1DO J=1,I1G=G+2.*A(I,J)*A(I,J)END DOEND DOS1=SQRT(G)S2=EPS/FLOAT(N)*S1S3=S1L=050 S3=S3/FLOAT(N)60 DO 130 IQ=2,NIQ1=IQ-1DO 130 IP=1,IQ1IF(ABS(A(IP,IQ)).LT.S3) GOTO 130L=1V1=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 I=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)*CTS(I,IP)=GEND DOA(IP,I)=A(I,IP)A(IQ,I)=A(I,IQ)END DOG=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)130 CONTINUEIF (L-1) 150,140,150140 L=0GOTO 60150 IF (S3.GT.S2) GOTO 50RETURNENDc********************************************************************** C SUBROUTINE FUNCTION * C * C THIS SUBROUTINE PROVIDES A SERIES OF EIGENVALUES * C FROM MAX TO MIN * C********************************************************************** SUBROUTINE ARRANG(MNH,A,ER,S)DIMENSION A(MNH,MNH),ER(MNH,4),S(MNH,MNH)TR=0.0DO I=1,MNHTR=TR+A(I,I)ER(I,1)=A(I,I)END DOMNH1=MNH-1DO K1=MNH1,1,-1DO K2=K1,MNH1IF(ER(K2,1).LT.ER(K2+1,1)) THENC=ER(K2+1,1)ER(K2+1,1)=ER(K2,1)ER(K2,1)=CDO I=1,MNHC=S(I,K2+1)S(I,K2+1)=S(I,K2)S(I,K2)=CEND DOEND IFEND DOEND DOER(1,2)=ER(1,1)ER(I,2)=ER(I-1,2)+ER(I,1)END DODO I=1,MNHER(I,3)=ER(I,1)/TRER(I,4)=ER(I,2)/TREND DORETURNENDC********************************************************************** C THIS SUBROUTINE PROVIDES STANDARD EIGENVECTORS * C (M.GE.N, SAVED IN S; M.LT.N, SAVED IN F) AND ITS TIME COEFFICENTS * C SERIES (M.GE.N, SAVED IN F; M.LT.N, SAVED IN S) * C********************************************************************** SUBROUTINE TCOEFF(KVT,N,M,MNH,S,F,V,ER)DIMENSION S(MNH,MNH),F(N,M),V(MNH),ER(MNH,4)DO J=1,MNHC=0.0DO I=1,MNHC=C+S(I,J)*S(I,J)END DOC=SQRT(C)DO I=1,MNHS(I,J)=S(I,J)/CEND DOEND DOIF (M.GE.N) THENDO J=1,MDO I=1,NV(I)=F(I,J)F(I,J)=0.0END DODO IS=1,KVTDO I=1,NF(IS,J)=F(IS,J)+V(I)*S(I,IS)END DOEND DOEND DOELSEDO I=1,NDO J=1,MV(J)=F(I,J)F(I,J)=0.0END DODO JS=1,KVTDO J=1,MF(I,JS)=F(I,JS)+V(J)*S(J,JS) END DO END DO END DODO JS=1,KVTDO J=1,MS(J,JS)=S(J,JS)*SQRT(ER(JS,1)) END DODO I=1,NF(I,JS)=F(I,JS)/SQRT(ER(JS,1)) END DO END DOEND IF RETURN END 结果输出:所用资料为NCEP/NCAR （0N -90N ，60E -120W ） 59年（1948年-2006年）逐年8月的850hPa 高度场资料，方差贡献及典型场分布图1 EOF 第一特征向量场图2 第一特征向量场对应的时间系数序列图3 EOF第二特征向量场图4 第二特征向量场对应的时间系数序列图5 EOF第三特征向量场图6 第三特征向量场对应的时间系数序列第一模态分析第一模态表现出大体以东经130o为界,西低东高的形势,即在东亚大陆气压较低而北太平洋地区气压较高,另从时间系数序列中可以看出这个形势总体趋势在减弱。