傲娇师兄开贴帮师妹调程序（暑假归来，继续）

离岛

有幸看到师兄的帖子是等待
EOF作业，程序出来的东西跟别人的不一样
再上传一个别人出的图

离岛

等一下，为神马下载下来要金钱啊！

离岛

我把程序贴一下
c-----*----------------------------------------------------6---------7--
C     EMPIRICAL ORTHOGONAL FUNCTIONS (EOF's)
c     This subroutine applies the EOF approach to analysis time series
c       of meteorological field f(m,n).
c     input: m,n,mnl,f(m,n),ks
c       m: number of grid-points
c       n: lenth of time series
c       mnl=min(m,n)
c       f(m,n): the raw spatial-temporal seires
c       ks: contral parameter
c           ks=-1: self; ks=0: depature; ks=1: normalized depature   !原值，距平值，标准化值
c     output: egvt,ecof,er
c       egvt(m,mnl): array of eigenvactors                           !特征向量
c       ecof(mnl,n): array of time coefficients for the respective eigenvectors
c       er(mnl,1): lamda (eigenvalues), its sequence is from big to small.
c       er(mnl,2): accumulated eigenvalues from big to small
c       er(mnl,3): explained variances (lamda/total explain) from big to small
c       er(mnl,4): accumulated explaned variances from big to small
c     work variables:
c     Last updated by Dr. Jianping Li on October 20, 2005.
c-----
        program eof_1
          implicit none
          integer i,j,k,kk,kkk,ii,jj,mnl,t
        integer,parameter::i_grid=26,j_grid=101,n=159,m_low=17,time=n
        integer,parameter::mm=i_grid*j_grid,m=2626,job=-1,num=6
c     job=-1: self; job=0: depature; job=1: normalized depature
        real,parameter::undef=32767.0
        real temp(mm,time),ersst(i_grid,j_grid,time)
        real ersst_vector(i_grid,j_grid,num)
        real nor_cof(num,time),temp_1,temp_2
        real vector_all(mm,num)
        real ff(m,time),f(m,n),er(time,4),egvt(m,time),ecof(time,n)
        real x(time),y(time),pc_low(num,time)

          open(4,file="f:\sy6\ersstjan.txt")
        open(1,file='f:\sy6\ersstjan.grd',form='binary')
                read(1) (((ersst(i,j,k),i=1,26),j=1,101),k=1,159)
          
          !print*, (((data1(i,j,k),i=1,26),j=1,101),k=1,1902)          
       
        !do k=1,time
        !t=0
        !do j=1,101
      !do i=1,26
       ! t=t+1
       !temp(t,k)=ersst(i,j,k)
        !print*, x(t,k)
      !end do
      !end do
        !end do
       


      do j=1,j_grid  
         do i=1,i_grid  
          temp(i_grid*(j-1)+i,:)=ersst(i,j,:)
         enddo
        enddo
      


        !open(2,file="f:\sy6\ersster.txt")
        !do i=1,mnl
        !write(2,"(4f30.8)") (er(i,j),j=1,4)
        !enddo
        !close(2)

      !处理缺测值
      do k=1,time
        kk=0
        do i=1,mm
        if(temp(i,k)/=undef) then
        kk=kk+1
        ff(kk,k)=temp(i,k)
        end if
        enddo
        enddo
        kkk=kk
        if(kkk>=time) then
        mnl=time
        else
        mnl=kkk
        endif
        !print*, ff
        write(4,'(159f13.3)') ((ff(kk,k),kk=1,m),k=1,n)
        call eof(kkk,time,mnl,ff,job,er,egvt,ecof)

        print*, er(1,4)
        k=0
        do j=1,j_grid
        do i=1,i_grid
        if(ersst(i,j,1)/=undef) then
        vector_all(i_grid*(j-1)+i,1:num)=egvt(i_grid*(j-1)+i-k,1:num)
        else
        k=k+1
        vector_all(i_grid*(j-1)+i,1:num)=undef
        end if
        enddo
        enddo
        do j=1,j_grid
        do i=1,i_grid
        ersst_vector(i,j,1:num)=vector_all(i_grid*(j-1)+i,1:num)
        enddo
        enddo

        do k=1,time
        open(3,file="f:\sy6\ersstecof.txt")
        write(3,rec=k) ecof(1:num,k)
        enddo
        close(3)

        do k=1,num
        open(2,file="f:\sy6\ersstegvt.grd",form='binary')
        write(2,rec=k) ersst_vector(:,:,k)
        enddo
        close(2)
        end
c-----*----------------------------------------------------6---------7--       
        subroutine eof(m,n,mnl,f,ks,er,egvt,ecof)
      dimension f(m,n),er(mnl,4),egvt(m,mnl),ecof(mnl,n)
      dimension cov(mnl,mnl),s(mnl,mnl),d(mnl),v(mnl) !work array
c---- Preprocessing data
      print *,"transf"
      call transf(m,n,f,ks)
c---- Crossed product matrix of the data f(m,n)
      call crossproduct(m,n,mnl,f,cov)
c---- Eigenvalues and eigenvectors by the Jacobi method
      call jacobi(mnl,cov,s,d,0.00001)
c---- Specified eigenvalues
      call arrang(mnl,d,s,er)
c---- Normalized eignvectors and their time coefficients
      call tcoeff(m,n,mnl,f,s,er,egvt,ecof)
      return
      end
c-----*----------------------------------------------------6---------7--
c     Preprocessing data to provide a field by ks.
c     input: m,n,f
c       m: number of grid-points
c       n: lenth of time series
c       f(m,n): the raw spatial-temporal seires
c       ks: contral parameter
c           ks=-1: self; ks=0: depature; ks=1: normalized depature
c     output: f
c       f(m,n): output field based on the control parameter ks.
c     work variables: fw(n)
      subroutine transf(m,n,f,ks)
      dimension f(m,n)
      dimension fw(n),wn(m)           !work array
        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 '
          write(*,*)' The original field has invalid data.'
          write(*,*)' There are totally ',i0,
     *    '  gridpionts with invalid data.'
               write(*,*)' The array WN stores the positions of those invalid'
        write(*,*)'   grid-points. You must pick off those invalid data'  
        write(*,*)'   from the orignal field and then reinput a new'
          write(*,*)'   field to calculate its EOFs.'   
          write(*,*)'****  FAULT  ****'
          stop
        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)
          do j=1,n
            f(i,j)=(f(i,j)-af)/sf
          enddo
        enddo
      endif
      return
      end
