关于温度的滤波问题

Jadeite

我的数据为1979-2012共34年，lon 70 140，lat 15 55的气温资料（单位为K），希望得到10-30天的波段。下面贴出程序及滤波前后的对比：
1、程序：      parameter(nx=29,ny=17,nt=12419,mw1=10,mw2=30,nnt=nt)
      dimension r(nx,ny,nt),olr1(nx,ny,nnt),h(nx,ny,nnt),h1(nx,ny,nnt) ,olr(nx,ny,nnt)
      dimension cx(nx,ny,nt),xy(nx,ny,nt),bx(nx,ny,nt)
c-----------input data------------------
      open(11,file='1979-20120.grd',form='binary')
      do 10 k=1,nt
10  read(11)((r(i,j,k),i=1,nx),j=1,ny)
      close(11)
        write(*,*)'ok'
         do k=1,nnt
        do j=1,ny
        do i=1,nx
         olr(i,j,k)=r(i,j,k)
      enddo
        enddo
        enddo
      call filter(nt,nx,ny,olr,bx,cx,xy,mw1,mw2)
        do 60 j=1,ny
        do 60 i=1,nx
      do 60 it=1,nt
  60  olr1(i,j,it)=bx(i,j,it)


c----writing olr1(nx,ny,nt) to file 'olrarealb.grd'-----
      open(12,file='1979-20121.grd',form='binary')
        do  k=1,nt
     
      write(12)((olr1(i,j,k),i=1,nx),j=1,ny)
        enddo
      close(12)
       
      end
      subroutine filter(mt,lx,ly,x,bx,cx,xy,mw1,mw2)
      dimension x(lx,ly,mt),bx(lx,ly,mt),cx(lx,ly,mt),xy(lx,ly,mt)
      w1=2*3.1415926/float(mw1)
      w2=2*3.1415926/float(mw2)
      w0=sqrt(w1*w2)
      dt=1.0
      ww=2*abs(sin(w1*dt)/(1+cos(w1*dt))
     $ -sin(w2*dt)/(1+cos(w2*dt)))
      ww0=4*sin(w1*dt)*sin(w2*dt)/((1+cos(w1*dt))
     $ *(1+cos(w2*dt)))
      a0=2*ww/(4+2*ww+ww0)
      b1=2*(ww0-4)/(4+2*ww+ww0)
      b2=(4-2*ww+ww0)/(4+2*ww+ww0)
        do j=1,ly
        do ix=1,lx
      do 89 i0=1,mt-2
      i=i0+2
      bx(ix,j,1)=0.0
      bx(ix,j,2)=0.0
      cx(ix,j,i)=a0*(x(ix,j,i)-x(ix,j,i-2))-b1*bx(ix,j,i-1)-
     *        b2*bx(ix,j,i-2)
  89  bx(ix,j,i)=cx(ix,j,i)
      enddo
        enddo
        do j=1,ly
        do ix=1,lx
      do 91 i=1,mt
  91  x(ix,j,i)=bx(ix,j,mt+1-i)
       enddo
        enddo
        do j=1,ly
        do ix=1,lx


      do 92 i0=1,mt-2
      i=i0+2
      xy(ix,j,1)=x(ix,j,1)
      xy(ix,j,2)=x(ix,j,2)
      cx(ix,j,i)=a0*(x(ix,j,i)-x(ix,j,i-2))-b1*xy(ix,j,i-1)
     *-b2*xy(ix,j,i-2)
  92  xy(ix,j,i)=cx(ix,j,i)
       enddo
        enddo
        do j=1,ly
        do ix=1,lx
      do 93 i=1,mt
  93  bx(ix,j,i)=xy(ix,j,mt+1-i)
      enddo
        enddo
      end
2、滤波前后对比
前：

后：

3、数据ctl
dset E:\today\ERA-Interim\1979-20120.grd
undef -9.99E+8
title E:\today\ERA-Interim\1979-20120
ydef 17 linear 15 2.5
xdef 29 linear 70 2.5
tdef 12419 linear 00Z01jan1979 1dy
zdef 1 linear 1 1
vars 1
no2Tsfc  0 167,1,0  ** surface 2 metre temperature K
ENDVARS
求助各位来看看滤波的结果对不对

小亭子ng

楼主，你的问题解决了没啊，我最近也在处理butterworth滤波，是对夏季高温资料进行直接滤波，结果出来有负值，想问楼主是直接对资料滤波还是距平后滤波的，结果有负值是怎么回事啊

stone

subroutine filter子程序的滤波公式能解释下吗