小波程序有错求帮忙改正 急急！！！！！

wysyhd

求大神帮忙改正，谢谢了，正在做毕业论文     急急！！！！！
    program wwave
c      this file use the 93' data in nh( 105.0-120.E, 4.762-20.N)
       parameter(n=50,dt=1.0,s0=2*dt,dj=0.25,m=18)
       dimension sst(n),recon(n),con(n),wte(n,m),sst1(n)
       dimension s(m),period(m),dof(m),fft_theo(m),signif(m)
       dimension glob_ws(m),wtm(n,m),wtr(n,m),glob_signif(m)
       complex wave(n,m)
       real lag1
c
        open(1,file='E:\wavelet\data.grd',form='binary',
     *recl=4)
c       read(1,*)sst
       do 1 i=1,n
       read(1) sst(i)
1     continue
       close(1)
       print*,'sst()',sst(1),sst(n-1),sst(n)
c
c      ---------canceled by hhy
c       call corr(sst,n,lag1,r1,r2)
c       write(*,*)'lag1=',lag1
c
       pi=4.0*atan(1.0)
       ibase2=nint(alog(float(n))/alog(2.0))+1
       npad=int(2.0**ibase2)
       print*,'step 1'
       write(*,*)npad
       call wavelet(n,sst,dt,s0,dj,m,wave,s,period,con,npad)
c  significance testing
c  local significance test
       isigtest=0
       siglvl=0.05
       call wave_signif(isigtest,n,sst,dt,dj,m,s,period,
     &      lag1,siglvl,dof,fft_theo,signif,ymean,variance,cdelta,psi0)
c  global wavelet spectrum & signficance test
       isigtest=1
       siglvl=0.05
       do 10 j=1,m
         glob_ws(j)=0.
         do 20 i=1,n
           glob_ws(j)=glob_ws(j)+cabs(wave(i,j))**2
20      continue
         glob_ws(j)=glob_ws(j)/(n*variance)
         dof(j)=n-s(j)
10    continue
       call wave_signif(isigtest,n,sst,dt,dj,m,s,period,
     &      lag1,siglvl,dof,fft_theo,glob_signif,ymean,
     &      variance,cdelta,psi0)
       delta=0.
       do 31 i=1,n
       do 31 j=1,m
         delta=delta+(cabs(wave(i,j))**2)/s(j)
31    continue
       delta=delta*dj*dt/(cdelta*n)
       print*,'delta,variance,ymean'
       write(*,*)delta,variance,ymean
c  reconstruct the time series
       recon_mean=0.
       recon_vari=0.
       do 32 i=1,n
         recon(i)=0.
         do 33 j=1,m
33        recon(i)=recon(i)+(real(wave(i,j)))/sqrt(s(j))
         recon(i)=dj*sqrt(dt)*recon(i)/(cdelta*psi0)        
         recon_vari=recon_vari+(sst(i)-ymean-recon(i))**2
         recon_mean=recon_mean+recon(i)
32    continue
       recon_mean=recon_mean/n
       recon_vari=sqrt(recon_vari/n)
       print*,'recon_mean,recon_vari'
       write(*,*)recon_mean,recon_vari

       do 81 i=1,m
       do 81 j=1,n
         wtm(j,i)=cabs(wave(j,i))**2/variance
         wtr(j,i)=real(wave(j,i))
         wte(j,i)=cabs(wave(j,i))
81    continue
c
       open(21,file='E:\wavelet\nhr_re.dat',
     *access='direct',form='binary',recl=n*m*4)
       open(61,file='E:\wavelet\nhr_mse.dat',
     *access='direct',form='binary',recl=n*m*4)
       open(31,file='E:\wavelet\nhr_mod.dat',
     *access='direct',form='binary',recl=n*m*4)
c       open(61,file='e:\data\nh82md.dat')
c       open(31,file='e:\data\nh90en.dat')
c       open(4,file='f:\data\nh90recn.dat')
c       open(5,file='f:\data\nh90glb_ws.dat')


        write(21,rec=1) ((wtr(j,i),j=1,n),i=1,m)
        write(31,rec=1) ((wtm(j,i),j=1,n),i=1,m)
        write(61,rec=1) ((wte(j,i),j=1,n),i=1,m)
80    continue
c       write(4,200)(recon(i),i=1,72)
c       do 40 i=1,m
c         write(5,300)period(i),s(i)
c         write(5,301)i,glob_ws(i),glob_signif(i)
c 40    continue

200   format(10f7.2)
300   format(5x,'period is',f7.3,3x,'scale is',f7.3)
301   format(4x,i5,2(6x,f8.3))

       stop
       end
c-------------------------------------------------------------
       subroutine fast(n,x,nd,ind)
       complex x(nd),temp,theta

       j=1
       do 140 i=1,n
         if(i.ge.j)goto 110
         temp=x(j)
         x(j)=x(i)
         x(i)=temp
110     m=n/2
120     if(j.le.m)goto 130
         j=j-m
         m=m/2
         if(m.ge.2)goto 120
130     j=j+m
140   continue
       kmax=1
