气象统计分析与预报经验正交函数分解

实验二经验正交函数分解一、目的和要求：经验正交函数分解（EOF）是统计天气分析中气象要素场最基础的研究模型，是必须理解和掌握的方法之一，是后续课程中许多气象要素场的计算结果的理解的基础理论，也是毕业设计和论文中的基本分析方法。该方法用个数较少的几个空间分布模态来描述环流形势，而且基本涵盖环流场的信息，既能作为天气分析模型，其方法的延拓又能作为天气预报模型，在实际工作中也有极强的实用意义。通过该实验，深刻理解气象要素场的统计模型的意义，掌握气象要素场分析的基本方法，为实际预报业务和科研工作打下一定的基础。二、实验的主要内容：对（0N-90N，60E-120W）850hPa高度场进行经验正交展开（EOF.FOR），输出分析主要参数指标；绘制环流型图和相应的时间系数序列图，并加以分析。三、步骤：3.1熟悉资料方法3.1.1资料提供的资料为NCEP/NCAR60年（1948年-2007年）逐年1～12月的850hPa高度场资料，资料范围为（90N-90S，0E-360E），网格距为2.5*2.5，纬向格点数为144，经向格点数为73。资料为NC格式，资料从南到北、自西向东排列，每月为一个记录，按年逐月排放，注意读取方式以及记录长度。本次实验应用NCEP/NCAR（0N-90N，60E-120W）58年（1948年-2005年）逐年7月的850hPa高度场资料，纬向格点数为73，经向格点数为37。3.1.2方法（经验正交函数分解EOF）EOF（经验正交函数分解）是针对气象要素场进行的，其基本原理是把包含p个空间点（变量）的场随时间变化进行分解。设抽取样本容量为n的资料.则场中任一空间点i和任一时间点j的距平观测值ijx可看成由p个空间函数ikv和时间函数kjy(k=1,2,…,p)的线性组合,表示成11221pijikkjijijippjkxvyvyvyvyEOF功能是从一个气象场多次观测资料中识别出主要空间型及其时间演变规律。EOF展开就是将气象变量场分解为空间函数（V）和时间函数（T）两部分的乘积之和：X=VT。应用步骤：1)资料预处理(距平或标准化处理)2)计算协方差矩阵3)用Jacobi方法或迭代法计算协方差矩阵的特征值与特征向量4)将特征值从大到小排列5)计算特征向量的时间系数6)计算每个特征向量的方差贡献7)结果输出3.2编写程序要求编写主程序，其中包括资料读入，范围截取，子程序调用。注意：EOF的资料输入，时间场一维，空间场一维。*********************（附程序，对关键部分标志出）**********************EOF程序C**********************************************************************C*CPROGRAMNOTES*C*CTHISPROGRAMUSESEOFTOANALYSISTIMESERIES*COFMETEOROLOGICALFIELD*C*C**********************************************************************C*C********ParameterTable**********C*CMt===LENTHOFTIMESERIES*CN===NUMBEROFGRID-POINTS(orSTATIONS)*CKS=-1,SELF;KS=0,DEPATURE;KS=1,STANDERDLIZEDDEPATURE*CKV=NUMBEROFEIGENVALUESWILLBEOUTPUT*CKVT=NUMBEROFEIGENVECTORSANDTIMESERIESWILLBEOUTPUT*CMNH=Minimum(Mt,N)*CEGVT===EIGENVECTORS,ECOF===TIMECOEFFICIENTSFOREGVT*CER(KV,1)====LAMDA;LAMDA===EIGENVALUE*CER(KV,2)====ACCUMULATELAMDA*CER(KV,3)====THESUMOFCOMPONENTSVECTORSPROJECTEDONTO*CEIGENVACTOR.*CER(KV,4)====ACCUMULATEER(KV,3)*C*C**********************************************************************PARAMETER(N=73*37,MT=58,MNH=58)PARAMETER(KS=1,KV=10,KVT=10)REALF(N,MT),AVF(N),DF(N),ER(MNH,4)REALA(MNH,MNH),S(MNH,MNH),V(MNH)c**************************************************************************cINFN－输入数据文件名；OUTERA－输出特征值及方差贡献、累积方差贡献的文件名（文本）;cOUTTC1－输出时间系数文件（文本）；OUTTC2－输出时间系数文件（二进制）；cOUTTEVT－输出特征向量文件（二进制）；c**************************************************************************CHARACTER*50INFN,OUTERA,OUTTC1,OUTTC2,OUTEVTDATAINFN/'hgt8501948-2005july.grd'/DATAOUTERA/'hgt_XT03ER3.DAT'/DATAOUTTC1/'hgt_XT03TC13.DAT'/DATAOUTTC2/'hgt_XT03TC23.DAT'/DATAOUTEVT/'hgt_XT03VT3.DAT'/C----------------ReadORIGINALDATA----------------------------write(*,*)'Nowisreadingprimativefield!'OPEN(8,FILE=INFN,FORM='UNFORMATTED',ACCESS='DIRECT',RECL=N)DOIT=1,MTREAD(8,REC=IT)(F(IS,IT),IS=1,N)ENDDOpauseC**************STARTTORUNEOFPROGRAM******************************WRITE(*,*)write(*,*)'FIRSTSTEP'write(*,*)'formingtheinitialmatrix(F)byusingTRANSF!'CALLTRANSF(N,Mt,F,AVF,DF,KS)WRITE(*,*)write(*,*)'STEP2'write(*,*)'achivingthecovariancematrixbyusingtheFORMA!'CALLFORMA(N,Mt,MNH,F,A)WRITE(*,*)write(*,*)'STEP3'write(*,*)'caculatingtheeigenvalueandeigenvectors'WRITE(*,*)'byusingJacobmethod!'CALLJCB(MNH,A,S,0.001)WRITE(*,*)write(*,*)'STEP4'write(*,*)'arrangetheeigenvalueandeigenvectors'WRITE(*,*)'byusingARRANG!'CALLARRANG(MNH,A,ER,S)WRITE(*,*)write(*,*)'STEP5'write(*,*)'thecaculationofstandardeigenvectors'WRITE(*,*)'byusingTCOEFF!'CALLTCOEFF(KVT,N,Mt,MNH,S,F,V,ER)write(*,*)write(*,*)'STEP6'write(*,*)'outputingeigenvalueandaccumulationusingOUTER!'CALLOUTER(MNH,ER,OUTERA)WRITE(*,*)WRITE(*,*)'STEP7'write(*,*)'outputingthetimecoefficientoftheeigenvecters!'CALLOUTVT1(KVT,N,Mt,MNH,S,F,OUTTC1,OUTTC2)WRITE(*,*)WRITE(*,*)'STEP8'write(*,*)'outputingtheeigenvecters!'CALLOUTVT3(KVT,N,Mt,MNH,S,F,OUTEVT)ENDC***************FINISHTHEMAINPROGRAM*****************************C*CSUBROUTINEFUNCTION*C*CTHISSUBROUTINEPRINTSARRAYER*CER(KV,1)FORSEQUENCEOFEIGENVALUEFROMBIGTOSMALL*CER(KV,2)FOREIGENVALUEFROMBIGTOSMALL*CER(KV,3)FORSMALLLO=(LAMDA/TOTALVARIANCE)*CER(KV,4)FORBIGLO=SUMOFSMALLLO)*C*********************************************************************C--------SAVINGTHEEIGENVALUEANDERROR---------------------------*SUBROUTINEOUTER(MNH,ER,OUTERA)DIMENSIONER(MNH,4)CHARACTER*50OUTERAopen(30,file=OUTERA)WRITE(30,510)WRITE(30,520)WRITE(30,530)(IS,(ER(IS,J),J=1,4),IS=1,MNH)CLOSE(30)510FORMAT(25X,'EIGENVALUEANDANALYSISERROR')520FORMAT(5X,'N',8X,'LAMDA',10X,'SLAMDA',11X,'PH',12X,'SPH')530FORMAT(I6,2E15.6,2F15.5)RETURNENDC**********************************************************************CSUBROUTINEFUNCTION*C*CTHISSUBROUTINEPRINTSSTANDARDEIGENVECTORS*CANDITSTIME-COEFFICENTSERIES*C**********************************************************************C-------------savetime-coeffivcentseriedofS.E.------------SUBROUTINEOUTVT1(KVT,N,M,MNH,S,F,OUTTC1,OUTTC2)DIMENSIONF(N,M),S(MNH,MNH)CHARACTER*50OUTTC1,OUTTC2OPEN(31,file=OUTTC1)OPEN(32,file=OUTTC2,FORM='UNFORMATTED',ACCESS='DIRECT',RECL=KVT)WRITE(31,400)WRITE(31,200)(IS,IS=1,KVT)DOJ=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)ELSEWR