c***********************************************************************
c   subroutine 'sfftfil' removes flagged values from grid, transforms
c   data, calls routines to apply filter, and writes filtered grid
c***********************************************************************
          subroutine sfftfil(a1,work[huge],nx,nri,itype,nxa,n1,n2,id2)
          dimension a1(2,n2,nri),work(id2),jflag(2,20),title(14)
          character flist(10)*6,help*4,quer*1,fopr*6
          character*50 coname,fname,fileo
          character*2 drive
          common/mane/id(14),pgm(2),nz,yo,xo,dx,dy,iw,kr,ny,nadd,ir
          common/basic/c1,el,em,en,bl,bm,bn
          common/parm1/fname,coname,iopt1,iopt2,idval,icoef,ddx,ddy
          common/parm2/den,dec,xinc,bmag,bdec,binc,z,thet1,thet2,istr,w1
     1    ,w2,w3,w4,drive
          common/tr/prad,p,q,th1p,th2p
          data help/'help'/,nlist,flist/10,'psdgrv','redpol','upcont',
     1   'dncont','2ndver','strike','banpas','psdmag','1stver','nofilt'/
          data iunita/10/,iunitb/12/,dval/1.0e38/,ddval/1.0e37/
          data gamma/6.67323/,pi2/6.28318531/,pi/3.141592654/,iunitc/15/
          ncase=0
      n22=2*n2
          if(icoef.ne.-1) go to 10
          ncase=1
         open(15,access='direct',status='old',form='unformatted',
     1   file='fftfil.cof',recl=n22*4)
         open(12,access='direct',status='unknown',form='unformatted',
     1   file=drive//'slave2.tmp',recl=n22*4)
          close(11)
          go to 20
c   remove mean using boundary values.
   10     if(iopt2.lt.0) go to 20
          call srmean(nx,ny,a1)
          fileo='fftfil.grd'
          open(11,status='old',form='unformatted',file=fileo)
          read(kr)title,pgm,ny,nx,nz,yo,dy,xo,dx
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=drive//'slave1.tmp',recl=n22*4)
      if(ncase.eq.1) go to 590
          if(idval.eq.1)open(14,status='unknown',form='unformatted',
     1    file='flag.loc')
          irr=1
          iwr=0
          jr=0
   40     jr=jr+1
          if(idval.lt.1) go to 310
c   idval>0: removal of flagged values.
          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).lt.ddval) go to 90
          do 60 i=nap1,n2a
   60     if(a1(1,i,1).lt.ddval) 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).lt.ddval) 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).lt.ddval) 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).lt.ddval) go to 220
          do 180 i=nap1,n2a
  180     if(a1(1,i,2).lt.ddval) go to 190
          i=nap1
          jflag(1,1)=n2a+1
          if(a1(1,nap1,3).ge.ddval)a1(1,nap1,2)=a1(1,nap1,1)
          if(a1(1,nap1,3).lt.ddval)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.ddval) 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).lt.ddval) 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).lt.ddval) go to 270
  250     if(a1(1,i,3).ge.ddval) 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.ddval) 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   idval.eq.0: read & check input data.
  310     if(jr.ne.1)irr=2
          read(kr)dum,(a1(1,i,irr),i=nap1,n2a)
          do 320 ii=nap1,n2a
  320     if(a1(1,ii,irr).ge.ddval) go to 330
          if(jr.ne.1) go to 390
          go to 350
  330     write(iw,340)jr,ii
  340     format(' #flagged value row=',i4,3x,'col=',i4,' may be more?')
          close(11)
          close(10)
          return
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
          if(idval.lt.1) go to 40
          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 if(iopt2.ne.0)close(11)
      if(iopt2.eq.0)close(11,status='delete')
      if(idval.eq.1)close(14)
  460     format(i5)
c  write headings
  470     if(ncase.eq.0) go to 590
  480 write(iw,485)
  485 format(' enter new operator: ',$)
          read(ir,490)fopr
  490     format(a6)
  500     if(fopr.eq.help) go to 530
          do 510 i=1,nlist
          if(fopr.ne.flist(i)) go to 510
          itype=i
          go to 560
  510     continue
          write(iw,520)
  520     format(' #invalid operator')
  530     write(iw,540)
  540     format(' quick list of options')
          write(iw,550)flist
  550     format(2x,a6)
          go to 480
  560 write(iw,565)
  565 format(' enter new output file name: ',$)
          read(ir,570)coname
  570     format(a)
      write(iw,575)
