气象统计分析与预报经验正交函数分解Word格式文档下载.docx
《气象统计分析与预报经验正交函数分解Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《气象统计分析与预报经验正交函数分解Word格式文档下载.docx（19页珍藏版）》请在冰豆网上搜索。
2）计算协方差矩阵
3）用Jacobi方法或迭代法计算协方差矩阵的特征值与特征向量
4）将特征值从大到小排列
5）计算特征向量的时间系数
6）计算每个特征向量的方差贡献
7）结果输出
3.2编写程序
要求编写主程序，其中包括资料读入，范围截取，子程序调用。
注意：
EOF的资料输入，时间场一维，空间场一维。
*********************（附程序，对关键部分标志出）**********************
EOF程序
C**********************************************************************
C*
CPROGRAMNOTES*
CTHISPROGRAMUSESEOFTOANALYSISTIMESERIES*
COFMETEOROLOGICALFIELD*
C********ParameterTable**********
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）*
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（*,*）
FIRSTSTEP'
formingtheinitialmatrix（F）byusingTRANSF!
CALLTRANSF（N,Mt,F,AVF,DF,KS）
STEP2'
achivingthecovariancematrixbyusingtheFORMA!
CALLFORMA（N,Mt,MNH,F,A）
STEP3'
caculatingtheeigenvalueandeigenvectors'
WRITE（*,*）'
byusingJacobmethod!
CALLJCB（MNH,A,S,0.001）
STEP4'
arrangetheeigenvalueandeigenvectors'
byusingARRANG!
CALLARRANG（MNH,A,ER,S）
STEP5'
thecaculationofstandardeigenvectors'
byusingTCOEFF!
CALLTCOEFF（KVT,N,Mt,MNH,S,F,V,ER）
write（*,*）
STEP6'
outputingeigenvalueandaccumulationusingOUTER!
CALLOUTER（MNH,ER,OUTERA）
WRITE（*,*）
STEP7'
outputingthetimecoefficientoftheeigenvecters!
CALLOUTVT1（KVT,N,Mt,MNH,S,F,OUTTC1,OUTTC2）
STEP8'
outputingtheeigenvecters!
CALLOUTVT3（KVT,N,Mt,MNH,S,F,OUTEVT）
END
C***************FINISHTHEMAINPROGRAM*****************************
C*
CSUBROUTINEFUNCTION*
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
CTHISSUBROUTINEPRINTSSTANDARDEIGENVECTORS*
CANDITSTIME-COEFFICENTSERIES*
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='
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.'
C---------savestandardeignvectors------------------
SUBROUTINEOUTVT3（KVT,N,M,MNH,S,F,OUTEVT）
CHARACTER*50OUTEVT
OPEN（33,file=OUTEVT,FORM='
DOJS=1,KVT
WRITE（33,REC=JS）（S（I,JS）,I=1,N）
WRITE（33,REC=JS）（F（I,JS）,I=1,N）
CLOSE（33）
C***********************************************************************
CTHISSUBROUTINEPROVIDESINITIALFBYKS（optionalparameter）*
Cks=-1,0,or1accordingtoprimativefield*
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
DOI=1,N
AVF（I）=AVF（I）+F（I,J）
AVF（I）=AVF（I）/M
F（I,J）=F（I,J）-AVF（I）
IF（KS.EQ.1）THEN
DF（I）=DF（I）+F（I,J）*F（I,J）
DF（I）=SQRT（DF（I）/M）
F（I,J）=F（I,J）/DF（I）
ENDIF
30CONTINUE
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）
A（J,I）=A（I,J）
50DOI=1,MNH
DOJS=1,M
A（I,J）=A（I,J）+F（I,JS）*F（J,JS）
c***********************************************************************
CTHISSUBROUTINECOMPUTSEIGENVALUESANDEIGENVECTORSOFA*
SUBROUTINEJCB（N,A,S,EPS）
DIMENSIONA（N,N）,S（N,N）
DOJ=1,N
IF（I.EQ.J）THEN
S（I,J）=1.0
S（I,J）=0.0
G=0.0
DOI=2,N
I1=I-1
DOJ=1,I1
G=G+2.*A（I,J）*A（I,J）
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）
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
A（IP,I）=A（I,IP）
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）
130CONTINUE
IF（L-1）150,140,150
140L=0
GOTO60
150IF（S3.GT.S2）GOTO50
c**********************************************************************
CTHISSUBROUTINEPROVIDESASERIESOFEIGENVALUES*
CFROMMAXTOMIN*
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）
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
C=S（I,K2+1）
S（I,K2+1）=S（I,K2）
S（I,K2）=C
ER（1,2）=ER（1,1）
DOI=2,MNH
ER（I,2）=ER（I-1,2）+ER（I,1）
ER（I,3）=ER（I,1）/TR
ER（I,4）=ER（I,2）/TR
CTHISSUBROUTINEPROVIDESSTANDARDEIGENVECTORS*
C（M.GE.N,SAVEDINS;
M.LT.N,SAVEDINF）ANDITSTIMECOEFFICENTS*
CSERIES（M.GE.N,SAVEDINF;
M.LT.N,SAVEDINS）*
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
C=C+S（I,J）*S（I,J）
C=SQRT（C）
DOI=1,M
升级会员 