帮忙忙看下这个程序，菜鸟一枚看不懂

yh1223

APO定义是北半球EOF的第一模态，其指数是指第一模态的时间系数。这个程序是我别人针对春季用的，我现在想算的是1984年至2013年夏季APO指数，但是改了时间路径后生成的文件一个里面全是0，另外一个里面有840个数。。。。实在不知道该怎么改了，麻烦各位前辈帮帮忙，这个指数算不出来，后面的相关啥都分析不了。。。原来的Fortran程序（EOF）：麻烦各位看下红色区域的do循环是干什么用的？没太看明白。。。。
      parameter(in=144,jn=73,jn0=25,ny=35,nm=12)
c     input
        real  t(in,jn,nm,0:ny)
c     output
      real  tvec(in*jn0,ny)
        real  tcof(ny,ny)
        real  ter(ny,4)
c     work array
        real  tf(in*jn0,ny)
c-----------------------------------------------------------------------
      open(2,file='tzon400-200-0172-1207.grd',form='binary')     ！grd是已经处理好的温度的纬向偏差
        read(2)((((t(i,j,k1,k2),i=1,in),j=1,jn),k1=1,nm),k2=0,ny)
        close(2)

        do k=1,ny
           ll=0
         do j=37,37+jn0-1  !0-60N
           do i=1,in
              ll=ll+1
              tf(ll,k)=(t(i,j,3,k)+t(i,j,4,k)+t(i,j,5,k))/3.0   !MAM           
            weight=cos((j-37)*2.5*3.1415/180.)
              tf(ll,k)=tf(ll,k)*weight
           enddo
           enddo
        enddo

      call eof(in*jn0,ny,ny,tf,0,ter,tvec,tcof)
c-----
      open(3,file='pdo-eof-spr-60.grd',form='binary')        
        do k=1,ny
         write(3)(tvec(i,k),i=1,in*jn0)
        enddo        
        close(3)
ccccc
      open(3,file='pdo-cof-spr-60.grd',form='binary')
        do k=1,ny
         write(3)(tcof(i,k),i=1,ny)
        enddo        
        close(3)
ccccc        
      write(*,*)(ter(i,3),i=1,ny)
        write(*,*)'-------------'
        write(*,*)(ter(i,1),i=1,ny)
        xx=ter(1,1)-ter(2,1)
        yy=ter(2,1)*(float(2)/float(32))**(0.5)
        write(*,*)'xx=',xx,yy
        end
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 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     work variables:
c     Last updated by Dr. Jianping Li on October 20, 2005.
c-----
      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
      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
      vx=vx/float(n)
      sx=sqrt(vx)
      return
      end

言深深

好长的程序，而且貌似这个是fortran程序吧···

程序太长，出现这样的问题可能性比较多，建议你住不的print出来看看，到底哪一步出粗，拿程序边跑边找错还是比较容易的。

yh1223

言深深 发表于 2014-4-25 14:12
好长的程序，而且貌似这个是fortran程序吧···

程序太长，出现这样的问题可能性比较多，建议你住不的p ...

我现在知道问题出在哪里了。但是我想问下，生成的grd文件的ctl如何编写？是要根据Fortran定义的来编写吗？

言深深

yh1223 发表于 2014-4-25 14:37
我现在知道问题出在哪里了。但是我想问下，生成的grd文件的ctl如何编写？是要根据Fortran定义的来编写吗 ...

是fortran输出的时候按照特定的格式和顺序输出。具体请在论坛找程序，至于ctl描述的话已经是后来的事情了，论坛也有的