菜鸟求帮助，帮看下程序谢谢

ayzqs

test.txt (1.56 KB, 下载次数: 0)

2012-8-21 13:06 上传

点击文件名下载附件
资料


! EOF：经验正交函数分解!  “实现时空场的分离”具体可以参考黄嘉佑老师的《气象统计分析与预报方法》
! 这是来自“lijianping”老师的EOF分解程序
! 相关参数说明
! m：格点数
! n：时间长度
! mnl：模态数   
! ks：[1]-1表示对原场EOF；
!  [2]0表示多距平场EOF；
!  [3]1表示对归一化场EOF。
! "year.mean.txt" ：需要分解的时空场
! "er.txt"：存放解释方差 [1]：er(mnl,1)特征向量—— [2]：er(mnl,2)累计特征向量—— [3]：er(mnl,3)解释方差—— [4]：er(mnl,4)累计解释方差
! "ecof.txt"：存放时间系数------每一列是一个模态
! "egvt.txt"：存放各个模态的空间向量场------每一列是一个模态
parameter(m=5,n=30,mnl=5,ks=0)
real x(m,n),egvt(m,mnl),ecof(mnl,n),er(mnl,4)
open(1,file="test.txt")
do i=1,m
read(1,*) (x(i,j),j=1,n)
enddo
close(1)
call eof(m,n,mnl,x,ks,er,egvt,ecof)
open(2,file="er.txt")
do i=1,n
write(2,"(4f30.8)") (er(i,j),j=1,4)
enddo
close(2)
open(3,file="ecof.txt")
do j=1,n
write(3,"(<mnl>f20.6)") (ecof(i,j),i=1,mnl)
enddo
close(3)
open(4,file="egvt.txt")
do i=1,m
write(4,"(<mnl>f20.6)") (egvt(i,j),j=1,mnl)
enddo
close(4)
end
!---------------------------------------------------------------
subroutine eof(m,n,mnl,f,ks,er,egvt,ecof)
  implicit none
  integer*4    :: m,n,mnl,ks
  real*4     :: f(m,n),er(mnl,4),egvt(m,mnl),ecof(mnl,n)
  real*4,allocatable  :: cov(:,:),s(:,:),d(:) !---work array
!--------------------preprocessing data
  call transf(m,n,f,ks)
  allocate(cov(mnl,mnl))
!--------------------crossed product matrix of the data f(m,n)
  call crossproduct(m,n,mnl,f,cov)
  
  allocate(s(mnl,mnl))
  allocate(d(mnl))
!--------------------eigenvalues and eigenvectors by the jacobi method
  call jacobi(mnl,cov,s,d,0.00001)
  
!--------------------specified eigenvalues
  call arrang(mnl,d,s,er)
!--------------------normalized eignvectors and their time coefficients
  call tcoeff(m,n,mnl,f,s,er,egvt,ecof)
  return
end
!-----------------------------------------------------------

!--------------------------------------------------------------
!---preprocessing data to provide a field by ks.
!---input: m,n,f
!------m: number of grid-points
!------n: lenth of time series
!------f(m,n): the raw spatial-temporal seires
!------ks: contral parameter:ks=-1: self; ks=0: depature; ks=1: normalized depature
!---output: f
!------f(m,n): output field based on the control parameter ks.
!------work variables: fw(n)
!----------------------------------------------------------------
subroutine transf(m,n,f,ks)
  implicit none
  integer*4     :: m,n,ks
  real*4      :: f(m,n)
  real*4,allocatable   :: fw(:),wn(:)     !---work array
  integer*4     :: i0,i,j
  real*4      :: af,sf,vf
  
  allocate(fw(n))
  allocate(wn(m))
  i0=0
  do i=1,m
   do j=1,n
    fw(j)=f(i,j)
   enddo
   call meanvar(n,fw,af,sf,vf)
   if(sf.eq.0.)then
    i0=i0+1
    wn(i0)=i
   endif
  enddo
  if(i0.ne.0)then
   write(*,*)'*************************  fault  ***********************************'
   write(*,*)'The program cannot go on because the original field has invalid data.'
   write(*,*)'There are totally ',i0,'gridpionts with invalid data.'
   write(*,*)'The array wn(m) stores the positions of those invalid grid-points.'
   write(*,*)'You must pick off those invalid data from the orignal field---$ '  
   write(*,*)'$---and then reinput a new field to calculate its eofs.'
   write(*,*)'************************  fault  ************************************'
  endif     
  if(ks.eq.-1)return
  if(ks.eq.0)then                !anomaly of f
   do i=1,m
    do j=1,n
     fw(j)=f(i,j)
    enddo
    call meanvar(n,fw,af,sf,vf)
    do j=1,n
     f(i,j)=f(i,j)-af
    enddo
   enddo
   return
  endif
  if(ks.eq.1)then                 !normalizing f
   do i=1,m
    do j=1,n
     fw(j)=f(i,j)
    enddo
    call meanvar(n,fw,af,sf,vf)
    if(sf==0.0)then
     do j=1,n
      f(i,j)=0.0
     enddo
    else
     do j=1,n
      f(i,j)=(f(i,j)-af)/sf
     enddo
    endif
   enddo
  endif
  return