575   format(' enter new title: ',$)
          read(ir,580)id
  580     format(14a4)
      write(iw,585)
  585 format(' parm change? if y enter parameters via &parms): ',$)
          read(ir,570)quer
          if(quer.ne.'y') go to 590
          call snamel(ir)
  590     ip=iopt1
          ih=iw
          if(ip.ne.0)ih=ip
          pid=pi/180.
          rn1=1./(float(n1)*dx)
          rn2=1./(float(n2)*dy)
          c2=pi2*z
          go to (620,640,670,690,710,730,780,600,890,595,930,910),itype
  595	  write(ih,597)
  597     format(1h1,'transformation of vertical-component magnetic ',
     1    'data to total-field data,',/,' by means of the fast fourier',
     2    ' transform and poissons equation.',/)
          go to 640
  600     write(ih,610)
  610     format(1h1,'transformation of the gravity field to the ',
     1    'pseudo-magnetic field,',/,' by means of the fast fourier',
     2    ' transform and poissons equation.',/)
          go to 640
  620     write(ih,630)
  630     format(1h1,'transformation of total magnetic intensity field',
     1    ' to pseudo-gravitational field,',/,' by means of the fast ',
     2    'fourier transform and the poisson equation.',/)
  640     rbdec=bdec*pid
          rbinc=binc*pid
          cosrb=cos(rbinc)
          bl=cosrb*cos(rbdec)
          bm=cosrb*sin(rbdec)
          bn=sin(rbinc)
          rdec=dec*pid
          rinc=xinc*pid
          cosr=cos(rinc)
          el=cosr*cos(rdec)
          em=cosr*sin(rdec)
          en=sin(rinc)
          if(itype.eq.2) go to 650
          c1=(gamma*den)/(pi2*bmag)
          go to 950
  650     write(ih,660)
  660     format(1h1,'transformation of the total magnetic field ',
     1    'intensity to the pole,',/,' by means of the fast fourier ',
     2    'transform and the poisson equation.',/)
          go to 950
  670     write(ih,680)
  680     format(1h1,'upward continuation of field values,',/,
     1    ' by means of the fast fourier transform.',/)
          go to 950
  690     write(ih,700)
  700     format(1h1,'downward continuation of field values,',/,
     1    ' by means of the fast fourier transform.',/)
          go to 950
  710     write(ih,720)
  720     format(1h1,'second vertical derivative operation,',/,
     1    ' by means of the fast fourier transform.',/)
          c4=4.*pi*pi
          go to 950
  730     write(ih,740)
  740     format(1h1,'strike sensitive filtering using the fast ',
     1    'fourier transform.',/)
          p=1.
          q=0.0
          if(istr.ne.1) go to 760
          if(thet1.lt.0.0.and.thet2.ge.0.0) go to 750
          go to 770
  750     p=0.0
          q=1.
          go to 770
  760     if(thet1.lt.0.0.and.thet2.ge.0.0) go to 770
          p=0.0
          q=1.
  770     th1p=180.-thet1
          th2p=180.-thet2
          if(th1p.gt.180.)th1p=th1p-180.
          if(th2p.gt.180.)th2p=th2p-180.
          prad=180./pi
          go to 950
  780     write(ih,790)
  790     format(1h1,'band-pass filter applied to data',/)
          if(w1.ne.0.0) go to 800
          f4=1.0e+30
          go to 810
  800     f4=1./w1
  810     if(w2.ne.0.0) go to 820
          f3=1.0e+30
          go to 830
  820     f3=1./w2
  830     if(w3.ne.1.0e+30) go to 840
          f2=0.
          go to 850
  840     f2=1./w3
  850     if(w4.ne.1.0e+30) go to 860
          f1=0.
          go to 870
  860     f1=1./w4
  870     if(w3.eq.w4) go to 880
          btan12=1./(f2-f1)
  880     if(w1.eq.w2) go to 950
          btan34=1./(f4-f3)
          go to 950
  890     write(ih,900)
  900     format(1h1,'first vertical derivative operation,',/,
     1    ' by means of the fast fourier transform.',/)
          c4=pi2
          go to 950
  910     write(ih,920)
  920     format(1h1,'no filter operation performed: fourier coeff.''s',
     1    ' are input and inverse transformed',/)
          go to 950
  930     write(ih,940)
  940     format(1h1,'no filter operation performed: only fourier ',
     1    'coef.s are output in fftfil.cof',/)
  950     write(ih,960)id,z,den,bmag,xinc,dec,binc,bdec,thet1,thet2,
     1     istr,w1,w2,w3,w4,drive(1:1),nadd,dx,dy,xo,yo,ny,nx,n2,n1,
     2     fname,coname
  960     format(' title:',6x,14a4,/' parameters: datum level (z)=',f7.3
     1,11x,'density (den)=',2x,f7.3/13x,'intensity of magnetization',15x
     1,'(bmag)=',f9.3/13x,'incl & decl of geomagnetic field     (xinc de
     1c)=',2f9.3/13x,'incl & decl of magnetization vector (binc bdec)=',
     12f9.3/13x,'azimuthal filter angles',11x,'(thet1 thet2)=',
     12f9.3/13x,'azimuthal filter option (+1=pass,-1=cut) (istr)=',
     1i5/13x,'bandpass wavelengths',14x,'(w1,w2,w3&w4)=',/35x,4g11.4/,
     113x,'(RAM) drive letter for temporary files  (drive)= ',2x,a2/13x,
     1'requested number of additional rows/cols (nadd)=',i5/,
     1' grid specs: grid interval (dx dy)= ',2g14.6/13x,
     1'origin of grid (xo yo)=',2g14.6/13x,'no. of columns & rows=',2x,
     12i5/13x,'extended columns & rows=',2i5/13x,
     8'input file name:  ',a49/13x,'output file name: ',a49,/)
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)
          if(itype.eq.12) go to 1110
          go to 1000
  980     isign=-1
          open(12,access='direct',status='unknown',form='unformatted',
     1    file=drive//'slave2.tmp',recl=n22*4)
          call sfftmg(iunita,iunitb,n2,n1,a1,nri,isign,work)
          open(15,access='direct',status='unknown',form='unformatted',
     1    file='fftfil.cof',recl=n22*4)
          jr=0
          do 990 i=1,n1
          jr=jr+1
          call srdda(iunita,jr,a1,n22)
  990     call swrda(iunitc,jr,a1,n22)
          if(itype.eq.11) go to 1250
      write(iw,995)
  995 format(' forward transform completed')
c  a1 array now holds complex fourier coefficients.
 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
          go to (1020,1020,1030,1030,1040,1050,1060,1020,1070,1020),
     1  itype
 1020     if(xysq.eq.0.0) go to 1080
          call spspol(itype,x,y,xysq,a1(1,i,1),a1(2,i,1))
          go to 1090
 1030     call scont(c2,xysq,a1(1,i,1),a1(2,i,1))
          go to 1090
 1040     call ssvert(c4,xysq,a1(1,i,1),a1(2,i,1))
          go to 1090
 1050     if(xysq.eq.0.0) go to 1090
          call strend(x,y,a1(1,i,1),a1(2,i,1))
          go to 1090
 1060     call sbpass(a1(1,i,1),a1(2,i,1),f1,f2,f3,f4,btan12,btan34,
     1    xysq)
          go to 1090
 1070     call sfvert(c4,xysq,a1(1,i,1),a1(2,i,1))
          go to 1090
 1080     a1(1,1,1)=0.0
          a1(2,1,1)=0.0
 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)
 1115 format(' completed applying filter to coefficients')
          call sfftmg(iunita,iunitb,n2,n1,a1,nri,isign,work)
      write(iw,1118)
 1118 format(' inverse transform completed')
          area=1./float(n1*n2)
