Convolution operator

# conv1dc
# convolution 1 dimensional circular boundary condition
# u and w are conjugate to each other
# axis: axis over which to convolve
# stagger:  0 normal conv   
#           1 shift by one grid point
#%
subroutine conv1dc(conj,u,t,k,w,sidx,  axis, scale,bnd)
implicit none
#include "commons"
integer            conj,  t,k,  sidx 
integer                                axis
real		  	 	             scale
WFIELD(u)
WFIELDLAY(u)
TENSOR(w)
TENSORLAY(w)
BOUND(bnd)
BOUNDLAY(bnd)
integer ishift

if (conj==0) {				#  the gradient
do ishift = nlo-stagger, nhi-stagger {
   w(sidx,:,:,:) = w(sidx,:,:,:) + 
        coef(ishift+stagger)*cshift(u(t,k,:,:,:),dim=axis,shift=ishift) * scale
   }
}
else{
do ishift = nlo-stagger, nhi-stagger {  #  the divergence
   u(t,k,:,:,:) = u(t,k,:,:,:) + 
       coef(ishift+stagger)*cshift(w(sidx,:,:,:),dim=axis,shift=ishift) * scale
   }
}
return
end

The finite difference operator is implemented as a 8 point convolution. The boundary conditions are circular, but a two line change will allow absorbing one-way boundary conditions. Even a spatially changing convolution operator (variable grid spacing) is mere one line change.