c-----*----------------------------------------------------6---------7--
c     Crossed product martix of the data. It is n times of
c       covariance matrix of the data if ks=0 (i.e. for anomaly).
c     input: m,n,mnl,f
c       m: number of grid-points
c       n: lenth of time series
c       mnl=min(m,n)
c       f(m,n): the raw spatial-temporal seires
c     output: cov(mnl,mnl)  
c       cov(m,n)=f*f' or f'f dependes on m and n.
c         It is a mnl*mnl real symmetric matrix.
      subroutine crossproduct(m,n,mnl,f,cov)
      dimension f(m,n),cov(mnl,mnl)
      if(n-m) 10,50,50
  10  do 20 i=1,mnl
      do 20 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)
  20  continue
      return
  50  do 60 i=1,mnl
      do 60 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)
  60  continue
      return
      end
c-----*----------------------------------------------------6---------7--
c     Computing eigenvalues and eigenvectors of a real symmetric matrix
c       a(m,m) by the Jacobi method.
c     input: m,a,s,d,eps
c       m: order of matrix
c       a(m,m): the covariance matrix
c       eps: given precision
c     output: s,d
c       s(m,m): eigenvectors
c       d(m): eigenvalues
      subroutine jacobi(m,a,s,d,eps)
      dimension a(m,m),s(m,m),d(m)
      do 30 i=1,m
      do 30 j=1,i
        if(i-j) 20,10,20
  10    s(i,j)=1.
        go to 30
  20    s(i,j)=0.
        s(j,i)=0.
  30  continue
      g=0.
      do 40 i=2,m
        i1=i-1
        do 40 j=1,i1
  40      g=g+2.*a(i,j)*a(i,j)
      s1=sqrt(g)
      s2=eps/float(m)*s1
      s3=s1
      l=0
  50  s3=s3/float(m)
  60  do 130 iq=2,m
        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)
        st=g/sqrt(2.*(1.+sqrt(1.-g*g)))
        ct=sqrt(1.-st*st)
        do 110 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
  110     s(i,ip)=g
        do 120 i=1,m
          a(ip,i)=a(i,ip)
  120     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)
  130 continue
      if(l-1) 150,140,150
  140 l=0
      go to 60
  150 if(s3.gt.s2) goto 50
      do 160 i=1,m
        d(i)=a(i,i)
  160 continue
      return
      end
