subroutine filt(savex,savek,work,nx,ny,xlong,ylong,iwant,par)
c
	implicit none
c
c All arguments are sent; savex is altered by the program.
c
	integer nx, ny, iwant(6)
	real*4 savex(nx,ny), xlong, ylong, par(13)
	complex savek(nx/2+1, ny), work(nx)
c
c savex,savek are declared and equivalenced in the calling program
c
c xlong, ylong are the lengths of the array in km
c
c iwant selects which one or more filters to do, and whether they
c are low-pass or high-pass filters.  Elements of iwant correspond
c to filters like this:
c	iwant(1) -- Cosine taper between two harmonic degrees
c	iwant(2) -- Gaussian filter with a sigma value
c	iwant(3) -- Wiener filter W2(k) as in Smith & Sandwell,
c		Nov 1994 JGR red.
c	iwant(4) -- Downward (+) or upward (-) continuation.
c	iwant(5) -- Two-layer flexural bathymetric prediction.
c	iwant(6) -- Kill wavelengths > ylong.
c If an element of iwant is non-zero, the corresponding filter will
c be applied to the data.  For the Cosine and Gaussian filters,
c the sign of iwant indicates lowpass (-) or highpass (+); the sign
c is also used to indicate up/down continuation.  The sign is ignored
c for Wiener filter (which is always lowpass) and for prediction.
c
c par is an array of parameters used by the various possible filters.
c values stored in par will be used as needed according to iwant.
c elements of par are interpreted like this:
c	par(1) -- harmonic degree h1 (for cosine)
c	par(2) -- harmonic degree h2 (for cosine; h2 > h1)
c	par(3) -- gaussian sigma (in km)
c	par(4) -- Wiener constant A (in km**4; 9500 in Paper-1)
c	par(5) -- z, depth in km for up/down continuation
c	par(6) -- lambda, wavelength used for flexure in km
c	par(7) -- thickness of layer 2, km
c	par(8) -- thickness of layer 3, km
c	par(9) -- density of seawater, SI
c	par(10) - density of layer 2, SI (also load density)
c	par(11) - density of layer 3, SI
c	par(12) - density of mantle, SI
c	par(13) - harmonic degree for kill
c
	integer n(2), ntot, i, j, k, nxhalf, nxh1, nyhalf, nyh1
	real*8 pi, rntot, rxl, ryl, cosk1, cosk2, cosfac, gausfac, gama
	real*8 wiener, zdown, xk, yk, rk, gamfac, z32, zm3, phi
	real*8 admit, rho32, rhom3, lambda, yksq, ff, ftot, wiendown
	real*8 kkill
	logical	nobody

c
c Make sure something needs to be done:
c
	k=1
	nobody = .true.
	do 10 while (nobody .and. k .le. 6)
	  nobody = (iwant(k) .eq. 0)
	  k = k+1
  10	continue
	if (nobody) then
	  write(*,901)
 901	  format('filt:  NON-FATAL ERROR.  iwant is null.')
	  return
	end if
c
c Set up some constants, etc.
c
	pi = 3.14159265358979323846
        gama = 6.67d-11
	n(1) = nx
	n(2) = ny
	ntot = nx * ny
	nxhalf = nx/2
	nyhalf = ny/2
	nxh1 = nxhalf + 1
	nyh1 = nyhalf + 1
	rntot = 1.0d+00/(ntot)
	rxl = 1.0d+00/xlong
	ryl = 1.0d+00/ylong
c
	if (iwant(1) .ne. 0) then
	  if (par(2) .le. par(1) ) then
	    write(*,902)
 902	    format('filt:  ERROR in cosine taper.  par(2) <= par(1).')
	    return
	  end if
	  cosk1 = par(1)/4.0d+04
	  cosk2 = par(2)/4.0d+04
	  cosfac = pi /(cosk2-cosk1)
	  if (iwant(1) .lt. 0) then
	    write(*,903) 1.0/cosk1, 1.0/cosk2
 903	    format('filt:  Low-pass cosine from ',f7.1,' to ', f7.1)
	  else
	    write(*,904) 1.0/cosk1, 1.0/cosk2
 904	    format('filt:  High-pass cosine from ',f7.1,' to ', f7.1)
	  end if
	end if
c
	if (iwant(2) .ne. 0) then
	  if (par(3) .le. 0.0) then
	    write(*,905)
 905	    format('filt:  ERROR.  par(3) (sigma) <= 0.0')
	    return
	  end if
	  gausfac = pi * par(3) * dsqrt(2.0d+00)
	  if (iwant(2) .lt. 0) then
	    write(*,906) par(3)
 906	    format('filt:  Low-pass gaussian, sigma = ',f7.2)
	  else
	    write(*,907) par(3)
 907	    format('filt:  High-pass gaussian, sigma =  ',f7.2)
	  end if
	end if
c
	if (iwant(3) .ne. 0) then
	  if (par(4) .le. 0.0) then
	    write(*,908)
 908	    format('filt:  ERROR.  par(4) (Wiener constant) <= 0.0')
	    return
	  end if
	  wiener = par(4)
	  write(*,909) par(4)
 909	  format('filt:  Wiener filter with constant = ',f7.2)
	end if
