module nhconest {           # Nonstationary convolution, adjoint is the filter.
use nhelix
  real, dimension(:), pointer                 :: x
  type( nfilter)                              :: aa
  integer                                     :: nhmax
  integer, private                            :: np  
#% _init( x, aa, nhmax)
    np = size( aa%hlx)
#% _lop( a( nhmax, np), y(:))
    integer                        :: ia, ix, iy, ip
    integer, dimension(:), pointer :: lag
    do iy = 1, size( y) {  if( aa%mis( iy)) cycle
        ip = aa%pch( iy); lag => aa%hlx( ip)%lag
    	do ia = 1, size( lag) {
		ix = iy - lag( ia);   if( ix < 1) cycle
		     if( adj)    a( ia, ip) +=  y( iy) * x( ix)
		     else        y( iy) +=  a( ia, ip) * x( ix)
		 }
        }
}