!/*@@
!
! $Id: inter3dp.f90,v 1.1 1996/12/16 17:51:47 scott Exp $
!
!@@*/
!
! INPUT:
!             h : Physical spatial discretization.
!      dx,dy,dz : The change in the x-tilde coordinate value for each point.
!       srcMask : The Mask for the uniform grid in x-tilde.
!       dstMask : The Mask for the uni uniform grid in x-tilde.
!       dstMask : The Mask for the uniform grid in x.
!
! OUTPUT:
!        shifts : Kind of point interpolation to be used.
!                 NOPOINT -> the point should not be filled.
!                 NOSHIFT -> the point can be interpolated without shifting
!                 positive -> ind interpolated without shifting
!                 positive -> index into shift arrays to get necessary shift
!                             value
!
!--------------------------------------------------------------------------
subroutine preshift2(h,dx,dy,dz,srcMask,dstMask,shifts,ex,shift_p1, &
                    shift_p2,shift_p3,num_int_points)
use input_params
  implicit none
  integer Nshifts
  parameter (Nshifts = 728)
  real*8,  dimension(3)                 :: h
  integer, dimension(3)                 :: ex
  real*8,  dimension(:,:,:) :: dx,dy,dz,srcMask,dstMask,shifts
  integer, dimension(:)    :: shift_p1,shift_p2,shift_p3 !1D array of points
  integer :: num_int_points  ! the total # of points to shift
  intent(IN) h,dx,dy,dz,srcMask,dstMask,ex
  !
  integer :: i,j,k,test,n
  integer, dimension(3) :: pnt,ind,shift
  real*8, dimension(3) :: hm1,dxs
  real*8 :: INVALID,ONE
  parameter(INVALID = 0.d0,ONE = 1.d0)
  integer :: NOPOINT,NOSHIFT
  parameter(NOPOINT = -1,NOSHIFT = 0)
  !
  hm1 = ONE/h
  num_int_points = 0
  !

  do k=1,ex(3)
    do j=1,ex(2)
      do i=1,ex(1)
	      if(dstMask(i,j,k) /= INVALID) then
	         pnt = (/ i, j, k /)
	         dxs = hm1* (/ dx(i,j,k), dy(i,j,k), dz(i,j,k) /)
	         ind = FLOOR(dxs) + pnt - 1
	         where(ind < 1) ind = 1       ! Find the starting point index for
	         where(ind + 3 > ex) ind = ex - 3  ! 4 point interpolation
            if ((ind(1)<1).or.(ind(2)<1).or.(ind(3)<1).or.(ind(1)+3>ex(1)).or. &
               (ind(2)+3>ex(2)).or.(ind(3)+3>ex(3)) ) then
               test = 1
            else
               test=int(sum(srcmask(ind(1):ind(1)+3,ind(2):ind(2)+3,ind(3):ind(3)+3) ) )
               if (test.eq.64) then
                  test = 1 
               else
                  test = 0
               end if
            end if

	         if(test == 1) then ! no shifting is required
               continue
	         else
               num_int_points = num_int_points + 1
               shift_p1(num_int_points) = i
               shift_p2(num_int_points) = j
               shift_p3(num_int_points) = k
            end if
         end if
      end do 
    end do
  end do
  !
end subroutine preshift2