c
	if (iwant(4) .ne. 0) then
	  if (par(5) .lt. 0.0) then
	    write(*,910)
 910	    format('filt:  ERROR in up/down.  par(5) (depth) <= 0.')
	    return
	  end if
	  if (iwant(4) .lt. 0) then
	    zdown = -pi * 2.0d+00 * par(5)
	    write(*,911) par(5)
 911	    format('filt:  Upward continuation, z = ',f4.2)
	  else
	    zdown = pi * 2.0d+00 * par(5)
	    write(*,912) par(5)
 912	    format('filt:  Downward continuation, z = ', f4.2)
	  end if
	end if
c
	wiendown = 2.0d+00 * dabs(zdown)
c
	if (iwant(5) .ne. 0) then
	  do 914 k=6,12
	    if (par(k) .le. 0.0) then
	      write(*,913) k
 913	      format('filt:  ERROR. Flexure Parameter ', i2, ' <= 0.0')
	      return
	    end if
 914	  continue
	  lambda = par(6)
	  gamfac = 1.0d05 * 2.* pi * gama * (par(10)-par(9))
	  rho32 = (par(11)-par(10)) / (par(12)-par(10))
	  rhom3 = (par(12)-par(11)) / (par(12)-par(10))
	  z32 = -2.0d+00 * pi * par(7)
	  zm3 = -2.0d+00 * pi * (par(7) + par(8))
	  write(*,915)par(6),par(7),par(8),par(9),par(10),par(11),par(12)
 915	  format('filt:  Flex. predict. w. ',f7.2,2f4.1,4f7.1)
	end if
c
	if (iwant(6) .ne. 0) then
	  if (par(13) .le. 0.0) then
	    write(*,916)
 916	    format('filt:  ERROR.  Bad kill harmonic.')
	    return
	  end if
	  kkill = par(13)/4.0d+04
	end if
c
c
c Get here when things are ready to go.  Call forward transform:
c
	call fourt(savex,n,2,-1,0,work,nx)
c
c Loop over the rows and y wavenumbers:
c
	do 100 i=1, ny
	  if (i .le. nyh1) then
	    yk = (i-1) * ryl
	  else
	    yk = (ny-i+1) * ryl
	  end if
	  yksq = yk * yk
c
c Loop over the columns and x wavenumbers.
c -- Clue for speed later:  Maybe compute once and store in work()
c
	  do 100 j=1, nxh1
	    xk = (j-1)*rxl
c
	    rk = dsqrt(xk*xk + yksq)
c
c If kill long, can jump out here and save time
c
	    if (iwant(6) .ne. 0 .and. rk .lt. kkill) then
	      ftot = 0.0d+00
	      go to 90
	    end if
c
c Initialize ftot, the total filter factor, for the normalization.
c Then for each filter we want to do, multiply ftot by that filter:
c
	    ftot = rntot
c
	    if (iwant(1) .ne. 0) then
c
c Cosine.  Make lowass, then switch if needed
c
	      if (rk .le. cosk1) then
		ff = 1.0
	      else if (rk .ge. cosk2) then
	        ff = 0.0
	      else
	        ff = 0.5d+00 * (1.0d+00 + dcos(cosfac*(rk-cosk1)))
	      end if
	      if (iwant(1) .gt. 0) ff = 1.0d+00 - ff
	      ftot = ftot * ff
	    end if
c
	    if (iwant(2) .ne. 0) then
c
c Gaussian.  Make lowass, then switch if needed
c
	      ff = dexp(-(gausfac*rk)**2)
	      if (iwant(2) .gt. 0) ff = 1.0d+00 - ff
	      ftot = ftot * ff
	    end if
c
	    if (iwant(3) .ne. 0) then
c
c Wiener.
c
	      ff = 1.0d+00/(1.0d+00 + dexp(wiendown*rk)*wiener*rk**4)
	      ftot = ftot * ff
	    end if
c
	    if (iwant(4) .ne. 0) then
c
c Upward or downward continuation.  Sign set in zdown.
c
	      ff = dexp(rk * zdown)
	      ftot = ftot * ff
	    end if

c
	    if (iwant(5) .ne. 0) then
c
c Flexurally compensated bathymetric prediction.
c
	      phi = 1.0d+00 / (1.0d+00 + (lambda * rk)**4)
	      admit = gamfac*(1.0d+00 - phi*(rho32*dexp(z32*rk) +
     +		rhom3*dexp(zm3*rk) ) )
	      if (admit .eq. 0.0) then
	        ff = 0.0d+00
	        if (rk .gt. 0.0) then
	          write (*, 991)
 991	          format('filt:  ERROR.  Inverse admittance is singular.')
	        end if
	      else
	        ff = 1.0d+00 / admit
	      end if
	      ftot = ftot * ff
	    end if
c
c Get here when ftot has the complete multiplying and normalizing factor.
c
  90	    savek(j,i) = savek(j,i) * ftot
c
c End of both loops:
c
 100	continue
c
c Now call the inverse transform and return
c
	call fourt(savex,n,2,1,-1,work,nx)
c
	return
	end

C23456789012345678901234567890123456789012345678901234567890123456789012