# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "<stdin>"
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