module int2 {                                      # generic 2-d interpolation
  integer :: nd
  integer, dimension (2) :: n, nf
  integer, dimension (nd,2),   allocatable :: nxy
  logical, dimension (nd),     allocatable :: m
  real,    dimension (nf(1),nf(2),nd), allocatable :: w
#%  _init (coord, o, d, n, interp, nf, nd, weight)
    real,    dimension (:,:), intent (in)           :: coord
    real,    dimension (:),  optional :: weight
    real, dimension (2) ::     o, d, x
    integer  id, nd, ix (2), stat
    interface
       function interp (x, w) result (stat)
         integer                              :: stat
         real, dimension (2), intent (in)     :: x
         real, dimension (:,:)                :: w
       end function interp
    end interface
     do id = 1, nd {
       x = (coord (id,:) - o)/d ; ix = x
       x = x - ix ; nxy (id,:) = ix + 1 - 0.5*nf
       m (id) = .true. ; w (:,:, id) = 0. 
       if (  all((nxy (id,:) + 1  >= 1)) _
          && all((nxy (id,:) + nf <= n))) {
          m (id) = .false. ; stat = interp (x, w (:,:,id))
          if (present (weight)) w (:,:,id) = w (:,:,id) * weight (id)
       }
    }
#%  _lop (x (n(1), n(2)), ord (:)) 
    integer                          :: i
    integer, dimension (2)           :: i1, i2 
  do i = 1, size (ord) { if (m (i)) cycle
       i1 = nxy (i,:) + 1
       i2 = nxy (i,:) + nf
       if( adj) 
          x (i1(1):i2(1),i1(2):i2(2)) = &
          x (i1(1):i2(1),i1(2):i2(2)) +                ord (i) * w (:,:,i)
       else 
          ord (i) = ord (i) + sum (x (i1(1):i2(1),i1(2):i2(2)) * w (:,:,i))
   }
}