气象统计分析与预报经验正交函数分解.docx
《气象统计分析与预报经验正交函数分解.docx》由会员分享，可在线阅读，更多相关《气象统计分析与预报经验正交函数分解.docx（19页珍藏版）》请在冰豆网上搜索。
气象统计分析与预报经验正交函数分解
实验二经验正交函数分解
一、目的和要求：
经验正交函数分解（EOF）是统计天气分析中气象要素场最基础的研究模型，是必须理解和掌握的方法之一，是后续课程中许多气象要素场的计算结果的理解的基础理论，也是毕业设计和论文中的基本分析方法。
该方法用个数较少的几个空间分布模态来描述环流形势，而且基本涵盖环流场的信息，既能作为天气分析模型，其方法的延拓又能作为天气预报模型，在实际工作中也有极强的实用意义。
通过该实验，深刻理解气象要素场的统计模型的意义，掌握气象要素场分析的基本方法，为实际预报业务和科研工作打下一定的基础。
二、实验的主要内容：
对（
-
，
-
）850hPa高度场进行经验正交展开（EOF.FOR），输出分析主要参数指标；绘制环流型图和相应的时间系数序列图，并加以分析。
三、步骤：
3.1熟悉资料方法
3.1.1资料
提供的资料为NCEP/NCAR60年（1948年-2007年）逐年1～12月的850hPa高度场资料，资料范围为（
-
，
-
），网格距为2.5*2.5，纬向格点数为144，经向格点数为73。
资料为NC格式，资料从南到北、自西向东排列，每月为一个记录，按年逐月排放，注意读取方式以及记录长度。
本次实验应用NCEP/NCAR（
-
，
-
）58年（1948年-2005年）逐年7月的850hPa高度场资料，纬向格点数为73，经向格点数为37。
3.1.2方法（经验正交函数分解EOF）
EOF（经验正交函数分解）是针对气象要素场进行的，其基本原理是把包含p个空间点
（变量）的场随时间变化进行分解。
设抽取样本容量为n的资料.则场中任一空间点i和任一时间点j的距平观测值
可看成由p个空间函数
和时间函数
（k=1,2,…,p）的线性组合,表示成
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,OUTEVT
DATAINFN/'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,MT
READ（8,REC=IT）（F（IS,IT）,IS=1,N）
ENDDO
pause
C**************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）
END
C***************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*50OUTERA
open（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）
RETURN
END
C**********************************************************************
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,OUTTC2
OPEN（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,M
IF（M.GE.N）THEN
WRITE（31,300）J,（F（IS,J）,IS=1,KVT）
WRITE（32,REC=J）（F（IS,J）,IS=1,KVT）
ELSE
WRITE（31,300）J,（S（J,IS）,IS=1,KVT）
WRITE（32,REC=J）（s（J,IS）,IS=1,KVT）
ENDIF
ENDDO
CLOSE（31）
200FORMAT（3X,10I15）
300FORMAT（I5,10E15.7）
400FORMAT（30X,'TIME-COEFFICENTSERIESOFS.E.'）
RETURN
END
C---------savestandardeignvectors------------------
SUBROUTINEOUTVT3（KVT,N,M,MNH,S,F,OUTEVT）
DIMENSIONF（N,M）,S（MNH,MNH）
CHARACTER*50OUTEVT
OPEN（33,file=OUTEVT,FORM='UNFORMATTED',ACCESS='DIRECT',RECL=N）
DOJS=1,KVT
IF（M.GE.N）THEN
WRITE（33,REC=JS）（S（I,JS）,I=1,N）
ELSE
WRITE（33,REC=JS）（F（I,JS）,I=1,N）
ENDIF
ENDDO
CLOSE（33）
RETURN
END
C***********************************************************************
CSUBROUTINEFUNCTION*
C*
CTHISSUBROUTINEPROVIDESINITIALFBYKS（optionalparameter）*
Cks=-1,0,or1accordingtoprimativefield*
C***********************************************************************
SUBROUTINETRANSF（N,M,F,AVF,DF,KS）
REALF（N,M）,AVF（N）,DF（N）
IF（KS）30,10,10
10DOI=1,N
AVF（I）=0.0
DF（I）=0.0
ENDDO
DOI=1,N
DOJ=1,M
AVF（I）=AVF（I）+F（I,J）
ENDDO
AVF（I）=AVF（I）/M
DOJ=1,M
F（I,J）=F（I,J）-AVF（I）
ENDDO
ENDDO
IF（KS.EQ.1）THEN
DOI=1,N
DOJ=1,M
DF（I）=DF（I）+F（I,J）*F（I,J）
ENDDO
DF（I）=SQRT（DF（I）/M）
DOJ=1,M
F（I,J）=F（I,J）/DF（I）
ENDDO
ENDDO
ENDIF
30CONTINUE
RETURN
END
C-----------------FORMA-----------------------------
SUBROUTINEFORMA（N,M,MNH,F,A）
REALF（N,M）,A（MNH,MNH）
IF（M-N）40,50,50
40DOI=1,MNH
DOJ=1,I
A（I,J）=0.0
DOIS=1,N
A（I,J）=A（I,J）+F（IS,I）*F（IS,J）
ENDDO
A（J,I）=A（I,J）
ENDDO
ENDDO
RETURN
50DOI=1,MNH
DOJ=1,I
A（I,J）=0.0
DOJS=1,M
A（I,J）=A（I,J）+F（I,JS）*F（J,JS）
ENDDO
A（J,I）=A（I,J）
ENDDO
ENDDO
RETURN
END
c***********************************************************************
CSUBROUTINEFUNCTION*
C*
CTHISSUBROUTINECOMPUTSEIGENVALUESANDEIGENVECTORSOFA*
c***********************************************************************
SUBROUTINEJCB（N,A,S,EPS）
DIMENSIONA（N,N）,S（N,N）
DOI=1,N
DOJ=1,N
IF（I.EQ.J）THEN
S（I,J）=1.0
ELSE
S（I,J）=0.0
ENDIF
ENDDO
ENDDO
G=0.0
DOI=2,N
I1=I-1
DOJ=1,I1
G=G+2.*A（I,J）*A（I,J）
ENDDO
ENDDO
S1=SQRT（G）
S2=EPS/FLOAT（N）*S1
S3=S1
L=0
50S3=S3/FLOAT（N）
60DO130IQ=2,N
IQ1=IQ-1
DO130IP=1,IQ1
IF（ABS（A（IP,IQ））.LT.S3）GOTO130
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）
DOI=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
S（I,IP）=G
ENDDO
DOI=1,N
A（IP,I）=A（I,IP）
A（IQ,I）=A（I,IQ）
ENDDO
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）
130CONTINUE
IF（L-1）150,140,150
140L=0
GOTO60
150IF（S3.GT.S2）GOTO50
RETURN
END
c**********************************************************************
CSUBROUTINEFUNCTION*
C*
CTHISSUBROUTINEPROVIDESASERIESOFEIGENVALUES*
CFROMMAXTOMIN*
C**********************************************************************
SUBROUTINEARRANG（MNH,A,ER,S）
DIMENSIONA（MNH,MNH）,ER（MNH,4）,S（MNH,MNH）
TR=0.0
DOI=1,MNH
TR=TR+A（I,I）
ER（I,1）=A（I,I）
ENDDO
MNH1=MNH-1
DOK1=MNH1,1,-1
DOK2=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
DOI=1,MNH
C=S（I,K2+1）
S（I,K2+1）=S（I,K2）
S（I,K2）=C
ENDDO
ENDIF
ENDDO
ENDDO
ER（1,2）=ER（1,1）
DOI=2,MNH
ER（I,2）=ER（I-1,2）+ER（I,1）
ENDDO
DOI=1,MNH
ER（I,3）=ER（I,1）/TR
ER（I,4）=ER（I,2）/TR
ENDDO
RETURN
END
C**********************************************************************
CTHISSUBROUTINEPROVIDESSTANDARDEIGENVECTORS*
C（M.GE.N,SAVEDINS;M.LT.N,SAVEDINF）ANDITSTIMECOEFFICENTS*
CSERIES（M.GE.N,SAVEDINF;M.LT.N,SAVEDINS）*
C**********************************************************************
SUBROUTINETCOEFF（KVT,N,M,MNH,S,F,V,ER）
DIMENSIONS（MNH,MNH）,F（N,M）,V（MNH）,ER（MNH,4）
DOJ=1,MNH
C=0.0
DOI=1,MNH
C=C+S（I,J）*S（I,J）
ENDDO
C=SQRT（C）
DOI=1,M
升级会员 