program nfilt_it
use sep 
use sep3d_struct_mod
use npolydiv
use nhelicon
use io_nfilt
implicit none

logical             :: adj,add,prec
real,allocatable    :: model(:),data(:),tmp(:)
type(sep3d)         :: mod,dat
character(len=128)  :: label
type(nfilter)       :: nfilt,nfilt2
integer             :: n,ierr

call sep_init(SOURCE)
call init_sep3d("data",dat,"SCRATCH")
call init_sep3d("model",mod,"SCRATCH")
call auxclose("data")
call auxclose("model")
call auxinout("model")
call auxinout("data")
call auxin("filter")
 write(0,*) "XX "

n=product(mod%n(1:sep3d_ndims(mod)))
if (n/=product(dat%n(1:sep3d_ndims(dat)))) call seperr("expecting model and data to be the same size")
if(.not. read_nfilt("filter",nfilt)) call seperr("trouble reading in filter")
if(.not. read_nfilt("filter2",nfilt2)) call seperr("trouble reading in filter")
if(size(nfilt%pch)/=n) call seperr("filter and data size not equiv")

call from_param("add",add,.false.)
call from_param("adj",adj,.false.)
call from_param("prec",prec,.false.)

allocate(model(n),data(n),tmp(n))
model=0
data=0
!call sreed("data",data,size(data)*4)
call sreed("model",model,size(model)*4)

!if(.not. prec) then
 call nhelicon_init(nfilt)
!else
  call npolydiv_init(n,nfilt2)
  ierr=npolydiv_lop(.false.,add,model,tmp)
  call srite("data",tmp,size(data)*4)
  ierr=nhelicon_lop(.false.,add,tmp,data)
!end if
  call srite("data",data,size(data)*4)
end program
