*--------------------------------------------------------------
* FFILT - convolution of Butterworth low-pass filter with
* real data in U(M*N). The cutoff wavenumber K0 is in
* grid units in frequency space.
*---------------------------------------------------------------
subroutine ffilt(u,m,n,k0,k1,work,nwork)
real*4 u(2*(m/2+1)*n),k0,k1
complex*8 work(nwork)
integer*4 nu(2)
data ndim/2/
nu(1)=m
nu(2)=n
* do the forward transform of a real 2d array
iform=0
isign=1
call fourt(u,nu,ndim,isign,iform,work,nwork)
m2=m/2+1
* filter the transformed array
call hc_filt(u,m2,n,k0,k1)
* inverse transform of a half complex array and renormalise
iform=-1
isign=-1
fn=float(m*n)
call fourt(u,nu,ndim,isign,iform,work,nwork)
do i=1,m*n
u(i)=u(i)/fn
end do
return
end
real*4 function hypot(x,y)
real*4 x,y
hypot=sqrt(x*x+y*y)
return
end
*---------------------------------------------------------------
* HC_FILT - apply real filter function to Fourier transformed
* array which is half-complex, i.e. derived from a real
* signal. The transformed array is in U and D0 is the
* filter cut-off wavenumber. D1 is the cutoff in y grid
* units and KV/KU scales with D1/D0 for non-square grids.
* IP is the order of the filter.
*---------------------------------------------------------------
subroutine hc_filt(u,m,n,d0,d1)
complex*8 u(m*n),t
real*4 d0,d1,ku,kv,kv2
data ip/4/
j2=n/2+1
kv2=-n/2
do j=1,n
j0=(j-1)*m
if(j.le.j2)then
kv=j-1
else
kv=kv2+j-j2
endif
kv=kv*(d1/d0)
do i=1,m
loc=j0+i
ku=i-1
d=hypot(ku,kv)
t=cmplx(blpf(d,d0,ip),0.)
u(loc)=u(loc)*t
end do
end do
return
end
*---------------------------------------------------------------
* BLPF - Butterworth low-pass filter with cut-off D0
* in frequency space. Filter is order IP.
*---------------------------------------------------------------
real*4 function blpf(d,d0,ip)
real*4 d,d0
integer*4 ip
blpf=1./(1. + 0.414*(d/d0)**(2*ip))
return
end