c   subroutine 'cfftfil' removes flagged values from grid, transforms
c   data, calls routines to apply filter, and writes filtered grid
c***********************************************************************
          subroutine cfftfil(a1[huge],work,nx,nri,nxa,n1,n2,id2)
          dimension a1(2,n2,nri),work(id2),jflag(2,20)
      character crfile*12
          character*50 coname,fname
       common/bmain[near]/id(14),pgm(2),nz,yo,xo,dx,dy,iw,kr,ny,nadd,ir
       common/parm1[near]/fname,coname,zlev(25),nlev,inocal,barlev,ifac
          data iunita/10/,iunitb/12/,dval/1.701412E+38/
          data pi2/6.28318531/,pi/3.141592654/,iunitc/15/
      nlevel=1
      z=barlev-zlev(1)
      coname='mchess10.tmp'
      w1=0.
      w2=0.
      if(ifac.lt.0) go to 2
      del=(dx+dy)/2.
      flim=del*2.5
      fac=(ifac+1.5)*del
      wavlen=z*1.25+fac
      if(wavlen.lt.flim)wavlen=flim
      wdel=.5*del+0.25*z
      w1=wavlen-wdel
      w2=wavlen+wdel
      f4=1./w1
      f3=1./w2
      btan34=1./(f4-f3)
      go to 4
2     f4=0.
4         ncase=0
      n22=2*n2
c  compute beginning and end column & row in original grid.
   20     do 30 i=1,2
          do 30 j=1,n2
   30     a1(2,j,i)=0.0
          nap1=nadd+1
          n2a=n2-nadd
          nxap1=nxa*0.5+1
          fxap1=float(nxa)/2.+1.
          if(abs(fxap1-float(nxap1)).gt.0.0001)nxap1=nxap1+1
          n1a=nxap1-1+nx
c  remove flagged values and add '2*nadd' rows & columns to
c  reduce the effect of gibbs phenomena.
         open(10,access='direct',status='unknown',form='unformatted',
     1   file='slave1.tmp',recl=n22*4)
          open(14,status='unknown',form='unformatted',file='flag.loc')
          irr=1
          iwr=0
          jr=0
   40     jr=jr+1
          if(jr.gt.1) go to 140
c   read 1st row & remove flagged values.
          read(kr)dum,(a1(1,i,1),i=nap1,n2a)
          read(kr)dum,(a1(1,i,2),i=nap1,n2a)
          nflag=1
          i=nap1
   50     jflag(1,1)=nap1
          if(a1(1,nap1,1).le.1.0e38) go to 90
          do 60 i=nap1,n2a
   60     if(a1(1,i,1).le.1.0e38) go to 70
          go to 1270
   70     jflag(1,1)=i
          aval=a1(1,i,1)
          do 80 ii=nap1,i-1
   80     a1(1,ii,1)=aval
   90     if(i.ge.n2a) go to 101
          do 100 ii=i+1,n2a
          if(a1(1,ii,1).le.1.0e38) go to 100
          jflag(2,nflag)=ii-1
          go to 110
  100     continue 
101     continue
          jflag(2,nflag)=n2a
          write(14)nflag,(jflag(1,j),jflag(2,j),j=1,nflag)
          go to 350
  110     aval=a1(1,ii-1,1)
          do 120 i=ii,n2a
          if(a1(1,i,1).le.1.0e38) go to 130
  120     a1(1,i,1)=aval
          write(14)nflag,(jflag(1,j),jflag(2,j),j=1,nflag)
          go to 350
  130     a1(1,i-1,1)=(a1(1,i-1,1)+a1(1,i,1))*0.5
          nflag=nflag+1
          jflag(1,nflag)=i
c   1st row completed: now extend grid & write to disk.
          go to 90
c   remove flagged values from remaining rows.
  140     if(jr.ge.nx) go to 150
          read(kr)dum,(a1(1,i,3),i=nap1,n2a)
          go to 170
  150     do 160 i=nap1,n2a
  160     a1(1,i,3)=dval
  170     nflag=1
          i=nap1
          jflag(1,1)=nap1
          if(a1(1,nap1,2).le.1.0e38) go to 220
          do 180 i=nap1,n2a
  180     if(a1(1,i,2).le.1.0e38) go to 190
          i=nap1
          jflag(1,1)=n2a+1
          if(a1(1,nap1,3).ge.1.0e38)a1(1,nap1,2)=a1(1,nap1,1)
          if(a1(1,nap1,3).le.1.0e38)a1(1,nap1,2)=(a1(1,nap1,1)+
     1    a1(1,nap1,3))*0.5
          go to 220
  190     jflag(1,1)=i
          ii=i
          do 210 irev=nap1,i-1
          ii=ii-1
          if(a1(1,ii,3).ge.1.0e38) go to 200
          a1(1,ii,2)=(a1(1,ii,1)+a1(1,ii+1,2)+a1(1,ii,3))*0.3333333
          go to 210
  200     a1(1,ii,2)=(a1(1,ii,1)+a1(1,ii+1,2))*0.5
  210     continue
  220     if(i.ge.n2a) go to 231
          do 230 ii=i+1,n2a
          if(a1(1,ii,2).le.1.0e38) go to 230
          jflag(2,nflag)=ii-1
          go to 240
  230     continue  