end
!--------------------------------------------------------

!----------------------------------------------------------
!---crossed product martix of the data. it is n times of covariance matrix of the data if ks=0 (i.e. for anomaly).
!---input: m,n,mnl,f
!---m: number of grid-points
!---n: lenth of time series
!---mnl=min(m,n)
!---f(m,n): the raw spatial-temporal seires
!---output: cov(mnl,mnl)  
!---cov(m,n)=f*f' or f'f dependes on m and n.
!---it is a mnl*mnl real symmetric matrix.
subroutine crossproduct(m,n,mnl,f,cov)
  implicit none
  integer*4  :: m,n,mnl
  real*4   :: f(m,n),cov(mnl,mnl)
  integer*4  :: i,j,is,js
  
   if((n-m)<0)then
    do i=1,mnl
     do j=i,mnl
      cov(i,j)=0.0
      do is=1,m
       cov(i,j)=cov(i,j)+f(is,i)*f(is,j)
      enddo
      cov(j,i)=cov(i,j)
     enddo
    enddo
   else
    do i=1,mnl
     do j=i,mnl
      cov(i,j)=0.0
      do js=1,n
       cov(i,j)=cov(i,j)+f(i,js)*f(j,js)
      enddo
      cov(j,i)=cov(i,j)
     enddo
    enddo
   endif
   
  return
end
!----------------------------------------------------------


!----------------------------------------------------------
!---computing eigenvalues and eigenvectors of a real symmetric matrix
!---a(m,m) by the jacobi method.
!---input: m,a,s,d,eps
!------m: order of matrix
!------a(m,m): the covariance matrix
!------eps: given precision
!---output: s,d
!------s(m,m): eigenvectors
!------d(m): eigenvalues
subroutine jacobi(m,a,s,d,eps)
  implicit none
  integer*4  :: m
  real*4   :: a(m,m),s(m,m),d(m),eps
  integer*4  :: i,j,i1,l,iq,iq1,ip
  real*4   :: g,s1,s2,s3,v1,v2,v3,u,st,ct
   
   s=0.0
   do i=1,m
    s(i,i)=1.0
   enddo
   g=0.
   do i=2,m
    i1=i-1
    do j=1,i1
     g=g+2.0*a(i,j)*a(i,j)
    enddo
   enddo
   s1=sqrt(g)
   s2=eps/float(m)*s1
   s3=s1
   l=0
50   s3=s3/float(m)
60   do iq=2,m
    iq1=iq-1     
    do ip=1,iq1
     if(abs(a(ip,iq)).lt.s3) exit
     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)
      
     do i=1,m
      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
     do i=1,m
      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)
    enddo
   enddo
   if((l-1)==0)then
    l=0
    go to 60
   else
    go to 150
   endif
150   if(s3.gt.s2) then
    goto 50
   end if
   do i=1,m
    d(i)=a(i,i)
   enddo
  return
end
!-----------------------------------------------------

!-----------------------------------------------------
!---provides a series of eigenvalues from maximuim to minimuim.
!---input: mnl,d,s
!------d(mnl): eigenvalues
!------s(mnl,mnl): eigenvectors
!---output: er
!------er(mnl,1): lamda (eigenvalues), its equence is from big to small.
!------er(mnl,2): accumulated eigenvalues from big to small
!------er(mnl,3): explained variances (lamda/total explain) from big to small
!------er(mnl,4): accumulated explaned variances from big to small
subroutine arrang(mnl,d,s,er)
  implicit none
  integer*4  :: mnl
  real*4   :: d(mnl),s(mnl,mnl),er(mnl,4)
  integer*4  :: i,mnl1,k1,k2
  real*4   :: c,tr
  
  tr=0.0
  do i=1,mnl
   tr=tr+d(i)
   er(i,1)=d(i)
  enddo
  mnl1=mnl-1
  do k1=mnl1,1,-1
   do k2=k1,mnl1
    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,mnl
      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)
  
  do i=2,mnl
   er(i,2)=er(i-1,2)+er(i,1)
  enddo
  do i=1,mnl
   er(i,3)=er(i,1)/tr
   er(i,4)=er(i,2)/tr
  enddo
  return
end
!--------------------------------------------

