module int2  
! generic 2-d interpolation
  use adj_mod
  implicit none
  integer, private :: nd 
  integer, private :: n1 
  integer, private :: n2 
  integer ,private,  dimension (2)  :: n 
  integer ,private,  dimension (2)  :: nf 
  integer ,private,  dimension (:,:),   allocatable  :: nxy 
  logical ,private,  dimension (:),     allocatable  :: m 
  real ,private,     dimension (:,:,:), allocatable  :: w 
  contains
  subroutine int2_init ( coord,o,d,n_in,interp,nf_in,nd_in,weight )
    integer ,   dimension (2) :: n_in 
    integer ,   dimension (2) :: nf_in 
    integer    :: nd_in 
    real,    dimension (:,:), intent (in)           :: coord
    real,    dimension (:),  optional :: weight
    real, dimension (2) ::     o, d, x
    integer  id, ix (2), stat
    interface  
      integer function interp (x, w)
        real, dimension (2), intent (in)     :: x
        real, dimension (:,:)                :: w
      end function
    end interface 
    n = n_in 
    nf = nf_in 
    nd = nd_in 
    if (.not. allocated(nxy)) then
      allocate(nxy ( nd,2)) 
    end if
    if (.not. allocated(m)) then
      allocate(m ( nd)) 
    end if
    if (.not. allocated(w)) then
      allocate(w ( nf(1),nf(2),nd)) 
    end if
    n1=n(1)
    n2=n(2)
    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))           .and. all((nxy (id,:&
        &) + nf <= n))) then
        m (id) = .false. 
        stat = interp (x, w (:,:,id))
        if (present (weight)) then
          w (:,:,id) = w (:,:,id) * weight (id)
        end if
      end if
    end do
  end subroutine
  function int2_lop ( adj, add, x, ord) result(stat)
    integer            :: stat 
    logical,intent(in) :: adj,add 
    real,dimension(:)  :: x,ord 
    call adjnull (adj,add,x,ord )
    call int2_lop2(adj,add,x,ord )
    stat=0
  end function 
  subroutine int2_lop2(adj,add,x,ord)
    logical,intent(in) :: adj,add 
    real, dimension (n1, n2)  :: x 
    real, dimension (:)  :: ord 
    integer                          :: i
    integer, dimension (2)           :: i1, i2 
    do i = 1, size (ord)
      if (m (i)) then
        cycle
      end if
      i1 = nxy (i,:) + 1
      i2 = nxy (i,:) + nf
      if ( adj) then
        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))
      end if
    end do
  end subroutine 
  subroutine int2_close()
    deallocate( nxy ,m,w)
  end subroutine 
end module 
