module npef  
! Estimate non-stationary PEF
  use nhconest
  use npolydiv2
  use cgstep_mod
  use solver_mod
  implicit none
  contains
  subroutine find_pef( dd, aa, rr, niter, eps, nh)
    integer,          intent( in)    :: niter, nh      
    ! number of iterations
    real,             intent( in)    :: eps	       ! epsilon
    type( nfilter)                   :: aa             
    ! estimated PEF output.
    type( nfilter),   intent( in)    :: rr             
    ! roughening filter.
    real,    dimension(:), pointer   :: dd             ! input data
    integer                          :: ip, ih, np, nr 
    ! filter lengths
    real, dimension (:), allocatable :: flt	       
    ! np*na filter coefs
    np = size( aa%hlx)			! data   length
    nr = np*nh
    allocate( flt( nr))
    call nhconest_init( dd, aa, nh)
    call npolydiv2_init( nr, rr)
    call solver_prec( nhconest_lop, cgstep, x= flt, dat= -dd, niter=&
      & niter,prec= npolydiv2_lop, nprec= nr, eps= eps)
    call   cgstep_close()
    call npolydiv2_close()
    call nhconest_close()
    do ip = 1, np
      do ih = 1, size( aa%hlx(ip)%lag)
        aa%hlx( ip)%flt( ih) = flt( (ip-1)*nh + ih)
      end do
    end do 
    deallocate( flt)
  end subroutine 
end module 
