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, x0) result(stat)
    integer            :: stat 
    logical,intent(in) :: adj,add 
    real,dimension(:)  :: xx,yy,x0
    call adjnull (adj,add,xx,yy )
    call opera_lop2(adj,add,xx,yy,x0 )
    stat=0
  end function 
!==================================================
  subroutine opera_lop2(adj,add,xx,yy,x0)
    logical,intent(in)   :: adj,add 
    real, dimension (:)  :: xx 
    real, dimension (:)  :: yy 
    real, dimension (:)  :: x0 
    integer i, j, k
    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) * x0((i-1)*nw+1:i*nw))
             xx((i-1)*nw+1:i*nw) = xx((i-1)*nw+1:i*nw) + sum(            yy(k+i) * x0((j-1)*nw+1:j*nw))
          else
             yy(k+i)            =  yy(k+i)             + sum(xx((i-1)*nw+1:i*nw) * x0((j-1)*nw+1:j*nw))
!             yy(k+i)            =  yy(k+i)             + sum(xx((j-1)*nw+1:j*nw) * x0((i-1)*nw+1:i*nw))
          end if
       end do
    end do
  end subroutine 
end module opera
