program Pef3
  use sep
  use helix
  use createhelixmod
  use bound
  use misinput
  use pef2
  use compress
  use helicon
  use print
 
  implicit none
  integer, dimension (:), allocatable :: n, n0, a, center, gap
  integer                             :: dim, n123, i, niter, na,ndim,ip,np123
  real,    dimension (:), pointer     :: d1,d2,diff,kk
  integer :: i1,i2,i3
  type( filter)                       :: a1,a2
  real :: ot,dt
  integer :: nt
  character(len=128) :: label
  logical :: doit

  call sep_init (SOURCE)
  dim = sep_dimension ()
  ndim=dim/2
  allocate (n (dim), n0 (ndim), center (ndim), a (ndim), gap (ndim))
  call from_history (n)
  call to_history ("n",n)
  call to_history("n",n,"diff")
  call from_par ("a", a)
  center = 1
  where (a (2:) > 1) center (:ndim-1) = a(:ndim-1)/2 + 1
  call from_par ("center", center, center)
  call from_par ("na", na, 0)
     gap = 0
     call from_par ("gap", gap, gap)
     a1 = createhelix( n(1:ndim), center, gap, a) ! allocate PEF
     a2 = createhelix( n(1:ndim), center, gap, a) ! allocate PEF
     n0 = n

  n123 = product (n(1:ndim))
  np123 = product (n(ndim+1:))
  allocate (d1 (n123), kk (n123)) 
  allocate(d2(n123),diff(n123))

  kk = 1.

  do i=1,dim
   call sep_get_data_axis_par("in",i,nt,ot,dt,label)
   call sep_put_data_axis_par("diff",i,nt,ot,dt,label)
  end do

  call boundn( n(1:ndim), n(1:ndim), a, a1)
  call boundn( n(1:ndim), n(1:ndim), a, a2)

!  call find_mask( (kk /= 0.), a1)   ! missing data
  call from_par( "niter", niter, size(a1%lag)*2 )

  

  do ip=1,np123
   if(mod(ip,20)==0) write(0,*) "working on",ip,np123
    call sreed("in",d1,size(d1)*4)
    call sreed("full",d2,size(d2)*4)
    diff=d1-d2
    doit=.true.
    if(sum(abs(diff)) <.0001) doit=.false.
    write(0,*) "CHECK",sum(d2),sum(d2),ip,doit
    if(doit) then
      call find_pef ( d2, a1, niter,i1,i2,i3)         ! estimate aa
!    call find_pef ( d1, a2, niter,i1,i2,i3)         ! estimate aa
      call helicon_init(a1)
      i=helicon_lop(.false.,.false.,d1,diff)
      call helicon_close()
      where(a1%mis) diff=0
    else 
      diff=0
    end if
    call srite("out",diff,size(d2)*4)
    if(doit) then
      call helicon_init(a1)
      i=helicon_lop(.false.,.false.,d2,diff)
      call helicon_close()
      where(a1%mis) diff=0
    else
      diff=0
    end if
    call srite("diff",diff,size(d2)*4)
  end do

  call sep_close ()
  deallocate (d1,d2,diff, kk);
  deallocate (n, n0, center, a, gap)
  call deallocatehelix( a1)
end program 