c  a1 array now holds computed anomaly on datum = z.
c  create standard output file and list output.
          if(iopt1.eq.0) go to 1120
          write(ip,1240)id,z
          f16=float(ny)/16.
          i16=ny/16
          if(abs(f16-float(i16)).gt.0.0001)i16=i16+1
 1120     open(11,status='unknown',form='unformatted',file=coname)
          write(kr)id,pgm,ny,nx,nz,yo,dy,xo,dx
          if(idval.eq.1)open(14,status='old',form='unformatted',
     1    file='flag.loc')
          if(idval.lt.0)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)
          if(idval.eq.0) go to 1170
          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)=2048
 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)
           if(iopt1.eq.0) go to 1210
          write(ip,1220)iwr
          ih=nap1-1
          do 1200 ii=1,i16
          iout=ih+1
          ih=iout+15
          if(ih.gt.n2a)ih=n2a
 1200     write(ip,1230)(a1(1,i,1),i=iout,ih)
 1210     continue
 1220     format(/,1x,'row ',i3)
 1230     format(3x,16f7.1)
 1240     format(1h1,14a4,//' calculated anomaly on level z=',f7.3,//)
          close(11)
      write(iw,1245)
 1245 format(' additional filtering to original grid?(y or n) ',$)
          read(ir,570)quer
          if(quer.ne.'y') go to 1250
      if(idval.ne.0)close(14)
          ncase=1
          go to 470
 1250     if(idval.ne.1) go to 1255
          if(icoef.eq.0)close(14,status='delete')
          if(icoef.ne.0)close(14)
      go to 1260
 1255 if(idval.eq.0) go to 1260
      close(14)
 1260     close(10,status='delete')
          close(12,status='delete')
      if(icoef.eq.0)close(15,status='delete')
      if(icoef.ne.0)close(15)
          go to 1290
 1270     write(6,1280)
 1280     format(' #entire first row is flagged values-program aborted')
          close(11)
          close(10)
          if(idval.eq.1)close(14)
 1290     return
          end

c***********************************************************************
       subroutine scont(c2,xysq,ar,ai)
c  subroutine 'cont' upward or downward continues input fourier
c  coefficients.
c***********************************************************************
      xy=sqrt(xysq)
      e=exp(c2*xy)
      ar=ar*e
      ai=ai*e
      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




