对OLR的逐日资料做10-90天的滤波

stefan_scofield

刚开始上手做毕业论文，老师说要对OLR的逐日资料做10-90天的滤波，我一窍不通，找来了一段fortran程序，我看不大懂，自己改了主程序并尝试运行多次，出来的结果都是NaN，是不是子程序运算的有错，以为子程序我实在看不懂，希望会的同学老师们提供帮助。底下是程序：
      program mian
      parameter(nx=72,ny=32,nt=5202,freqL=1/90.,freqU=1/10.)

      dimension olr(nx,ny,nt)
        dimension r2(nx,ny,nt)
        dimension mid(nx,ny,nt)
      dimension olrf(nx,ny,nt)
!      complex r1(nt),work(nt)
      real r1(nt),work(nt)

      OPEN(11,FILE='olr.79-12.6-10.grd',FORM='BINARY')
      PRINT*,'NOW READING ...'
      DO IT=1,NT
      DO J=1,NY
      DO I=1,NX
        READ(11) olr(I,J,IT)
      ENDDO
      ENDDO
      enddo
      close(11)


        do it=1,nt
        do j=1,ny
        do i=1,nx
                mid(i,j,it)=olr(i,j,it)
        enddo
        enddo
        enddo

      do 111 i=1,nx
      do 111 j=1,ny
        do 50 k=1,nt
         r1(k)=mid(i,j,k)
   50   continue
        call bandps(nt,r1,1.,freqL,freqU,0,work)
        do 60 k=1,nt
        r2(i,j,k)=real(r1(k))
   60   continue
  111   continue
      do k=1,nt
      do j=1,ny
      do i=1,nx
         olrf(i,j,k)=r2(i,j,k)
      enddo
      enddo
      enddo



      OPEN(12,FILE='bp-filtered-olr.grd',FORM='BINARY')
      PRINT*,'NOW WRITING ...'
      do k=1,nt
      do j=1,ny
      do i=1,nx
        if(abs(olrf(i,j,k))>5000.0) then
         olrf(i,j,k)=-9999
           write(12) olrf(i,j,k)
          else
           write(12) olrf(i,j,k)
          endif
      enddo
      enddo
      enddo
      close(12)
      print*,'NOW COMPLETED!'

      stop
      end

      subroutine bandps (npts, ts, dt, freqL, freqU, iform, work)
      integer  npts, iform
      real     dt, freqL, freqU
    !  complex  ts(npts), work(npts)
       real  ts(npts), work(npts)
      integer  j, jL, jU
      if (iform .lt. 0 .or. iform .gt. 1)
     x   stop 'BANDPS: iform is out of range.'

      jL = int(freqL*dt*npts + 0.5) + 1
      if(freqL .lt. 0) jL=-int(-freqL*dt*npts + 0.5)+1
c
      jU = int(freqU*dt*npts + 0.5) + 1
      if(freqU .lt. 0) jU=-int(-freqU*dt*npts + 0.5)+1
c
c     If iform = 0, freqL and/or freqU are unsuitable if freqL .lt. 0,
c     freqU .lt. freqL, and/or freqU .gt. 1/(2*dt).
      if (iform .eq. 0  .and.
     x   (jL .lt. 1  .or.  jU .lt. jL  .or.  jU .gt. 1 + npts/2))
     x   stop 'BANDPS: freqL and/or freqU are out of range.'
c
c     If iform = 1, freqL and/or freqU are unsuitable if
c     freqL .lt. -1/(2*dt), freqU .lt. freqL, and/or freqU .gt. 1/dt.
      if (iform .eq. 1  .and.
     x   (jL .lt. 1-npts/2  .or.  jU .lt. jL  .or.  jU .gt. npts))
     x   stop 'BANDPS: freqL and/or freqU are out of range.'


c
c-----
c Fourier transform the complex time series.
      call fourt (ts, npts, 1, -1, iform, work)
c
c     Normalize the transform so that the values represent
c     wave amplitudes.
      do 10 j  = 1, npts
         ts(j) = ts(j)/npts
   10 continue


