	parameter ( nf=15,nft=2*nf,nfp=nf+1,di=3.,dl=10.,dit=1. )
        parameter ( nx=144,ny=72,nfile=1 )
	parameter ( ndata=365 )
	dimension c(nfp),b(nf),d(nx,ny,nft+1),xx(nx,ny),xf(nx,ny),
     +   iday(nft+1)
	dimension nc(nx,ny),sum(nx,ny),iys(nfile),iye(nfile)
	character*20 fname(nfile),pname(nfile)
	data dmiss/99.900/
      open(1,file='../ts.hcl.85-89.gdat',status='old',
     +  form='unformatted',access='direct',recl=nx*ny)
      open(2,file='hcloud.ts.dl',status='unknown',
     +  form='unformatted',access='direct',recl=nx*ny)
c      open(7,file='',form='unformatted',
c     +    access='direct',recl=nx*ny,status='old')

        call haming(di,dl,dit,0.01,nf,c)

	b0=c(1)
	crtl=b0*0.1
	do 22 i=1,nf
	b(i)=c(i+1)
22	continue

	do 1000 nfl=1,nfile
	do 5 i=1,nx
        do 5 j=1,ny
	sum(i,j)=0.
	nc(i,j)=0
5	continue

	irec=0
	do 6 i=1,ndata
        irec=irec+1
        read(1,rec=irec)xx

	do 7 j=1,nx
        do 7 k=1,ny
	if (xx(j,k).ne.dmiss) then
	nc(j,k)=nc(j,k)+1
	sum(j,k)=sum(j,k)+xx(j,k)
	endif
7	continue
6	continue
	do 8 j=1,nx
        do 8 k=1,ny
	if (nc(j,k).ne.0) then
	sum(j,k)=sum(j,k)/float(nc(j,k))
	else
	sum(j,k)=dmiss
	endif
8	continue

	irec=0
	do 10 i=1,nft
        irec=irec+1
c	read(7,rec=irec) iday(i+1)
        read(1,rec=irec) xx
	do 20 j=1,nx
        do 20 k=1,ny
	if (xx(j,k).ne.dmiss) then
	d(j,k,i+1)=xx(j,k)-sum(j,k)
	else
	d(j,k,i+1)=dmiss
	endif
20	continue
10	continue
c
	do 888 icon=1,ndata-nft
c	read(7,rec=irec)id
        irec=irec+1
        read(1,rec=irec)xx

	do 30 i=1,nft
c	iday(i)=iday(i+1)
	do 30 j=1,nx
        do 30 k=1,ny
	d(j,k,i)=d(j,k,i+1)
30	continue
c
	do 40 j=1,nx
        do 40 k=1,ny
        if (xx(j,k).ne.dmiss) then
        d(j,k,nft+1)=xx(j,k)-sum(j,k)
        else
        d(j,k,nft+1)=dmiss
        endif
c	iday(nft+1)=id
40	continue
c 	idate=iday(nfp)
c
	do 45 j=1,nx
        do 45 k=1,ny
	xf(j,k)=0.
45	continue
c
	do 50 j=1,nx
        do 50 k=1,ny
	if (d(j,k,nfp).eq.dmiss) then
	xf(j,k)=dmiss
 	goto 50
	endif
	xf(j,k)=xf(j,k)+b0*d(j,k,nfp)
	do 60 i=1,nf
	dm1=d(j,k,nfp-i)
	dp1=d(j,k,nfp+i)
	if (d(j,k,nfp-i).eq.dmiss.and.abs(b(i)).gt.crtl) then
	xf(j,k)=dmiss
	goto 50
	endif
	if (d(j,k,nfp+i).eq.dmiss.and.abs(b(i)).gt.crtl) then
        xf(j,k)=dmiss
        goto 50
        endif
	if (dm1.eq.dmiss) dm1=0.
	if (dp1.eq.dmiss) dp1=0.
	xf(j,k)=xf(j,k)+b(i)*( dm1+dp1 )
60	continue
50	continue
	write(2,rec=icon) xf
888	continue
1000	continue
	stop
	end
cc
	subroutine haming(tb,ta,dt,delta,n,c)
c	tb must be less than  ta
cref) hamming,r.w.,1983; digital filter , ch.9  kaiser windows & opimization
	dimension c(n+1),h(1000,1)
c	character*4 name(1)
	data pi/3.141592653/
c
	nfreq=ifix(1./delta/2.)+1
	fb=dt/tb
	fa=dt/ta
	deltaf=dt/(fa-fb)
	a=-alog10(delta)*20.
	newn=ifix((a-7.95)/deltaf/28.72+0.5)
	if (n.lt.newn) write(*,11)
11	format(1x,'$ you had better change your "n" value :: ')
	df=(a-7.95)/float(n)/28.72
	if(a.le.21.) alpha=0.
	if(a.gt.21 .and. a.lt.50.) alpha=0.5842*(a-21.)**(0.4)+
     &	0.07886*(a-21.)
	if(a.ge.50.) alpha=0.1102*(a-8.7)
	const=func(alpha,0,n)
	c(1)=2.*(fb-fa)
	do 10 k=1,n
	angle1=2.*pi*float(k)*fb
	angle2=2.*pi*float(k)*fa       
	c(k+1)=(sin(angle1)-sin(angle2))/float(k)/pi
	c(k+1)=c(k+1)*func(alpha,k,n)/const
10	continue
	nn=n+1
c	do 20 k=1,nn
c	kk=k-1
c	write(6,2) kk,c(k)
c20	continue
c2	format(10x,'c(',i2,')=',f20.10)
	f=-1.*delta
	do 30 j=1,nfreq
	f=f+delta
	h(j,1)=c(1)
	do 25 i=1,n
	const=float(i)
	angle=2.*pi*const*f
	h(j,1)=h(j,1)+2.*c(i+1)*cos(angle)
25	continue
c	write(7,rec=j) h(j,1)
30	continue
	return
	end
c
c
	function func(alpha,k,n)
	y=1.-(float(k)/float(n))**2
	x=alpha*sqrt(y)/2.
	func=1.
	do 10 j=1,20
	fac=1.
	do 9 i=1,j
	fac=fac*float(i)
9	continue
	func=func+(x**j/fac)**2
10	continue
	return
	end