c-----*----------------------------------------------------6---------7--
c     Provides a series of eigenvalues from maximuim to minimuim.
c     input: mnl,d,s
c       d(mnl): eigenvalues
c       s(mnl,mnl): eigenvectors
c     output: er
c       er(mnl,1): lamda (eigenvalues), its equence is from big to small.
c       er(mnl,2): accumulated eigenvalues from big to small
c       er(mnl,3): explained variances (lamda/total explain) from big to small
c       er(mnl,4): accumulated explaned variances from big to small
      subroutine arrang(mnl,d,s,er)
      dimension d(mnl),s(mnl,mnl),er(mnl,4)
      tr=0.0
      do 10 i=1,mnl
        tr=tr+d(i)
        er(i,1)=d(i)
  10  continue
      mnl1=mnl-1
      do 20 k1=mnl1,1,-1
      do 20 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 15 i=1,mnl
            c=s(i,k2+1)
            s(i,k2+1)=s(i,k2)
            s(i,k2)=c
  15      continue
        endif
  20  continue
      er(1,2)=er(1,1)
      do 30 i=2,mnl
        er(i,2)=er(i-1,2)+er(i,1)
  30  continue
      do 40 i=1,mnl
        er(i,3)=er(i,1)/tr
        er(i,4)=er(i,2)/tr
  40  continue
      return
      end
c-----*----------------------------------------------------6---------7--
c     Provides standard eigenvectors and their time coefficients
c     input: m,n,mnl,f,s,er
c       m: number of grid-points
c       n: lenth of time series
c       mnl=min(m,n)
c       f(m,n): the raw spatial-temporal seires
c       s(mnl,mnl): eigenvectors
c       er(mnl,1): lamda (eigenvalues), its equence is from big to small.
c       er(mnl,2): accumulated eigenvalues from big to small
c       er(mnl,3): explained variances (lamda/total explain) from big to small
c       er(mnl,4): accumulated explaned variances from big to small
c     output: egvt,ecof
c       egvt(m,mnl): normalized eigenvectors
c       ecof(mnl,n): time coefficients of egvt
      subroutine tcoeff(m,n,mnl,f,s,er,egvt,ecof)
      dimension f(m,n),s(mnl,mnl),er(mnl,4),egvt(m,mnl),ecof(mnl,n)
      dimension v(mnl)  !work array
      do j=1,mnl
        do i=1,m
          egvt(i,j)=0.
        enddo
        do i=1,n
          ecof(j,i)=0.
        enddo
      enddo
c-----Normalizing the input eignvectors s     
      do 10 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
  10  continue
c-----
      if(m.le.n) then
        do js=1,mnl
        do i=1,m
          egvt(i,js)=s(i,js)
        enddo
        enddo
        do 30 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
  30    continue
      else
        do 40 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
  40    continue
        do 50 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
  50    continue
      endif
      return
      end
