module cdoubint	 
! double integration
  use adj_mod
  use chain_mod
  implicit none
  real ,private,  dimension (:), allocatable  :: tmp 
  contains
  subroutine cdoubint_init ( nx )
    integer :: nx
    if (.not. allocated(tmp)) then
      allocate(tmp ( nx)) 
    end if
  end subroutine
  function cdoubint_lop ( adj, add, x, y) result(stat)
    integer            :: stat 
    logical,intent(in) :: adj,add 
    real,dimension(:)  :: x,y 
    call adjnull (adj,add,x,y )
    call cdoubint_lop2(adj,add,x,y )
    stat=0
  end function 
  subroutine cdoubint_lop2(adj,add,x,y)
    logical,intent(in) :: adj,add 
    real, dimension (:)  :: x 
    real, dimension (:)  :: y 
    if (adj) then
      call normal (causint_lop, .true., y, x, tmp)
    else
      call normal (causint_lop, .true., x, y, tmp)
    end if
  end subroutine 
  subroutine cdoubint_close()
    deallocate( tmp )
  end subroutine 

  function causint_lop ( adj, add, xx, yy) result(stat)
    integer            :: stat 
    logical,intent(in) :: adj,add 
    real,dimension(:)  :: xx,yy 
    call adjnull (adj,add,xx,yy )
    call causint_lop2(adj,add,xx,yy )
    stat=0
  end function 
  subroutine causint_lop2(adj,add,xx,yy)
    logical,intent(in) :: adj,add 
    real, dimension (:)  :: xx 
    real, dimension (:)  :: yy 
    integer i, n       
    real t
    n = size (xx)
    t = 0.
    if ( adj) then
      do i= n, 1, -1  
        t = t + yy(i)
        xx(i) = xx(i) + t
      end do
    else
      do i= 1, n  
        t   = t + xx(i)
        yy(i) = yy(i) + t
      end do
    end if
  end subroutine 
  subroutine causint_close()
  end subroutine 


end module 

