      subroutine GrdLowpFilter(dt,hd,nlat,nlon,rst,knd,nm)
      !对格网dt(nlat,nlon)进行低通滤波
      !scale-滤波尺度(分)
      !knd-0滑动平均滤波法,1高斯低通滤波法,2指数型低通滤波,3巴特沃斯滤波法
      !nm-指数大小1~10。当knd=1时无效。
!-------------------------------------------------------------
      implicit none
	include 'fftw3.f'
	integer::nlat,nlon,knd,mm,nm,i,j,ki,kj,kk,astat(8)
	integer::row,col,lp,rp,scale,ii,jj
	real*8::hd(6),dt(nlat,nlon),rst(nlat,nlon),pi,RAD,val,ds,dn
	complex*16,allocatable::cgr(:,:),fgr(:,:),tmp(:,:)
      integer*8::fwd=0,bwd=0
!-----------------------------------------------------------------------
  	pi=datan(1.d0)*4.d0;RAD=pi/180.d0
      if(knd==0)then !滑动平均滤波
	  do i=1,nlat
	    do j=1,nlon
	      kk=0
	      do ki=i-nm,i+nm
	        do kj=j-nm,j+nm
	          if(ki<1.or.ki>nlat.or.kj<1.or.kj>nlon)goto 9101
	          if(dt(ki,kj)>9900.d0)goto 9101
	          rst(i,j)=rst(i,j)+dt(ki,kj);kk=kk+1
9101	          continue
	        enddo
	      enddo
	      if(kk>0)rst(i,j)=rst(i,j)/real(kk)
	      if(kk==0)rst(i,j)=9999.d0
	    enddo
        enddo
        return
      endif
	lp=nlat/2;rp=nlon/2
	row=nlat+2*lp; col=nlon+2*rp
 	allocate(cgr(row,col), stat=astat(1))
 	allocate(fgr(row,col), stat=astat(2))
 	allocate(tmp(row,col), stat=astat(3))
	if (sum(astat(1:3)) /= 0) then
          rst=dt;goto 9002
	endif
      call dfftw_plan_dft_2d(fwd,row,col,cgr,tmp,FFTW_FORWARD,FFTW_ESTIMATE)
      call dfftw_plan_dft_2d(bwd,row,col,fgr,tmp,FFTW_BACKWARD,FFTW_ESTIMATE)
	cgr=0.d0;cgr(lp+1:lp+nlat,rp+1:rp+nlon)=dt
      call dfftw_execute(fwd)
      cgr=0.d0
	do i=1,row
	   do j=1,col
	      if(knd==1)val=exp(-real((i-1)**2+(j-1)**2)/real(row*col)/2.d0)
	      if(knd==2)val=exp(-(real((i-1)**2+(j-1)**2)/sqrt(sqrt(real(row*col)**nm))))
	      if(knd>0.and.knd<3)cgr(i,j)=val
            if(knd==3)then
	        val=sqrt(real((i-1)**2+(j-1)**2)/real(row*col))
	        cgr(i,j)=1.d0/(1.d0+val**nm)
            endif
         enddo
	enddo
	fgr=tmp*cgr
      call dfftw_execute(bwd)
	rst=real(tmp(lp+1:lp+nlat,rp+1:rp+nlon))/dble(row*col)
      call dfftw_destroy_plan(fwd)
      call dfftw_destroy_plan(bwd)
	deallocate(cgr,fgr,tmp)
9002	return
      end