c-----*----------------------------------------------------6---------7--
c     Computing the mean ax, standard deviation sx         
c       and variance vx of a series x(i) (i=1,...,n).
c     input: n and x(n)
c       n: number of raw series
c       x(n): raw series
c     output: ax, sx and vx
c       ax: the mean value of x(n) 平均值
c       sx: the standard deviation of x(n) 方差
c       vx: the variance of x(n) 标准差
c     By Dr. LI Jianping, May 6, 1998.
      subroutine meanvar(n,x,ax,sx,vx)
      dimension x(n)
      ax=0.
      vx=0.
      sx=0.
      do 10 i=1,n
        ax=ax+x(i)
  10  continue
      ax=ax/float(n)
      do 20 i=1,n
        vx=vx+(x(i)-ax)**2
  20  continue
      !print*,vx
      vx=vx/float(n)
      sx=sqrt(vx)
      return
      end

离岛

原始场做出来的是第一模态占百分之99

心语

楼主好人  可以帮我调试一下求标准化距平的程序么   查了好久还是有错

lqouc

心语 发表于 2014-7-12 14:49
楼主好人  可以帮我调试一下求标准化距平的程序么   查了好久还是有错

直接贴程序和报错截图上来。

心语

program stjp
integer ix,iy,it
parameter (x=90,y=62,t=53)
real sm(x,y,t),s(x,y),ave(x,y),ano(x,y,t),sd(x,y),anosd(x,y,t),su(x,y),sum1(x,y),ave1(x,y)
open (10,file='d:\fortran\stnew.grd',form='binary')
!!!   读数据
do iy=1,y
do ix=1,x
  do it=1,t
  read(10) sm(ix,iy,it)
  enddo
enddo
enddo
su=0
s=0
!!!   气候场
do iy=1,y
do ix=1,x
  do it=1,t
  su(ix,iy)=su(ix,iy)+sm(ix,iy,it)
  enddo
   ave(ix,iy)=(su(ix,iy))/53.0
enddo
enddo

!!!   距平场
do it=1,t
do iy=1,y
  do ix=1,x
    ano(ix,iy,it)=sm(ix,iy,it)-ave(ix,iy)

  enddo
enddo
enddo

!!!   标准差
do iy=1,y
  do ix=1,x
   do it=1,t
    s(ix,iy)=s(ix,iy)+ano(ix,iy,it)**2
    enddo
sd(ix,iy)=sqrt(s(ix,iy)/t)
enddo
enddo  

!!!   标准化距平

do iy=1,y
  do ix=1,x
   do it=1,t
   anosd(ix,iy,it)=ano(ix,iy,it)/sd(ix,iy)
sum1(ix,iy)=sum1(ix,iy)+anosd(ix,iy,it)  
  enddo
ave1(ix,iy)=sum1(ix,iy)/t  !!!!!!验证结果的正确性，因为标准化距平的平均值为0
enddo
enddo
print *,ave1

open (14,file='d:\fortran\biaozhunhua.grd',form='binary')
open (15,file='d:\fortran\biaozhunhua.txt')
do it=1,t
write(14)  ((anosd(ix,iy,it),ix=1,x),iy=1,y)
write(15,100)  ((anosd(ix,iy,it),ix=1,x),iy=1,y)
enddo

100 format(61f7.3)
end program

用的的数据stnew.grd为格点资料，x方向有90个格点，y方向有29个格点，t=53为时间长度，
PS:标准化距平是指:距平/标准差
程序没有语法上的错误，但验证的标准化距平的平均值ave1不为0，不知道哪里错了，求指点~

虫儿飞

心语 发表于 2014-7-12 15:39
program stjp
integer ix,iy,it
parameter (x=90,y=62,t=53)

对 累加 的量，在进行累加之前要赋初值，比如 0.0

虫儿飞

之前是我没仔细看，应该是这个错误，求标准偏差，开根之前需要除以样本量
!!!   标准差
do iy=1,y
  do ix=1,x
   do it=1,t
    s(ix,iy)=s(ix,iy)+ano(ix,iy,it)**2
    enddo
sd(ix,iy)=sqrt(s(ix,iy)/t)
enddo
enddo  

心语

虫儿飞 发表于 2014-7-12 18:48
之前是我没仔细看，应该是这个错误，求标准偏差，开根之前需要除以样本量
!!!   标准差
do iy=1,y

话说我开根号前除了样本量，标准差这块程序师兄你改了？