module opera  
!
  use adj_mod

  implicit none
  integer, private				:: nw,nx,nz,nreal

  contains

!==================================================
  subroutine opera_init(nw_in,nx_in,nz_in,nreal_in)
    integer	:: nw_in,nx_in,nz_in,nreal_in
    nw=nw_in
    nx=nx_in
    nz=nz_in
    nreal=nreal_in
  end subroutine
!==================================================
  function opera_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 opera_lop2(adj,add,xx,yy )
    stat=0
  end function 
!==================================================
  subroutine opera_lop2(adj,add,xx,yy)
    logical,intent(in)   :: adj,add 
    real, dimension (:)  :: xx 
    real, dimension (:)  :: yy 
    real pi
    integer i, j, k, iw
    real theta
  
    do j = 1, nx*nz
       k = (j-1)*nx*nz
       do i = 1, nx*nz
          if ( adj) then
             xx((j-1)*nw+1:j*nw) = xx((j-1)*nw+1:j*nw) + sum(            yy(k+i) * xx((i-1)*nw+1:i*nw))
          else
             yy(k+i)            =  yy(k+i)             + sum(xx((j-1)*nw+1:j*nw) * xx((i-1)*nw+1:i*nw))
          end if
       end do
    end do
write(0,*) "from opera xx,yy = ",minval(xx),maxval(xx),minval(yy),maxval(yy)
  end subroutine 
!          do iw=1,nw
!             if ( adj) then
!                xx((j-1)*nw+iw)=xx((j-1)*nw+iw)+yy(k) * xx((i-1)*nw+iw)
!             else
!                yy(k)=yy(k)+xx((j-1)*nw+iw) * xx((i-1)*nw+iw)
!          end do
!=================================================
end module opera
