      program four
      parameter(nx=49,ny=25,mx=73,my=37)
      parameter(p1=3.,p2=18.,p=1.,ndy=3796,l=nx*ny,n=ndy+p2)
      dimension wkin(l,n),wkout(l,n),tavg(l),avg(l),wk(l,15)
      dimension icount(l),indx(l)
      dimension wkout1(l,n),sst(nx,ny),sstm(nx,ny,12),var(mx,my)
      dmiss=9.999E+20
c---------------------------------------------------------
      open(14,file='/data2/NCEP/G_PEN_4899.gdat'
     *  ,status='old',access='direct',recl=mx*my,form='unformatted')
      open(20,file='/u1p/users/yklim/ncep/g200/G_bwfilt1590.gdat',
     *form='unformatted',status='unknown',access='direct',recl=nx*ny)
c---------------------------------------------------------
      irec=0
      do 10 i=1,ndy
       irec=(i-1)*5+5
       print*, i
      read(14,rec=irec) var
      ic=0
      do jj=1,ny
      do ii=1,nx
       ic=ic+1
       wkin(ic,i)=var(ii+12,jj+8)
      enddo
      enddo
  10  continue
c---------------------------------------------------------
      do 100 k=1,l
      print*,k
      avg(k)=0.
      tavg(k)=0.
      icount(k)=0
      do 20 j=1,ndy
      if (wkin(k,j) .ne. dmiss) then
       tavg(k)=tavg(k)+wkin(k,j)
       icount(k)=icount(k)+1
      endif
   20 continue
      avg(k)=tavg(k)/float(icount(k)) 
      do 30 j=1,ndy
      wkin(k,j)=wkin(k,j)-avg(k)
   30 continue
      do 40 j=ndy+1,n
      wkin(k,j)=0.
   40 continue
      call filtr4(p1,p2,p,wkin,wkout,wk,l,n)
      do 101 i1=1,n
  101 wkout1(k,i1)=wkout(k,n-i1+1)
      call filtr4(p1,p2,p,wkout1,wkout,wk,l,n)
  100 continue
      irec=0
      do 99 i=n,n-ndy+1,-1
      irec=irec+1
      ij=0
      do 400 jj=1,ny
      do 400 ii=1,nx
      ij=ij+1
      sst(ii,jj)=wkout(ij,i)
  400 continue
      write(20,rec=irec) sst
   99 continue
      stop
      end
c--------------------------------------------------
      subroutine filtr4(p1,p2,dt,ai,ao,wk,l,n)
c
c      fourth order butterworth band-pass filter
c
c      p1  lower limit of period range
c      p2  upper limit of period range (p1.lt.p2)
c      dt  time interval between two data points. dt should have
c           same units as p1 and p2
c      ai(l,n)  input time series
c      ao(l,n)  output time series
c      wk(l,15)  working space
c      l=nx*ny   number of time series to be filtered
c      n   length of each time series
c
c      before calling butter4 remove the time mean from
c      each series and extend the series by adding zeros
c      (usually the numer of zeros should correspond
c      roughly to p2). filter the extended series
c      and obtain the final series by retaining the first
c      n numbers only.
c
c      example 1: to filter two series of length 100 through
c      the period range 10 to 20 days. dt=2.0 days (every other day)
c      note that 10 data points correspond p2=20 days
c      in this case
c
c   example 2: for a time series of 100 days, with data each day
c   requesting a 10-20 day filter, use the following call
c   call filtr4(10.,20.,1.,a,b,wk,1,100) where
c   array to be filtered is in a, array returned is in b
c   main program must dimension wk to(l,15)
c   where l is the number of time series to be filtered
c
c      a(2,100) is the series to be filtered
c      ai(2,110),ao(2,110),wk(2,15)
c      p1=10.0
c      p2=20.0
c      dt=2.0
c      l=2
c      n=110
c
c      remove the mean
c
c      do 12 i=1,l
c      ave=0.0
c      do 11 kt=1,100
c   11 ave=ave+a(i,kt)/100.0
c      do 12 kt=1,100
c   12 ai(i,kt)=ai(i,kt)-ave
c
c      add zeros at the end
c
c      do 13 kt=101,110
c   13 ai(i,kt)=0.0
c      call filtr4(p1,p2,dt,ai,ao,l,n)
c
c      keep first 100 numbers only
c
c      do 14 i=1,l
c      do 14 kt=1,100
c   14 a(i,kt)=ao(i,kt)
c
c
      dimension ai(l,n),ao(l,n),wk(l,15)
      dimension aa(4),bb(3,4),cc(4),dd(3,4)
      complex a,b,c,x(4),y(2,4),z(4)
      pi=4.0*atan(1.0)
      w1=-2.0*cos(5.0*pi/8.0)
      w2=-2.0*cos(7.0*pi/8.0)
      f1=1.0/p2
      f2=1.0/p1