c
c-----
c Apply the frequency filter.
c     Ordinates indexed by jL through jU, inclusive, will be retained.
c     All other ordinates will be set to zero.
      if(jL .le. 0 .and. jU .ge. 1) then
         do 30 j  = jU+1, jL+npts-1, 1
            ts(j) = 0.0
   30    continue
      else
c        This covers the case (jL.le.0 .and. jU.le.0)
c        .or. (jL.ge.1 .and. jU.ge.1).
         if (jL .le. 0) then
            jL = jL + npts
            jU = jU + npts
         endif
         do 130 j  = 1, jL-1, 1
            ts (j) = 0.0
  130    continue
         if (iform .eq. 0) then
c           The time series is real (the imaginary part is zero),
c           so we want to retain both frequency ranges (-freqU, -freqL)
c           and (freqL, freqU) so that the inverse transform
c           will also be real.
            do 150 j  = jU+1, npts+1-jU, 1
               ts (j) = 0.0
  150       continue
            do 160 j  = npts+3-jL, npts, 1
               ts (j) = 0.0
  160       continue
         else
c           The time series is complex, so we want to retain only the
c           frequency range (freqL, freqU).
            do 180 j  = jU+1, npts, 1
               ts (j) = 0.0
  180       continue
         endif
      endif
c
c-----
c Compute the inverse transform.
      call fourt (ts, npts, 1, 1, 1, work)

      return
      end


*********************************************
      subroutine undefin(rain,nx,ny,nt)
      real rain(nx,ny,nt)

      do 21 k=1,nt
      do 20 j=1,ny
      do 20 i=1,nx

      if(rain(i,j,k).lt.0.0)then
      rain(i,j,k)=-9999
      endif

20    continue
21    continue  

      return
      end


*********************************************
      subroutine fourt (data, nn, ndim, isign, iform, work)
      real  data(*), work(*)
      integer  nn(*),   ifact(32)
      data twopi/6.2831853071796/,rthlf/0.70710678118655/
c     twopi = 2 pi, rthlf = sqrt(1./2.)
c
c     Check input parameters for propriety.  Return if non-positive
c     dimensions are encountered.
      if(ndim-1)920,1,1
1     ntot=2
      do 2 idim=1,ndim
      if(nn(idim))920,920,2
2     ntot=ntot*nn(idim)
      np1=2
      do 910 idim=1,ndim
      n=nn(idim)
      np2=np1*n
      if(n-1)920,900,5
c

5     m=n
      ntwo=np1
      if=1
      idiv=2
10    iquot=m/idiv
      irem=m-idiv*iquot
      if(iquot-idiv)50,11,11
11    if(irem)20,12,20
12    ntwo=ntwo+ntwo
      ifact(if)=idiv
      if=if+1
      m=iquot
      go to 10
20    idiv=3
      inon2=if
30    iquot=m/idiv
      irem=m-idiv*iquot
      if(iquot-idiv)60,31,31
31    if(irem)40,32,40
32    ifact(if)=idiv
      if=if+1
      m=iquot
      go to 30
40    idiv=idiv+2
      go to 30
50    inon2=if
      if(irem)60,51,60
51    ntwo=ntwo+ntwo
      go to 70
60    ifact(if)=m
c
c
70    icase=1
      ifmin=1
      i1rng=np1
      if(idim-4)71,100,100
71    if(iform)72,72,100
72    icase=2
      i1rng=npo*(1+nprev/2)
      if(idim-1)73,73,100
73    icase=3
      i1rng=np1
      if(ntwo-np1)100,100,74
74    icase=4
      ifmin=2
      ntwo=ntwo/2
      n=n/2
      np2=np2/2
      ntot=ntot/2
      i=1
      do 80 j=1,ntot
      data(j)=data(i)
80    i=i+2
c
c
100   if(ntwo-np2)200,110,110
110   np2hf=np2/2
      j=1
      do 150 i2=1,np2,np1
      if(j-i2)120,130,130