150   if(kmax.ge.n)return
       istep=kmax*2
       do 170 k=1,kmax
         theta=cmplx(0.0,3.141593*float(ind*(k-1))/float(kmax))
         do 160 i=k,n,istep
           j=i+kmax
           temp=x(j)*cexp(theta)
           x(j)=x(i)-temp
           x(i)=x(i)+temp
160     continue
170   continue
       kmax=istep
       goto 150

       return
       end
c-------------------------------------------------------------
       subroutine corr(sst,n,r,r1,r2)
       dimension sst(n)
       a=0.0
       s=0.0
       do 1 i=1,n
         a=a+sst(i)
1     continue
       a=a/float(n)
       do 2 i=1,n
         s=s+(sst(i)-a)**2
2     continue  
       s=sqrt(s/float(n))
       r1=0.
       do 4 i=1,n-1
         r1=r1+(sst(i)-a)*(sst(i+1)-a)
4     continue  
       r1=r1/(float(n-1)*s*s)
       r2=0.
       do 5 i=1,n-2
         r2=r2+(sst(i)-a)*(sst(i+2)-a)
5     continue  
       r2=r2/(float(n-2)*s*s)
       r=(r1+sqrt(r2))/2.0

       return
       end
c-----------------------------------------------
       subroutine wavelet(n,y,dt,s0,dj,jtot,
     &                    wave,s,period,coi,npad)
c       parameter(npad=128)
       real y(n),dt,s0,dj,s(jtot),period(jtot),coi(n)
       complex wave(n,jtot)
       real ymean,freq1,pi,period1,coi1,kwave(npad)
       complex yfft(npad),daughter(npad)

       pi=4.0*atan(1.0)
       if(npad.lt.n)then
         print*,'wavelet:"npad" must be greater than or equal to "n"'
         return
       endif
c  find the time-series mean & remove it      
       ymean=0.0
       do 10 i=1,n
         ymean=ymean+y(i)
10    continue
       ymean=ymean/n
       do 15 i=1,n
         yfft(i)=y(i)-ymean
15    continue

c  if desired,pad with extra zeroes
       do 20 i=n+1,npad
         yfft(i)=0.0
20    continue

c  find the FFT of the time series [Eqn(3)]                    
       call fast(npad,yfft,npad,-1)
       do 30 k=1,npad
         yfft(k)=yfft(k)/npad
30    continue

c  construct the wavenumber array [Eqn(5)]
       freq1=2.0*pi/(float(npad)*dt)
       kwave(1)=0.0
       do 40 i=2,npad/2+1
         kwave(i)=(float(i)-1.0)*freq1
40    continue
       do 50 i=npad/2+2,npad
         kwave(i)=-kwave(npad-i+2)
50    continue

c  main wavelet loop
       do 100 j=1,jtot
         s(j)=s0*(2.0**(float(j-1)*dj))
         call wave_function(npad,dt,s(j),kwave,period1,coi1,daughter)
         period(j)=period1
c   multiply the daughter by the time-series FFT
         do 60 k=1,npad
           daughter(k)=daughter(k)*yfft(k)
60      continue
c   inverse FFT [Eqn(4)]
         call fast(npad,daughter,npad,1)
c   store the wavelet tranform, discard,zero-padding at end
         do 70 i=1,n
           wave(i,j)=daughter(i)
70      continue
100   continue
       print*,'step2  s(1),s(17)='
       write(*,*)s(1),s(17)

c  construct the cone of influence
       do 110 i=1,(n+1)/2
         coi(i)=coi1*dt*(float(i)-1.0)
         coi(n-i+1)=coi(i)
110   continue

       return
       end
c-----------------------------------------
       subroutine wave_function(nk,dt,scale1,kwave,
     &                         period1,coi1,daughter)
       real dt,kwave(nk),scale1,period1,coi1
       complex daughter(nk)
       real expnt,pi,norm,fourier_factor

       pi=4.0*atan(1.0)
       norm=sqrt(2.*pi*scale1/dt)*(pi**(-0.25))
       do 10 k=1,nk/2+1
         expnt=-0.5*(scale1*kwave(k)-6.0)**2
         daughter(k)=cmplx(norm*exp(expnt))
10    continue
       do 20 k=nk/2+2,nk
         daughter(k)=cmplx(0.0)
20    continue
       fourier_factor=(4.*pi)/(6.0+sqrt(2.+6.**2))
       period1=scale1*fourier_factor
       coi1=fourier_factor/sqrt(2.0)

       return
       end                  
c-----------------------------------------------------
       subroutine wave_signif(isigtest,n,y,dt,dj,jtot,s,period,lag1,
     &        siglvl,dof,fft_theo,signif,ymean,variance,cdelta,psi0)
       real y(n),dt,dj,s(jtot),period(jtot)
       real lag1,siglvl,dof(jtot),fft_theo(jtot),signif(jtot)
       real ymean,variance,cdelta,psi0
       real pi,freq1,dofmin,dj0,savg,smid
       real fft_theor1

       pi=4.0*atan(1.0)
       siglvl=0.05
       if(lag1.le.0.)lag1=0.0
       dofmin=2.
       cdelta=0.776
       gammafac=2.32
       dj0=0.60
       psi0=pi**(-0.25)
