module array_ncnhconest_mod  

! non-causal non-stationary convolution, for an array operator. Adjoint is the filter.

  use ncnhconest1
  use ncnhconest2

  use row_mod

  implicit none

    integer, private                  :: ncoef,ncoef1
    real,private                      :: mu
    real,private                      :: eps1,eps2,eps3
    integer, private                  :: ndat

  contains

  subroutine array_ncnhconest_init (x_in,y_in,aa_in,bb_in,nhmax1_in,nhmax2_in, &
          & ncoeffs1_in,ncoeffs2_in,mu_in,eps1_in,eps2_in,eps3_in)

    real, dimension(:), pointer        :: x_in,y_in
    type(nfilter)                      :: aa_in,bb_in 
    integer                            :: nhmax1_in,nhmax2_in
    integer                            :: ncoeffs1_in,ncoeffs2_in
    real                               :: mu_in
    real                               :: eps1_in,eps2_in,eps3_in

    !! set auxiliary global variables
    ncoef1=ncoeffs1_in
    ncoef=ncoeffs1_in+ncoeffs2_in
    mu=mu_in
    eps1=eps1_in
    eps2=eps2_in
    eps3=eps3_in
    ndat=size(x_in)

    !! initialize the convolutional matrices for multiples and primaries
    call ncnhconest1_init(x_in,aa_in,nhmax1_in)
    call ncnhconest2_init(y_in,bb_in,nhmax2_in)

    !! initialize the row2 operator that computes M+muP by calling
    !! the operators ncnhconest1_lop and ncnhconest2_lop
    call row2_init(ncoef1,ncoef,mu) 

  end subroutine

  function array_ncnhconest_lop ( adj, add, a, y) result(stat)
    integer            :: stat 
    logical,intent(in) :: adj,add 
    real,dimension(:)  :: a,y 


    call adjnull1 (adj,add,a,y )

    !! call the actual array operator
    call row3(nrow2_lop,ncnhconest1_lop,ncnhconest2_lop,adj,add,a,y, &
         & eps1,eps2,eps3,ndat,ndat,ndat)

    stat=1

  end function array_ncnhconest_lop

  subroutine array_ncnhconest_close()

  end subroutine array_ncnhconest_close

end module array_ncnhconest_mod