120   i1max=i2+np1-2
      do 125 i1=i2,i1max,2
      do 125 i3=i1,ntot,np2
      j3=j+i3-i2
      tempr=data(i3)
      tempi=data(i3+1)
      data(i3)=data(j3)
      data(i3+1)=data(j3+1)
      data(j3)=tempr
125   data(j3+1)=tempi
130   m=np2hf
140   if(j-m)150,150,145
145   j=j-m
      m=m/2
      if(m-np1)150,140,140
150   j=j+m
      go to 300
c

200   nwork=2*n
      do 270 i1=1,np1,2
      do 270 i3=i1,ntot,np2
      j=i3
      do 260 i=1,nwork,2
      if(icase-3)210,220,210
210   work(i)=data(j)
      work(i+1)=data(j+1)
      go to 230
220   work(i)=data(j)
      work(i+1)=0.
230   ifp2=np2
      if=ifmin
240   ifp1=ifp2/ifact(if)
      j=j+ifp1
      if(j-i3-ifp2)260,250,250
250   j=j-ifp2
      ifp2=ifp1
      if=if+1
      if(ifp2-np1)260,260,240
260   continue
      i2max=i3+np2-np1
      i=1
      do 270 i2=i3,i2max,np1
      data(i2)=work(i)
      data(i2+1)=work(i+1)
270   i=i+2
c

300   if(ntwo-np1)600,600,305
305   np1tw=np1+np1
      ipar=ntwo/np1
310   if(ipar-2)350,330,320
320   ipar=ipar/4
      go to 310
330   do 340 i1=1,i1rng,2
      do 340 k1=i1,ntot,np1tw
      k2=k1+np1
      tempr=data(k2)
      tempi=data(k2+1)
      data(k2)=data(k1)-tempr
      data(k2+1)=data(k1+1)-tempi
      data(k1)=data(k1)+tempr
340   data(k1+1)=data(k1+1)+tempi
350   mmax=np1
360   if(mmax-ntwo/2)370,600,600
370   lmax=max0(np1tw,mmax/2)
      do 570 l=np1,lmax,np1tw
      m=l
      if(mmax-np1)420,420,380
380   theta=-twopi*float(l)/float(4*mmax)
      if(isign)400,390,390
390   theta=-theta
400   wr=cos(theta)
      wi=sin(theta)
410   w2r=wr*wr-wi*wi
      w2i=2.*wr*wi
      w3r=w2r*wr-w2i*wi
      w3i=w2r*wi+w2i*wr
420   do 530 i1=1,i1rng,2
      kmin=i1+ipar*m
      if(mmax-np1)430,430,440
430   kmin=i1
440   kdif=ipar*mmax
450   kstep=4*kdif
      if(kstep-ntwo)460,460,530
460   do 520 k1=kmin,ntot,kstep
      k2=k1+kdif
      k3=k2+kdif
      k4=k3+kdif
      if(mmax-np1)470,470,480
470   u1r=data(k1)+data(k2)
      u1i=data(k1+1)+data(k2+1)
      u2r=data(k3)+data(k4)
      u2i=data(k3+1)+data(k4+1)
      u3r=data(k1)-data(k2)
      u3i=data(k1+1)-data(k2+1)
      if(isign)471,472,472
471   u4r=data(k3+1)-data(k4+1)
      u4i=data(k4)-data(k3)
      go to 510
472   u4r=data(k4+1)-data(k3+1)
      u4i=data(k3)-data(k4)
      go to 510
480   t2r=w2r*data(k2)-w2i*data(k2+1)
      t2i=w2r*data(k2+1)+w2i*data(k2)
      t3r=wr*data(k3)-wi*data(k3+1)
      t3i=wr*data(k3+1)+wi*data(k3)
      t4r=w3r*data(k4)-w3i*data(k4+1)
      t4i=w3r*data(k4+1)+w3i*data(k4)
      u1r=data(k1)+t2r
      u1i=data(k1+1)+t2i
      u2r=t3r+t4r
      u2i=t3i+t4i
      u3r=data(k1)-t2r
      u3i=data(k1+1)-t2i
      if(isign)490,500,500
490   u4r=t3i-t4i
      u4i=t4r-t3r
      go to 510