c  find the time-series variance
       ymean=0.0
       do 10 i=1,n
         ymean=ymean+y(i)
10    continue
       ymean=ymean/n
       variance=0.0
       do 15 i=1,n
         variance=variance+(y(i)-ymean)**2
15    continue
       variance=variance/(float(n))

c  construct theoretical red(white)-noise power spectrum [Eqn(16)]
       do 20 j=1,jtot
         freq1=dt/period(j)
         fft_theo(j)=(1.-lag1**2)/(1.-2.*lag1*cos(freq1*2.*pi)
     &                               +lag1**2)
20    continue
       if(isigtest.eq.0)then
c----no smoothing, dof=dofmin----
c  see Eqn(18)
        chisqr2=5.99
        do 30 j=1,jtot
          signif(j)=fft_theo(j)*chisqr2/dofmin
30     continue
       elseif(isigtest.eq.1)then
c----time-averaged, dof depend on scale----
c  see Eqn(23)
        do 40 j=1,jtot
          dof(j)=dofmin*sqrt(1+(n*dt/(gammafac*s(j)))**2)
          chi1=dof(j)+0.85+1.645*sqrt(2*dof(j)-1.0)
          signif(j)=fft_theo(j)*chi1/dof(j)
40     continue
      elseif(isigtest.eq.2)then
c----scale-averaged, dof depend on scale----
        javg1=0
        javg2=0
        do 50 j=1,jtot
          if((s(j).ge.dof(1)).and.(javg1.eq.0))javg1=j
          if(s(j).le.dof(2))javg2=j
50     continue
        if((javg1.eq.0).or.(javg2.eq.0).or.(javg1.gt.javg2))then
          print*,'***wave_signif: Scales in "dof(1)" & "dof(2)"'//
     &            'are out of range'
        return
        end if
        navg=javg2-javg1+1
c  see Eqn(25)
        savg=0.
        do 60 j=javg1,javg2
          savg=savg+1./s(j)
60     continue
        savg=savg+1./savg
c  see Eqn(27)
        fft_theor1=0.
        do 70 j=javg1,javg2
          fft_theor1=fft_theor1+fft_theo(j)/s(j)
70     continue
        fft_theo(1)=savg*fft_theor1
c  see Eqn(28)
        smid=exp(0.5*(alog(s(javg1))+alog(s(javg2))))
        dof(1)=(dofmin*navg*savg/smid)*sqrt(1+(navg*dj/dj0)**2)
        chi=dof(1)+0.85+1.645*sqrt(2*dof(1)-1.0)
c  see Eqn(26)
        signif(1)=(dj*dt/(cdelta*savg))*fft_theo(1)*chi/dof(1)
        do 80 j=2,jtot
          dof(j)=0.0
          fft_theo(j)=0.0
          signif(j)=0.0
80     continue
        endif
        end

river

这么长的程序，谁会没事儿一行行去看。有什么错误总有错误提示吧，总得有个方向吧。提问总得提到点子上吧，这样的问题不会有人想回答的。你还不如直接在论坛搜索小波分析，下载一个你觉得没有错的程序呢

wysyhd

river 发表于 2013-5-18 17:13

                               
登录/注册后可看大图

这么长的程序，谁会没事儿一行行去看。有什么错误总有错误提示吧，总得有个方向吧。提问总得提到点子上吧， ...

好吧，还是谢了，那你有没有错的程序？非常感谢

bestluyf

程序太长了，问问题还是具体点好吧

wysyhd

bestluyf 发表于 2013-5-18 18:39

                               
登录/注册后可看大图

程序太长了，问问题还是具体点好吧

   第8行有错，也就是下面这行
open(1,file='E:\wavelet\data.grd',form='binary',
     *recl=4)
麻烦帮忙解决下，很急谢谢了

bestluyf

wysyhd 发表于 2013-5-18 18:42

                               
登录/注册后可看大图

第8行有错，也就是下面这行
open(1,file='E:\wavelet\data.grd',form='binary',
     *recl=4)

不好意思啊，我不是很懂。应该不是这一行的问题吧。

wysyhd

bestluyf 发表于 2013-5-18 18:47

                               
登录/注册后可看大图

不好意思啊，我不是很懂。应该不是这一行的问题吧。

运行的结果是说那句有错，你再看看，谢谢啊

river

wysyhd 发表于 2013-5-18 18:42

                               
登录/注册后可看大图

第8行有错，也就是下面这行
open(1,file='E:\wavelet\data.grd',form='binary',
     *recl=4)

这句有错的话open(1,file='E:\wavelet\data.grd',form='binary',
     *recl=4)  那就把后面那个 *recl=4)都去掉吧，直接写 open(1,file='E:\wavelet\data.grd',form='binary')

wysyhd

river 发表于 2013-5-18 19:37

                               
登录/注册后可看大图

这句有错的话open(1,file='E:\wavelet\data.grd',form='binary',
     *recl=4)  那就把后面那个 *recl= ...

嗯，知道了，谢谢，我试试