231     continue
          jflag(2,nflag)=n2a
          write(14)nflag,(jflag(1,j),jflag(2,j),j=1,nflag)
          go to 390
  240     do 290 i=ii,n2a
          if(i.eq.n2a) go to 250
          if(a1(1,i+1,2).le.1.0e38) go to 270
  250     if(a1(1,i,3).ge.1.0e38) go to 260
          a1(1,i,2)=(a1(1,i,1)+a1(1,i-1,2)+a1(1,i,3))*0.3333333
          go to 290
  260     a1(1,i,2)=(a1(1,i,1)+a1(1,i-1,2))*0.5
          go to 290
  270     if(a1(1,i,3).ge.1.0e38) go to 280
          a1(1,i,2)=(a1(1,i,1)+a1(1,i-1,2)+a1(1,i+1,2)+a1(1,i,3))*0.25
          go to 300
  280     a1(1,i,2)=(a1(1,i,1)+a1(1,i-1,2)+a1(1,i+1,2))*0.3333333
          go to 300
  290     continue
          write(14)nflag,(jflag(1,j),jflag(2,j),j=1,nflag)
          go to 390
  300     nflag=nflag+1
          i=i+1
          jflag(1,nflag)=i
          go to 220
c   extend row & write to disk.
  350     do 360 i=1,nap1
  360     a1(1,i,1)=a1(1,nap1,1)
          do 370 i=n2a,n2
  370     a1(1,i,1)=a1(1,n2a,1)
  380     iwr=iwr+1
          call swrda(iunita,iwr,a1,n22)
          if(iwr.eq.nxap1) go to 40
          go to 380
  390     do 400 i=1,nap1
  400     a1(1,i,2)=a1(1,nap1,2)
          do 410 i=n2a,n2
  410     a1(1,i,2)=a1(1,n2a,2)
          iwr=iwr+1
          call swrda(iunita,iwr,a1(1,1,2),n22)
          if(iwr.ge.n1a) go to 430
          do 420 i=nap1,n2a
          a1(1,i,1)=a1(1,i,2)
  420     a1(1,i,2)=a1(1,i,3)
          go to 40
  430     if(iwr.ge.n1) go to 450
  440     iwr=iwr+1
          call swrda(iunita,iwr,a1(1,1,2),n22)
          if(iwr.lt.n1) go to 440
  450 close(11)
      close(14)
  460     format(i5)
c  write headings
  470     if(ncase.eq.0) go to 590
      nlevel=nlevel+1
      if(nlevel.ne.inocal) go to 520
472   w1=0.
      w2=0.
       write(6,1115)nlevel,w1,w2
      ncase=1
      go to 470
520   z=barlev-zlev(nlevel)
      if(ifac.lt.0.or.z.lt.0.) go to 521
      wavlen=z*1.25+fac
      if(wavlen.lt.flim)wavlen=flim
      wdel=0.5*del+0.25*z
      w1=wavlen-wdel
      w2=wavlen+wdel
      f4=1./w1
      f3=1./w2
      btan34=1./(f4-f3)
521   no=nlevel+9
      write(crfile,560)no
560    format('mchess',i2,'.tmp')
      coname=crfile
590       pid=pi/180.
          rn1=1./(float(n1)*dx)
          rn2=1./(float(n2)*dy)
          c2=pi2*z
c  obtain complex fourier coefficients (f.f.t. a1).
          if(ncase.eq.0) go to 980
          jr=0
          do 970 i=1,n1
          jr=jr+1
          call srdda(iunitc,jr,a1,n22)
  970     call swrda(iunita,jr,a1,n22)
          go to 1000
  980     continue
  601     print*,'Compute .cof file (1), or supply your own (2):'
          print 600
  600     format(' *',$)
          read*,ifsup
          isign=-1
          open(12,access='direct',status='unknown',form='unformatted',
     1    file='slave2.tmp',recl=n22*4)
          open(15,access='direct',status='unknown',form='unformatted',
     1    file='fftfil.cof',recl=n22*4)
          if (ifsup.eq.2)go to 602
          call sfftmg(iunita,iunitb,n2,n1,a1,nri,isign,work)
          jr=0
          do 990 i=1,n1
          jr=jr+1
          call srdda(iunita,jr,a1,n22)
  990     call swrda(iunitc,jr,a1,n22)
  602     continue
c  a1 array now holds complex fourier coefficients.
      if(zlev(1).eq.barlev) go to 472
 1000     n1p1=n1+1
          n2p1=n2+1
          nnh1=(n1/2)+1
          nnh2=(n2/2)+1