500   u4r=t4i-t3i
      u4i=t3r-t4r
510   data(k1)=u1r+u2r
      data(k1+1)=u1i+u2i
      data(k2)=u3r+u4r
      data(k2+1)=u3i+u4i
      data(k3)=u1r-u2r
      data(k3+1)=u1i-u2i
      data(k4)=u3r-u4r
520   data(k4+1)=u3i-u4i
      kdif=kstep
      kmin=4*(kmin-i1)+i1
      go to 450
530   continue
      m=m+lmax
      if(m-mmax)540,540,570
540   if(isign)550,560,560
550   tempr=wr
      wr=(wr+wi)*rthlf
      wi=(wi-tempr)*rthlf
      go to 410
560   tempr=wr
      wr=(wr-wi)*rthlf
      wi=(tempr+wi)*rthlf
      go to 410
570   continue
      ipar=3-ipar
      mmax=mmax+mmax
      go to 360
c
c
600   if(ntwo-np2)605,700,700
605   ifp1=ntwo
      if=inon2
      np1hf=np1/2
610   ifp2=ifact(if)*ifp1
      j1min=np1+1
      if(j1min-ifp1)615,615,640
615   do 635 j1=j1min,ifp1,np1
      theta=-twopi*float(j1-1)/float(ifp2)
      if(isign)625,620,620
620   theta=-theta
625   wstpr=cos(theta)
      wstpi=sin(theta)
      wr=wstpr
      wi=wstpi
      j2min=j1+ifp1
      j2max=j1+ifp2-ifp1
      do 635 j2=j2min,j2max,ifp1
      i1max=j2+i1rng-2
      do 630 i1=j2,i1max,2
      do 630 j3=i1,ntot,ifp2
      tempr=data(j3)
      data(j3)=data(j3)*wr-data(j3+1)*wi
630   data(j3+1)=tempr*wi+data(j3+1)*wr
      tempr=wr
      wr=wr*wstpr-wi*wstpi
635   wi=tempr*wstpi+wi*wstpr
640   theta=-twopi/float(ifact(if))
      if(isign)650,645,645
645   theta=-theta
650   wstpr=cos(theta)
      wstpi=sin(theta)
      j2rng=ifp1*(1+ifact(if)/2)
      do 695 i1=1,i1rng,2
      do 695 i3=i1,ntot,np2
      j2max=i3+j2rng-ifp1
      do 690 j2=i3,j2max,ifp1
      j1max=j2+ifp1-np1
      do 680 j1=j2,j1max,np1
      j3max=j1+np2-ifp2
      do 680 j3=j1,j3max,ifp2
      jmin=j3-j2+i3
      jmax=jmin+ifp2-ifp1
      i=1+(j3-i3)/np1hf
      if(j2-i3)655,655,665
655   sumr=0.
      sumi=0.
      do 660 j=jmin,jmax,ifp1
      sumr=sumr+data(j)
660   sumi=sumi+data(j+1)
      work(i)=sumr
      work(i+1)=sumi
      go to 680
665   iconj=1+(ifp2-2*j2+i3+j3)/np1hf
      j=jmax
      sumr=data(j)
      sumi=data(j+1)
      oldsr=0.
      oldsi=0.
      j=j-ifp1
670   tempr=sumr
      tempi=sumi
      sumr=twowr*sumr-oldsr+data(j)
      sumi=twowr*sumi-oldsi+data(j+1)
      oldsr=tempr
      oldsi=tempi
      j=j-ifp1
      if(j-jmin)675,675,670
675   tempr=wr*sumr-oldsr+data(j)
      tempi=wi*sumi
      work(i)=tempr-tempi
      work(iconj)=tempr+tempi
      tempr=wr*sumi-oldsi+data(j+1)
      tempi=wi*sumr
      work(i+1)=tempr+tempi
      work(iconj+1)=tempr-tempi
680   continue
      if(j2-i3)685,685,686
685   wr=wstpr
      wi=wstpi
      go to 690
686   tempr=wr
      wr=wr*wstpr-wi*wstpi
      wi=tempr*wstpi+wi*wstpr
