module igrad2 {                       # 2-D gradient with adjoint,  r= grad( p)
integer :: n1, n2
#%_init   (n1, n2)
#%_lop ( p(n1, n2),  r(n1,n2,2))
integer i,j
do i= 1, n1-1 {
do j= 1, n2-1 {
	if( adj) {
		p(i+1,j  )  +=  r(i,j,1) 
                p(i  ,j  )  -=  r(i,j,1) 
                p(i  ,j+1)  +=  r(i,j,2) 
                p(i  ,j  )  -=  r(i,j,2)
		}
	else {  r(i,j,1)  +=  ( p(i+1,j) - p(i,j)) 
                r(i,j,2)  +=  ( p(i,j+1) - p(i,j)) 
		}
	}}
}