求助！EOF分析海温场fortran程序调试始终不对！

月是故乡明

这是eof程序，采用李建平版本： fortran在compile和builds时未提示任何错误，但在运行时，提示在红色句子左侧出现黄色箭头提示错误如下图

        program maineof
        implicit none

        integer temp,i,j,k !nth,r,nn,tmp
        integer,parameter :: nx=72,ny=36,nt=62 !没用到nx:xsize ny:ysize nt:时次
        integer,parameter :: m=nx*ny
        integer,parameter :: mnl=62 !min(m,nt)
        integer,parameter :: nnn=3 !模态数
        integer temp1(nx,ny)

        real data1(nx,ny,nt),data2(nx,ny,nt) !没用到ave(nx,ny)
        real,allocatable :: data0(:,:)
        real,allocatable :: egvt(:,:)
        real,parameter :: undef=-9.99e+33
        integer,parameter :: ks=0
        real er(mnl,4),ecof(mnl,nt),ecof1(nnn,nt)

        !没用到 character(len=6)::method
        character(len=100)::file1
        character(len=100)::file2
        character(len=100)::file3
        character(len=100)::file4
        character(len=100)::file5
        file1='d:\sst\sstave.dat'
        file2='d:\sst\sst_time.txt' !时间系数
        file3='d:\sst\sst_mode.dat'        !模态
        file4='d:\sst\sst_var.txt'  !方差贡献
        file5='d:\sst\sst_time.dat'

        open (55,file=file1,form='binary')
        read(55) data1
        close(55)
        
        data0=0
        temp1=1
        allocate(data0(1207,nt))   !1207为格点数
        allocate(egvt(1207,mnl))
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        do k=1,nt
        do i=1,nx
        do j=1,ny
                if(data1(i,j,k).eq.undef)temp1(i,j)=0
        enddo
        enddo
        enddo

        do k=1,nt
        temp=0
        do i=1,nx
        do j=1,ny
                if(temp1(i,j).eq.1)then
                temp=temp+1
                data0(temp,k)=data1(i,j,k)
                endif
        enddo
        enddo
        write(*,*) 'k=',k,temp
        enddo

        write(*,*) '开始进行EOF分解：'

        call eof(temp,nt,mnl,data0,ks,er,egvt,ecof)
        open(21,file=file2,form='formatted')
        do i=1,nt
        write(21,100) ecof(1,i),ecof(2,i),ecof(3,i),ecof(4,i),ecof(5,i)
        enddo
        close(21)
100        format(5f15.8)

        do i=1,nnn
        do j=1,nt
        ecof1(i,j)=ecof(i,j)
        enddo
        enddo
        open(30,file=file5,form='unformatted')
        open(31,file='d:\sst\time.txt')
        do i=1,nt
        write(30)ecof1(:,i)
        write(31,*)ecof1(:,i)
        enddo
        close(30)
        close(31)

        open(22,file=file4,form='formatted')
        write(22,*) 'lamda (eigenvalues), its sequence is from
        &big to small.'
        write(22,*) (er(i,1),i=1,mnl)
        write(22,*) 'accumulated eigenvalues from big to small'
        write(22,*) (er(i,2),i=1,mnl)
        write(22,*) 'explained variances(lamda/total explain)
        &from big to small'
        write(22,*) er(:,3)*100
        write(22,*) 'accumulated explaned variances from big to small'
        write(22,*) er(:,4)*100
        close(22)
        
        do k=1,nt
        temp=0
        do i=1,nx
        do j=1,ny
                if(temp1(i,j).eq.1)then
                temp=temp+1
                data2(i,j,k)=egvt(temp,k)
                else
                data2(i,j,k)=undef
                endif
        enddo
        enddo
        !write(*,*) 'k=',k,temp
        enddo
        
        open(23,file=file3,form='binary')
        write(23) data2(:,:,1:5)
        close(23)

        stop
        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, i.e., for the raw time series;
c           ks=0: departure, i.e., for the anomaly time series from climatological mean;
c           ks=1: normalized departure, i.e., for the normalized time series;
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     Citation: Li, J., 2005: EOF subroutine, http://web.lasg.ac.cn/staff/ljp/subroutine/EOF.f.

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.0001)
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, i.e., for the raw time series;
c           ks=0: departure, i.e., for the anomaly time series from climatological mean;
c           ks=1: normalized departure, i.e., for the normalized time series;
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
     *the original field has invalid data.'
          write(*,*)' There are totally ',i0,
     *    '  gridpionts with invalid data.'
               write(*,*)' These invalid data are those whose standard variance
     *equal zero.'
               write(*,*)' The array WN stores the positions of those invalid
     *grid-points. You must pick off those invalid data from the orignal
     *field and then reinput a new 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.
        goto 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
      goto 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)/(mnl*1.)
         er(i,1)=d(i)/(mnl*1.)
  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)*c+v(j)*s(j,js)*c !原来没有乘以C
          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)*c/sqrt(abs(er(js,1)))        !原来没有乘以C
          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

月是故乡明

Effy_yue 发表于 2014-4-2 20:02
同学，你调试好了么？最近我也在做这个，想请教一下你啊

程序没问题，最后发现是自己的初始资料出现问题，程序你可以去http://ljp.lasg.ac.cn/dct/page/65539下载，这里有许多常用的程序。

才仔

毕业论文的程序就靠你了

月是故乡明

才仔 发表于 2014-3-13 12:57

                               
登录/注册后可看大图

毕业论文的程序就靠你了

卡起了，头都大了

冬日的阳光

我也想知道，海表温度怎么做eof分析？？？

Effy_yue

同学，你调试好了么？最近我也在做这个，想请教一下你啊

s7udY

这么复杂。。。

infocert

EOF正交分解,cdo提供直接的函数可以处理。

阿勒邱

请问有与win7系统 64位兼容的Fortran版本吗？

月是故乡明

阿勒邱 发表于 2014-4-11 16:30
请问有与win7系统 64位兼容的Fortran版本吗？

有的，你去论坛搜索一下，有很多这方面的帖子