c  operating on a1 fourier coefficients.
          jr=0
          jmr=n1+1
 1010     jr=jr+1
          if(jr.gt.nnh1) go to 1110
          jj=jr-1
          x=float(jj)*rn1
          xsq=x*x
          call srdda(iunita,jr,a1,n22)
          do 1090 i=1,n2
          ii=i-1
          if(i.gt.nnh2)ii=-(n2p1-i)
          y=float(ii)*rn2
          xysq=xsq+y*y
 1030     call scont(c2,xysq,a1(1,i,1),a1(2,i,1),f3,f4,btan34)
1090  continue
          call swrda(iunita,jr,a1,n22)
          if(jr.eq.1) go to 1010
          jmr=jmr-1
          if(jr.eq.jmr) go to 1110
          a1(2,1,1)=-a1(2,1,1)
          iwr=n2p1
          do 1100 i=2,nnh2
          iwr=iwr-1
          x=a1(2,i,1)
          y=a1(1,i,1)
          a1(2,i,1)=-a1(2,iwr,1)
          a1(1,i,1)=a1(1,iwr,1)
          a1(1,iwr,1)=y
          a1(2,iwr,1)=-x
 1100     continue
          call swrda(iunita,jmr,a1,n22)
          go to 1010
c  a1 array now holds complex coefficients operated on by filter
c  compute desired anomaly = inv. f.f.t. of a1
 1110     isign=1
      write(iw,1115)nlevel,w1,w2
 1115 format(' completed continuing to level ',i3,' with w1&w2=',2f10.2)
          call sfftmg(iunita,iunitb,n2,n1,a1,nri,isign,work)
          area=1./float(n1*n2)
c  a1 array now holds computed anomaly on datum = z.
c  create standard output file
 1120     open(11,status='unknown',form='unformatted',file=coname)
          write(kr)id,pgm,ny,nx,nz,yo,dy,xo,dx
          open(14,status='old',form='unformatted',
     1    file='flag.loc')
          dum=0.0
          iwr=0
          do 1210 jr=nxap1,n1a
          iwr=iwr+1
          call srdda(iunita,jr,a1,n22)
          read(14)nflag,(jflag(1,j),jflag(2,j),j=1,nflag)
          j=1
c   restore shape of original grid (i.e. insert flagged values).
          do 1160 i=nap1,n2a
          if(i.lt.jflag(1,j)) go to 1130
          if(i.eq.jflag(2,j)) go to 1140
          go to 1150
 1130     a1(1,i,1)=dval
          go to 1160
 1140     j=j+1
          if(j.le.nflag) go to 1150
          j=j-1
          jflag(1,j)=1028
 1150     a1(1,i,1)=a1(1,i,1)*area
 1160     continue
          go to 1190
 1170     do 1180 i=nap1,n2a
 1180     a1(1,i,1)=a1(1,i,1)*area
 1190     write(kr)dum,(a1(1,i,1),i=nap1,n2a)
 1210     continue
          close(11)
          ncase=1
      if(nlevel.eq.nlev) go to 1250
      close(14)
          go to 470
1250       close(14,status='delete')
      close(15,status='delete')
        close(10,status='delete')
          close(12,status='delete')
          go to 1290
 1270     write(6,1280)
 1280     format(' #entire first row is flagged values-program aborted')
          close(11)
          close(10)
          close(14)
 1290     return
          end
c***********************************************************************
       subroutine scont(c2,xysq,ar,ai,f3,f4,btan34)
c  subroutine 'cont' upward or downward continues input fourier
c  coefficients.
c***********************************************************************
      f=1.
      rad=sqrt(xysq)
      e=exp(c2*rad)
      if(f4.eq.0..or.c2.lt.0.) go to 50
      if(rad.gt.f4) go to 10
      go to 20
10    ar=0.
      ai=0.
      go to 60
20    if(rad.le.f3) go to 50
      f=(f4-rad)*btan34
50    ar=ar*e*f
      ai=ai*e*f
60    return
      end
c***********************************************************************
      subroutine srdbin(no,dat[huge],n)
c       subroutine 'srdbin' reads standard grids (binary).
c***********************************************************************
      dimension dat(n)
      read(no)dum,dat
      return
      end
c***********************************************************************
      subroutine swrbin(no,dat[huge],n)
c       subroutine 'swrbin' writes standard grids (binary).
c***********************************************************************
      dimension dat(n)
      data dum/0.0/
      write(no)dum,dat
      return
      end
c***********************************************************************
       subroutine srdda(no,ipos,dat[huge],n)
c       subroutine 'srrda' reads keyed sequential files.
c***********************************************************************
      dimension dat(n)
      read(no,rec=ipos)dat
      return
      end
c***********************************************************************
      subroutine swrda(no,ipos,dat[huge],n)
c       subroutine 'swrda' writes keyed sequential files.
c***********************************************************************
      dimension dat(n)
      write(no,rec=ipos)dat
      return
      end
