module igrad2  
! 2-D gradient with adjoint,  r= grad( p)
  use adj_mod
  implicit none
  integer, private :: n1 
  integer, private :: n2 
  contains
  subroutine igrad2_init ( n1_in,n2_in )
    integer    :: n1_in 
    integer    :: n2_in 
    n1 = n1_in 
    n2 = n2_in
  end subroutine
  function igrad2_lop ( adj, add, p, r) result(stat)
    integer            :: stat 
    logical,intent(in) :: adj,add 
    real,dimension(:)  :: p,r 
    call adjnull (adj,add,p,r )
    call igrad2_lop2(adj,add,p,r )
    stat=0
  end function 
  subroutine igrad2_lop2(adj,add,p,r)
    logical,intent(in) :: adj,add 
    real, dimension (n1, n2)  :: p 
    real, dimension (n1,n2,2)  :: r 
    integer i,j
    do i= 1, n1-1  
      do j= 1, n2-1  
        if ( adj) then
          p(i+1,j  )  =		p(i+1,j  )  +  r(i,j,1) 
          p(i  ,j  )  =                p(i  ,j  )  -  r(i,j,1) 
          p(i  ,j+1)  =                p(i  ,j+1)  +  r(i,j,2) 
          p(i  ,j  )  =                p(i  ,j  )  -  r(i,j,2)
        else
          r(i,j,1)  =   r(i,j,1)  +  ( p(i+1,j) - p(i,j)) 
          r(i,j,2)  =                r(i,j,2)  +  ( p(i,j+1) - p(i,j))
        end if
      end do
    end do
  end subroutine 
  subroutine igrad2_close()
  end subroutine 
end module