#$=head1 NAME
#$
#$igrad2 - 2D gradient operator
#$
#$=head1 SYNOPSIS
#$
#$C<ierr=igrad2_init(n1,n2)>
#$
#$C<ierr=igrad2_lop(adj,add,xx,yy)>
#$
#$=head1 PARAMETERS
#$
#$=over 4
#$
#$=item n1,n2  - integer   
#$
#$      dimensions of dataset
#$
#$=item adj,add,xx,yy - 
#$
#$      standard operator interface
#$
#$=back
#$
#$=head1 DESCRIPTION
#$
#$ Calculate 2-D gradient
#$
#$=head1 SEE ALSO
#$
#$L<igrad1>
#$
#$=head1 LIBRARY
#$
#$B<geef90>
#$
#$=cut
#$

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)) 
		}
	}}
}