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

1、实验二 经验正交函数分解一、目的和要求：经验正交函数分解（eof）是统计天气分析中气象要素场最基础的研究模型，是必须理解和掌握的方法之一，是后续课程中许多气象要素场的计算结果的理解的基础理论，也是毕业设计和论文中的基本分析方法。该方法用个数较少的几个空间分布模态来描述环流形势，而且基本涵盖环流场的信息，既能作为天气分析模型，其方法的延拓又能作为天气预报模型，在实际工作中也有极强的实用意义。通过该实验，深刻理解气象要素场的统计模型的意义，掌握气象要素场分析的基本方法，为实际预报业务和科研工作打下一定的基础。二、实验的主要内容：对（-，-）850hpa高度场进行经验正交展开（eof.for），输出

2、分析主要参数指标；绘制环流型图和相应的时间系数序列图，并加以分析。三、步骤：3.1 熟悉资料方法3.1.1 资料提供的资料为ncep/ncar 60年（1948年-2007年）逐年112月的850hpa高度场资料，资料范围为（-，-），网格距为2.5*2.5，纬向格点数为144，经向格点数为73。资料为nc格式，资料从南到北、自西向东排列，每月为一个记录，按年逐月排放，注意读取方式以及记录长度。本次实验应用ncep/ncar（-，-） 58年（1948年-2005年）逐年7月的850hpa高度场资料，纬向格点数为73，经向格点数为37。3.1.2 方法（经验正交函数分解eof）eof（经验正交

3、函数分解）是针对气象要素场进行的，其基本原理是把包含p个空间点（变量）的场随时间变化进行分解。设抽取样本容量为n的资料.则场中任一空间点i和任一时间点j的距平观测值可看成由p个空间函数和时间函数(k=1,2,p)的线性组合,表示成eof功能是从一个气象场多次观测资料中识别出主要空间型及其时间演变规律。eof展开就是将气象变量场分解为空间函数（v）和时间函数（t）两部分的乘积之和：x=vt。应用步骤：1) 资料预处理(距平或标准化处理)2) 计算协方差矩阵3) 用jacobi方法或迭代法计算协方差矩阵的特征值与特征向量4) 将特征值从大到小排列5) 计算特征向量的时间系数6) 计算每个特征向量的

4、方差贡献7) 结果输出3.2 编写程序要求编写主程序，其中包括资料读入，范围截取，子程序调用。注意：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

5、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)=accumulat

6、e 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输出特征值及方差贡献、累积方差贡献的文件名（文本

7、）;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

8、is reading primative field ! open (8,file=infn,form=unformatted,access=direct,recl=n) do it=1,mt read (8,rec=it)(f(is,it),is=1,n) end do pause c* 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) writ

9、e(*,*) 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(*,*) arran

10、ge the eigenvalue and eigenvectorswrite(*,*) 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 accum

11、ulation 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 * fin

12、ish 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 eigenva

13、lue and error -* subroutine outer(mnh,er,outera) dimension er(mnh,4)character*50 outera 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) 510 format(25x,eigenvalue and analysis error) 520 format(5x,n,8x,lamda,10x,slamda,11x,ph,12x,sph) 530 forma

14、t(i6,2e15.6,2f15.5)return endc*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,outtc2 open(31,file=

15、outtc1) open(32,file=outtc2,form=unformatted,access=direct,recl=kvt) write(31,400) write(31,200) (is,is=1,kvt) do j=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) endifend do close(3

16、1) 200 format(3x,10i15) 300 format(i5,10e15.7) 400 format(30x,time-coefficent series of s. e.) return endc - save standard eignvectors - subroutine outvt3(kvt,n,m,mnh,s,f,outevt) dimension f(n,m),s(mnh,mnh)character*50 outevt open(33,file=outevt,form=unformatted,access=direct,recl=n) do js=1,kvt if(

17、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) endifend do close(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

18、),avf(n),df(n) if (ks) 30,10,1010 do i=1,n avf(i)=0.0 df(i)=0.0 end do do i=1,n do j=1,m avf(i)=avf(i)+f(i,j) end do avf(i)=avf(i)/m do j=1,m f(i,j)=f(i,j)-avf(i) end do end do if (ks.eq.1) then do i=1,n do j=1,m df(i)=df(i)+f(i,j)*f(i,j) end do df(i)=sqrt(df(i)/m) do j=1,m f(i,j)=f(i,j)/df(i) end d

19、o end do end if 30 continuereturn endc - forma - subroutine forma(n,m,mnh,f,a) real f(n,m),a(mnh,mnh) if (m-n) 40,50,50 40 do i=1,mnh do j=1,i a(i,j)=0.0 do is=1,n a(i,j)=a(i,j)+f(is,i)*f(is,j) end do a(j,i)=a(i,j) end doend do return 50 do i=1,mnh do j=1,i a(i,j)=0.0 do js=1,m a(i,j)=a(i,j)+f(i,js)

20、*f(j,js) end do a(j,i)=a(i,j) end doend doreturn endc*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,n do j=1,n if (i.eq.j) then s(i,j)=1.0 else s(i,j)=0.0 end if end do end do g=0.0 do i=2,n i1=i-1 d

21、o j=1,i1 g=g+2.*a(i,j)*a(i,j) end doend do 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

22、) st=g/sqrt(2.*(1.+sqrt(1.-g*g) ct=sqrt(1.-st*st) do 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 s(i,ip)=g end do do i=1,n a(ip,i)=a(i,ip) a(iq,i)=a(i,iq) end do g=2.*v2*st*ct a(ip,ip)=v1*ct*ct+v3*st*st-g a(iq,iq)=v1*st*

23、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 goto 60 150 if (s3.gt.s2) goto 50return endc*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

24、(mnh,mnh),er(mnh,4),s(mnh,mnh) tr=0.0 do i=1,mnh tr=tr+a(i,i) er(i,1)=a(i,i)end do mnh1=mnh-1 do k1=mnh1,1,-1 do 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 i=1,mnh c=s(i,k2+1) s(i,k2+1)=s(i,k2) s(i,k2)=c end do end if end doend do er(1,2)=er(1,1) do i=2,

25、mnh er(i,2)=er(i-1,2)+er(i,1)end do do i=1,mnh er(i,3)=er(i,1)/tr er(i,4)=er(i,2)/trend doreturn endc*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

26、,mnh,s,f,v,er) dimension s(mnh,mnh),f(n,m),v(mnh),er(mnh,4) do j=1,mnh c=0.0 do i=1,mnh c=c+s(i,j)*s(i,j) end do c=sqrt(c) do i=1,mnh s(i,j)=s(i,j)/c end doend do if (m.ge.n) then do j=1,m do i=1,n v(i)=f(i,j) f(i,j)=0.0 end do do is=1,kvt do i=1,n f(is,j)=f(is,j)+v(i)*s(i,is) end do end do end do else do i=1,n do j=1,m v(j)=f(i,j) f(i,j)=0.0 end do do js=1,kvt do j=1,m f(i,js)=f(i,js)+v(j)*s(j,js) end do end do end dodo js=1,kvt d