690   twowr=wr+wr
      i=1
      i2max=i3+np2-np1
      do 695 i2=i3,i2max,np1
      data(i2)=work(i)
      data(i2+1)=work(i+1)
695   i=i+2
      if=if+1
      ifp1=ifp2
      if(ifp1-np2)610,700,700
c
c
700   go to (900,800,900,701),icase
701   nhalf=n
      n=n+n
      theta=-twopi/float(n)
      if(isign)703,702,702
702   theta=-theta
703   wstpr=cos(theta)
      wstpi=sin(theta)
      wr=wstpr
      wi=wstpi
      imin=3
      jmin=2*nhalf-1
      go to 725
710   j=jmin
      do 720 i=imin,ntot,np2
      sumr=(data(i)+data(j))/2.
      sumi=(data(i+1)+data(j+1))/2.
      difr=(data(i)-data(j))/2.
      difi=(data(i+1)-data(j+1))/2.
      tempr=wr*sumi+wi*difr
      tempi=wi*sumi-wr*difr
      data(i)=sumr+tempr
      data(i+1)=difi+tempi
      data(j)=sumr-tempr
      data(j+1)=-difi+tempi
720   j=j+np2
      imin=imin+2
      jmin=jmin-2
      tempr=wr
      wr=wr*wstpr-wi*wstpi
      wi=tempr*wstpi+wi*wstpr
725   if(imin-jmin)710,730,740
730   if(isign)731,740,740
731   do 735 i=imin,ntot,np2
735   data(i+1)=-data(i+1)
740   np2=np2+np2
      ntot=ntot+ntot
      j=ntot+1
      imax=ntot/2+1
745   imin=imax-2*nhalf
      i=imin
      go to 755
750   data(j)=data(i)
      data(j+1)=-data(i+1)
755   i=i+2
      j=j-2
      if(i-imax)750,760,760
760   data(j)=data(imin)-data(imin+1)
      data(j+1)=0.
      if(i-j)770,780,780
765   data(j)=data(i)
      data(j+1)=data(i+1)
770   i=i-2
      j=j-2
      if(i-imin)775,775,765
775   data(j)=data(imin)+data(imin+1)
      data(j+1)=0.
      imax=imin
      go to 745
780   data(1)=data(1)+data(2)
      data(2)=0.
      go to 900
c
c
800   if (i1rng-np1)805,900,900
805   do 860 i3=1,ntot,np2
      i2max=i3+np2-np1
      do 860 i2=i3,i2max,np1
      imin=i2+i1rng
      imax=i2+np1-2
      jmax=2*i3+np1-imin
      if(i2-i3)820,820,810
810   jmax=jmax+np2
820   if(idim-2)850,850,830
830   j=jmax+npo
      do 840 i=imin,imax,2
      data(i)=data(j)
      data(i+1)=-data(j+1)
840   j=j-2
850   j=jmax
      do 860 i=imin,imax,npo
      data(i)=data(j)
      data(i+1)=-data(j+1)
860   j=j-npo
c
c
  900 npo=np1
      np1=np2
910   nprev=n
920   return
      end

大气Tan

{:eb302:}求大神

lqouc

这老长一大段，你自己又没有提供什么有用的信息，估计不会有人给你看的。
建议自己先检查，从原始数据，到数据读取一步步的检查，一般子程序是不会出错的。
如果你自己实在解决不了，至少先把问题大概定位，贴出原因，大家再帮你。

stefan_scofield

lqouc 发表于 2014-3-21 14:16

                               
登录/注册后可看大图

这老长一大段，你自己又没有提供什么有用的信息，估计不会有人给你看的。
建议自己先检查，从原始数据，到 ...

我就是不知道哪里错了，改了运行试了好多次，除了报错，输出的都是NaN。主要是子程序长，那请您帮我看下主程序有没有错吧

stefan_scofield

主要是子程序长，大家先可忽略子程序，先看下主程序有没有问题

chaijing27

楼主，你好，不知道你的问题有没有解决，解决的了的话求指导啊，我也遇到了同样的问题。