!--------------------------------------------
!---provides standard eigenvectors and their time coefficients
!---input: m,n,mnl,f,s,er
!------m: number of grid-points
!------n: lenth of time series
!------mnl=min(m,n)
!------f(m,n): the raw spatial-temporal seires
!------s(mnl,mnl): eigenvectors
!------er(mnl,1): lamda (eigenvalues), its equence is from big to small.
!------er(mnl,2): accumulated eigenvalues from big to small
!------er(mnl,3): explained variances (lamda/total explain) from big to small
!------er(mnl,4): accumulated explaned variances from big to small
!---output: egvt,ecof
!------egvt(m,mnl): normalized eigenvectors
!------ecof(mnl,n): time coefficients of egvt
subroutine tcoeff(m,n,mnl,f,s,er,egvt,ecof)
  implicit none
  integer*4     :: m,n,mnl
  real*4      :: f(m,n),s(mnl,mnl),er(mnl,4),egvt(m,mnl),ecof(mnl,n)
  real*4,allocatable   :: v(:)  !---work array
  integer*4     :: i,j,js,is
  real*4      :: c
  
  allocate(v(mnl))
  do j=1,mnl
   do i=1,m
    egvt(i,j)=0.
   enddo
   do i=1,n
    ecof(j,i)=0.
   enddo
  enddo
!---normalizing the input eignvectors s
  do j=1,mnl
   c=0.
   do i=1,mnl
    c=c+s(i,j)*s(i,j)
   enddo
   c=sqrt(c)
   do i=1,mnl
    s(i,j)=s(i,j)/c
   enddo
  enddo
!-----------------------------------------------
  if(m.le.n) then
   do js=1,mnl
    do i=1,m
     egvt(i,js)=s(i,js)
    enddo
   enddo
   do j=1,n
    do i=1,m
     v(i)=f(i,j)
    enddo
    do is=1,mnl
     do i=1,m
      ecof(is,j)=ecof(is,j)+v(i)*s(i,is)
     enddo
     enddo
   enddo
  else
   do i=1,m
    do j=1,n
     v(j)=f(i,j)
    enddo
    do js=1,mnl
     do j=1,n
      egvt(i,js)=egvt(i,js)+v(j)*s(j,js)
     enddo
    enddo
   enddo
   do js=1,mnl
    do j=1,n
     ecof(js,j)=s(j,js)*sqrt(abs(er(js,1)))
    enddo
    do i=1,m
     egvt(i,js)=egvt(i,js)/sqrt(abs(er(js,1)))
    enddo
   enddo
  endif
  
  return
end
!---------------------------------------------
!---------------------------------------------
!---computing the mean ax, standard deviation sx and variance vx of a series x(i) (i=1,...,n).
!---input: n and x(n)
!------n: number of raw series
!------x(n): raw series
!---output: ax, sx and vx
!------ax: the mean value of x(n)
!------sx: the standard deviation of x(n)
!------vx: the variance of x(n)
!------by dr. li jianping, may 6, 1998.
subroutine meanvar(n,x,ax,sx,vx)
  implicit none
  integer*4  :: n
  real*4   :: x(n)
  real*4   :: ax,vx,sx
  integer*4  :: i
  
  ax=0.
  vx=0.
  sx=0.
  do i=1,n
   ax=ax+x(i)
  enddo   
  ax=ax/float(n)
  do i=1,n
   vx=vx+(x(i)-ax)**2
  enddo
  vx=vx/float(n)
  sx=sqrt(vx)
  
  return
end

ayzqs

本人很菜，运行老是出现Array bounds exceeded

mofangbao

你不觉得贴这么多没啥必要么，还不会自己找错的话就先学着找错吧，贴这么一大段没人有闲工夫会给你看的啊，说的直接了一些，但事实的确是这样，所以楼主学会先把错误定位准确了，在来讨论吧。你先把整个程序的逻辑搞清楚，执行的流程是什么，在然后一步步的暂停调试，直到错误出现

做个霸气的木头

程序好长啊~不过这类数组越界，一般是数组定义长度与实际数据长度不对应，或者一些循环次数与数组长度不对应造成的，楼主要搞清程序的思路，一个个循环的检查。

ayzqs

谢谢，fortran都交给老师了，重新拾起感觉什么都不会了

wode_q_x

不着急，慢慢来，一行一行的查看

ayzqs

mofangbao 发表于 2012-8-21 13:57

                               
登录/注册后可看大图

你不觉得贴这么多没啥必要么，还不会自己找错的话就先学着找错吧，贴这么一大段没人有闲工夫会给你看的啊， ...

最近想用EOF分析方法分析点东西，要是会fortran就好了，掉了在自学可不容易呀

言深深

额，这个不是我之前贴的程序么，请楼主看看帖子前面有个介绍，之前有朋友反映程序27行有个错误，楼主不妨回到原帖子看看吧！我已经添加说明了

ayzqs

言深深 发表于 2012-8-21 20:23

                               
登录/注册后可看大图

额，这个不是我之前贴的程序么，请楼主看看帖子前面有个介绍，之前有朋友反映程序27行有个错误，楼主不妨回 ...

谢谢版主提醒！

言深深

ayzqs 发表于 2012-8-21 20:30

                               
登录/注册后可看大图

谢谢版主提醒！

我刚刚重新更新了一下帖子里面的程序，运行时没有问题了。挺抱歉的···希望你成功