! Copyright 1999 Stanford University
!All rights reserved
!
!Author: Sergey Fomel and Robert Clapp
module npolydiv_patch  
! Helix polynomial division
  use adj_mod
  use nhelix
  implicit none
  integer, private :: nd 
  type( nfilter), private :: aa 
  real ,private,  dimension (:), allocatable  :: tt 
  integer ,private,  dimension(:), pointer  :: lag 
  real ,private,     dimension(:), pointer  :: flt 
  logical, private, dimension(:), pointer    :: local_bounds
  contains
  subroutine npolydiv_patch_init ( nd_in,aa_in,bound)
    integer    :: nd_in 
    type( nfilter)    :: aa_in 
    logical,optional,pointer :: bound(:)

    if(present(bound))  then
      local_bounds=>bound
    else 
      nullify(local_bounds)
    end if
    nd = nd_in 
    aa = aa_in 
    if (.not. allocated(tt)) then
      allocate(tt ( nd)) 
    end if
  end subroutine
  function npolydiv_patch_op ( adj, add, xx, yy,index) result(stat)
    integer            :: stat 
    logical,intent(in) :: adj,add 
    real,dimension(:)  :: xx,yy 
    integer :: index(:)
    call adjnull (adj,add,xx,yy )
    call npolydiv_patch_lop2(adj,xx,yy,index )
    stat=0
  end function 
  subroutine npolydiv_patch_lop2(adj,xx,yy,index)
    logical,intent(in) :: adj 
    real, dimension (:)  :: xx 
    real, dimension (:)  :: yy 
    integer :: index(:)
    integer                        :: ia, ix, iy, ip
    tt = 0.
    if ( adj) then
      tt = yy
      do iy= nd, 1, -1  
        if (associated(aa%mis)) then
          if (aa%mis(index(iy))) then
            tt(iy)=0.
            cycle
          end if
        end if
        if (associated(local_bounds)) then
          if (local_bounds(iy)) then
            tt(iy)=0.
            cycle
          end if
        end if
        ip = aa%pch(index( iy))
        lag => aa%hlx( ip)%lag
        flt => aa%hlx( ip)%flt
        do ia = 1, size( lag)
          ix = iy - lag( ia)
          if ( ix < 1) then
            cycle
          end if
          tt( ix) =			tt( ix) -  flt( ia) * tt( iy)
        end do
      end do 
      xx =	xx + tt
    else
      tt = xx
      do iy= 1, nd  
        if (associated(local_bounds)) then
          if (local_bounds(index(iy))) then
            tt(iy)=0.
            cycle
          end if
        end if
        if (associated(aa%mis)) then
          if (aa%mis(iy)) then
            tt(iy)=0.
            cycle
          end if
        end if
        ip = aa%pch( index(iy))
        lag => aa%hlx( ip)%lag
        flt => aa%hlx( ip)%flt
        do ia = 1, size( lag)
          ix = iy - lag( ia)
          if ( ix < 1) then
            cycle
          end if
          tt( iy) =			tt( iy) -  flt( ia) * tt( ix)
        end do
      end do 
      yy =	yy + tt
    end if 
  end subroutine 
  subroutine npolydiv_patch_close()
    deallocate( tt )
  end subroutine 
end module 