c
c   compute coefficients
c
      pi=4.0*atan(1.0)
      f=2.0*pi*sqrt(f1*f2)
      df=2.0*pi*(f2-f1)
      a=(1.0,0.0)
      b=cmplx(w1,0.0)
      c=cmplx(1.0,0.0)
      x(1)=0.5*(-b+csqrt(b*b-4.0*a*c))/a
      x(2)=0.5*(-b-csqrt(b*b-4.0*a*c))/a
      a=(1.0,0.0)
      b=cmplx(w2,0.0)
      c=cmplx(1.0,0.0)
      x(3)=0.5*(-b+csqrt(b*b-4.0*a*c))/a
      x(4)=0.5*(-b-csqrt(b*b-4.0*a*c))/a
      do 16 k=1,4
      a=(1.0,0.0)
      b=-cmplx(df,0.0)*x(k)
      c=cmplx(f*f,0.0)
      y(1,k)=0.5*(-b+csqrt(b*b-4.0*a*c))/a
      y(2,k)=0.5*(-b-csqrt(b*b-4.0*a*c))/a
   16 continue
      z(1)=y(1,1)
      z(2)=y(2,1)
      z(3)=y(1,3)
      z(4)=y(2,3)
      do 13 k=1,4
      cca=-2.0*real(z(k))
      dda=cabs(z(k))**2
      bb(1,k)=(4.0+2.0*dt*cca+dt*dt*dda)/(2.0*df*dt)
      bb(2,k)=(dt*dt*dda-4.0)/(df*dt)
      bb(3,k)=(4.0-2.0*dt*cca+dt*dt*dda)/(2.0*df*dt)
   13 continue
c
c   coef complete
c
      do 14 k=1,4
      kr=5-k
c     print *,bb(1,k),bb(2,k)/bb(1,k),bb(3,k)/bb(1,k)
      aa(k)=1.0/bb(1,k)
      bb(2,k)=bb(2,k)/bb(1,k)
      bb(3,k)=bb(3,k)/bb(1,k)
      cc(kr)=aa(k)
      dd(1,kr)=bb(1,k)
      dd(2,kr)=bb(2,k)
      dd(3,kr)=bb(3,k)
   14 continue
      do 11 k=1,15
      do 11 i=1,l
      wk(i,k)=0.0
   11 continue
      do 12 k=1,n
      do 12 i=1,l
      wk(i,1)=ai(i,k)
      wk(i,4)=aa(1)*(wk(i,1)-wk(i,3))-bb(2,1)*wk(i,5)-bb(3,1)*wk(i,6)
      wk(i,7)=aa(2)*(wk(i,4)-wk(i,6))-bb(2,2)*wk(i,8)-bb(3,2)*wk(i,9)
      wk(i,10)=aa(3)*(wk(i,7)-wk(i,9))
     1 -bb(2,3)*wk(i,11)-bb(3,3)*wk(i,12)
      wk(i,13)=aa(4)*(wk(i,10)-wk(i,12))
     1 -bb(2,4)*wk(i,14)-bb(3,4)*wk(i,15)
      wk(i,15)=wk(i,14)
      wk(i,14)=wk(i,13)
      wk(i,12)=wk(i,11)
      wk(i,11)=wk(i,10)
      wk(i,9)=wk(i,8)
      wk(i,8)=wk(i,7)
      wk(i,6)=wk(i,5)
      wk(i,5)=wk(i,4)
      wk(i,3)=wk(i,2)
      wk(i,2)=wk(i,1)
      ao(i,k)=wk(i,13)
   12 continue
      do 21 k=1,15
      do 21 i=1,l
      wk(i,k)=0.0
   21 continue
      do 22 kr=1,n
      k=n+1-kr
      do 22 i=1,l
      wk(i,1)=ao(i,k)
      wk(i,4)=aa(1)*(wk(i,1)-wk(i,3))-bb(2,1)*wk(i,5)-bb(3,1)*wk(i,6)
      wk(i,7)=aa(2)*(wk(i,4)-wk(i,6))-bb(2,2)*wk(i,8)-bb(3,2)*wk(i,9)
      wk(i,10)=aa(3)*(wk(i,7)-wk(i,9))
     1 -bb(2,3)*wk(i,11)-bb(3,3)*wk(i,12)
      wk(i,13)=aa(4)*(wk(i,10)-wk(i,12))
     1 -bb(2,4)*wk(i,14)-bb(3,4)*wk(i,15)
      wk(i,15)=wk(i,14)
      wk(i,14)=wk(i,13)
      wk(i,12)=wk(i,11)
      wk(i,11)=wk(i,10)
      wk(i,9)=wk(i,8)
      wk(i,8)=wk(i,7)
      wk(i,6)=wk(i,5)
      wk(i,5)=wk(i,4)
      wk(i,3)=wk(i,2)
      wk(i,2)=wk(i,1)
      ao(i,k)=wk(i,13)
   22 continue
      return
      end
