!/*@@
!
! $Id: shift2.f90,v 1.7 1997/09/12 23:33:00 mijan Exp $
!
!@@*/
!--------------------------------------------------------------------------
! Interpolation routines for Causal differencing.
!--------------------------------------------------------------------------
!
! Precompute the masked interpolation shifts for interpolating a grid of
! points. Computes the interpolation stencils for the interpolator: make_xgrid.
! This routine will correctly compute stencils for interior points as well as
! mask boundary points. But it *will not* compute stencils correctly for 
! situations where new points come into the domain where the source and
! destination masks are different. An additional routine is required for this
! to preserve the existing structure where mask boundary points can be taken
! out of the interpolation.
!
! 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 preshift(h,dx,dy,dz,srcMask,dstMask,shifts,ex,shift_p1, &
           shift_p2,shift_p3,shift_b,params,num_int_points,m,nl)
use input_params
  implicit none
  integer Nshifts
  parameter (Nshifts = 1330)
  integer, dimension(Nshifts),save :: shiftx,shifty,shiftz
  integer, dimension(:) :: m
  integer               :: nl




data shiftx/1,-1,0,0,0,0,1,-1,1,-1,1,-1,0,0,1,-1,0,0,1,-1,1,-1,1,-1,1,-1,2,-2,&
     0,0,0,0,2,-2,1,-1,2,-2,1,-1,2,-2,0,0,1,-1,0,0,2,-2,0,0,1,-1,0,0,&
     2,-2,1,-1,2,-2,1,-1,1,-1,1,-1,2,-2,1,-1,2,-2,1,-1,1,-1,1,-1,2,-2,2,-2,&
     2,-2,0,0,2,-2,0,0,3,-3,0,0,2,-2,2,-2,2,-2,1,-1,2,-2,1,-1,0,2,-2,2,&
     -2,2,-2,1,-1,2,-2,1,-1,0,3,-3,1,-1,3,-3,1,-1,3,-3,0,0,1,-1,0,0,3,-3,&
     0,0,1,-1,0,0,3,-3,1,-1,3,-3,1,-1,1,-1,1,-1,3,-3,1,-1,3,-3,1,-1,1,-1,&
     1,-1,2,-2,2,-2,2,-2,2,-2,3,-3,2,-2,3,-3,2,-2,3,-3,0,0,2,-2,0,0,3,-3,&
     0,0,2,-2,0,0,3,-3,2,-2,3,-3,2,-2,3,-3,1,-1,3,-3,1,-1,2,-2,1,-1,2,-2,&
     1,-1,3,-3,2,-2,3,-3,2,-2,3,-3,1,-1,3,-3,1,-1,2,-2,1,-1,2,-2,1,-1,4,-4,&
     0,0,0,0,4,-4,1,-1,4,-4,1,-1,4,-4,0,0,3,-3,2,-2,3,-3,2,-2,2,-2,2,-2,&
     1,-1,0,0,4,-4,0,0,3,-3,2,-2,3,-3,2,-2,2,-2,2,-2,1,-1,0,0,3,-3,3,-3,&
     4,-4,1,-1,4,-4,1,-1,3,-3,0,0,1,-1,1,-1,4,-4,1,-1,4,-4,1,-1,3,-3,0,0,&
     1,-1,1,-1,3,-3,3,-3,3,-3,1,-1,3,-3,1,-1,3,-3,3,-3,3,-3,1,-1,3,-3,1,-1,&
     4,-4,2,-2,4,-4,2,-2,4,-4,0,0,2,-2,0,0,4,-4,0,0,2,-2,0,0,4,-4,2,-2,&
     4,-4,2,-2,4,-4,1,-1,4,-4,1,-1,2,-2,1,-1,2,-2,1,-1,4,-4,2,-2,4,-4,2,-2,&
     4,-4,1,-1,4,-4,1,-1,2,-2,1,-1,2,-2,1,-1,3,-3,3,-3,3,-3,2,-2,3,-3,2,-2,&
     3,-3,3,-3,3,-3,2,-2,3,-3,2,-2,4,-4,2,-2,4,-4,2,-2,2,-2,2,-2,4,-4,2,-2,&
     4,-4,2,-2,2,-2,2,-2,5,-5,4,-4,3,-3,0,4,-4,3,-3,0,4,-4,0,0,3,-3,0,0,&
     0,4,-4,0,0,3,-3,0,0,0,5,-5,1,-1,5,-5,1,-1,5,-5,4,-4,3,-3,0,4,-4,3,&
     -3,0,4,-4,1,-1,4,-4,1,-1,3,-3,1,-1,3,-3,1,-1,1,-1,0,0,5,-5,4,-4,3,-3,&
     0,4,-4,3,-3,0,4,-4,1,-1,4,-4,1,-1,3,-3,1,-1,3,-3,1,-1,1,-1,0,0,5,-5,&
     1,-1,5,-5,1,-1,3,-3,3,-3,1,-1,1,-1,5,-5,1,-1,5,-5,1,-1,3,-3,3,-3,1,-1,&
     1,-1,5,-5,2,-2,5,-5,2,-2,5,-5,4,-4,3,-3,0,4,-4,3,-3,0,4,-4,2,-2,4,-4,&
     2,-2,3,-3,2,-2,3,-3,2,-2,2,-2,0,0,5,-5,4,-4,3,-3,0,4,-4,3,-3,0,4,-4,&
     2,-2,4,-4,2,-2,3,-3,2,-2,3,-3,2,-2,2,-2,0,0,5,-5,2,-2,5,-5,2,-2,5,-5,&
     1,-1,5,-5,1,-1,2,-2,1,-1,2,-2,1,-1,5,-5,2,-2,5,-5,2,-2,5,-5,1,-1,5,-5,&
     1,-1,2,-2,1,-1,2,-2,1,-1,4,-4,4,-4,4,-4,0,0,4,-4,0,0,4,-4,4,-4,5,-5,&
     2,-2,5,-5,2,-2,4,-4,1,-1,4,-4,1,-1,2,-2,2,-2,4,-4,4,-4,5,-5,2,-2,5,-5,&
     2,-2,4,-4,1,-1,4,-4,1,-1,2,-2,2,-2,5,-5,3,-3,5,-5,3,-3,5,-5,4,-4,3,-3,&
     0,4,-4,3,-3,0,3,-3,3,-3,3,-3,0,0,5,-5,4,-4,3,-3,0,4,-4,3,-3,0,3,-3,&
     3,-3,3,-3,0,0,5,-5,3,-3,5,-5,3,-3,5,-5,1,-1,5,-5,1,-1,3,-3,1,-1,3,-3,&
     1,-1,5,-5,3,-3,5,-5,3,-3,5,-5,1,-1,5,-5,1,-1,3,-3,1,-1,3,-3,1,-1,4,-4,&
     4,-4,4,-4,2,-2,4,-4,2,-2,4,-4,4,-4,4,-4,2,-2,4,-4,2,-2,5,-5,3,-3,5,-5,&
     3,-3,5,-5,2,-2,5,-5,2,-2,3,-3,2,-2,3,-3,2,-2,5,-5,3,-3,5,-5,3,-3,5,-5,&
     2,-2,5,-5,2,-2,3,-3,2,-2,3,-3,2,-2,5,-5,4,-4,5,-5,4,-4,4,-4,4,-4,5,-5,&
     4,-4,3,-3,0,4,-4,3,-3,0,4,-4,0,0,4,-4,4,-4,5,-5,4,-4,3,-3,0,4,-4,3,&
     -3,0,4,-4,0,0,5,-5,4,-4,5,-5,4,-4,5,-5,1,-1,5,-5,1,-1,4,-4,1,-1,4,-4,&
     1,-1,5,-5,4,-4,5,-5,4,-4,5,-5,1,-1,5,-5,1,-1,4,-4,1,-1,4,-4,1,-1,5,-5,&
     3,-3,5,-5,3,-3,3,-3,3,-3,5,-5,3,-3,5,-5,3,-3,3,-3,3,-3,5,-5,4,-4,5,-5,&
     4,-4,5,-5,2,-2,5,-5,2,-2,4,-4,2,-2,4,-4,2,-2,5,-5,4,-4,5,-5,4,-4,5,-5,&
     2,-2,5,-5,2,-2,4,-4,2,-2,4,-4,2,-2,4,-4,4,-4,4,-4,4,-4,5,-5,5,-5,5,-5,&
     4,-4,5,-5,4,-4,5,-5,3,-3,5,-5,3,-3,5,-5,4,-4,3,-3,0,4,-4,3,-3,0,5,-5,&
     4,-4,5,-5,4,-4,5,-5,3,-3,5,-5,3,-3,5,-5,4,-4,3,-3,0,4,-4,3,-3,0,5,-5,&
     5,-5,5,-5,1,-1,5,-5,1,-1,5,-5,5,-5,5,-5,1,-1,5,-5,1,-1,5,-5,5,-5,5,-5,&
     2,-2,5,-5,2,-2,5,-5,5,-5,5,-5,2,-2,5,-5,2,-2,5,-5,4,-4,5,-5,4,-4,4,-4,&
     4,-4,5,-5,4,-4,5,-5,4,-4,4,-4,4,-4,5,-5,5,-5,5,-5,3,-3,5,-5,3,-3,5,-5,&
     5,-5,5,-5,3,-3,5,-5,3,-3,5,-5,5,-5,5,-5,4,-4,5,-5,4,-4,5,-5,5,-5,5,-5,&
     4,-4,5,-5,4,-4,5,-5,5,-5,5,-5,5,-5/
data shifty/0,0,1,-1,0,0,1,1,-1,-1,0,0,1,-1,0,0,1,-1,1,1,-1,-1,1,1,-1,-1,0,0,&
     2,-2,0,0,1,1,2,2,-1,-1,-2,-2,0,0,2,-2,0,0,1,-1,0,0,2,-2,0,0,1,-1,&
     1,1,2,2,-1,-1,-2,-2,1,1,-1,-1,1,1,2,2,-1,-1,-2,-2,1,1,-1,-1,2,2,-2,-2,&
     0,0,2,-2,0,0,2,-2,0,0,3,-3,2,2,-2,-2,1,1,2,2,-1,-1,-2,-2,0,2,2,-2,&
     -2,1,1,2,2,-1,-1,-2,-2,0,1,1,3,3,-1,-1,-3,-3,0,0,3,-3,0,0,1,-1,0,0,&
     3,-3,0,0,1,-1,1,1,3,3,-1,-1,-3,-3,1,1,-1,-1,1,1,3,3,-1,-1,-3,-3,1,1,&
     -1,-1,2,2,-2,-2,2,2,-2,-2,2,2,3,3,-2,-2,-3,-3,0,0,3,-3,0,0,2,-2,0,0,&
     3,-3,0,0,2,-2,2,2,3,3,-2,-2,-3,-3,1,1,3,3,-1,-1,-3,-3,1,1,2,2,-1,-1,&
     -2,-2,2,2,3,3,-2,-2,-3,-3,1,1,3,3,-1,-1,-3,-3,1,1,2,2,-1,-1,-2,-2,0,0,&
     4,-4,0,0,1,1,4,4,-1,-1,-4,-4,0,0,4,-4,2,2,3,3,-2,-2,-3,-3,2,2,-2,-2,&
     0,0,1,-1,0,0,4,-4,2,2,3,3,-2,-2,-3,-3,2,2,-2,-2,0,0,1,-1,3,3,-3,-3,&
     1,1,4,4,-1,-1,-4,-4,0,0,3,-3,1,1,-1,-1,1,1,4,4,-1,-1,-4,-4,0,0,3,-3,&
     1,1,-1,-1,3,3,-3,-3,1,1,3,3,-1,-1,-3,-3,3,3,-3,-3,1,1,3,3,-1,-1,-3,-3,&
     2,2,4,4,-2,-2,-4,-4,0,0,4,-4,0,0,2,-2,0,0,4,-4,0,0,2,-2,2,2,4,4,&
     -2,-2,-4,-4,1,1,4,4,-1,-1,-4,-4,1,1,2,2,-1,-1,-2,-2,2,2,4,4,-2,-2,-4,-4,&
     1,1,4,4,-1,-1,-4,-4,1,1,2,2,-1,-1,-2,-2,3,3,-3,-3,2,2,3,3,-2,-2,-3,-3,&
     3,3,-3,-3,2,2,3,3,-2,-2,-3,-3,2,2,4,4,-2,-2,-4,-4,2,2,-2,-2,2,2,4,4,&
     -2,-2,-4,-4,2,2,-2,-2,0,0,3,3,4,4,5,-3,-3,-4,-4,-5,0,0,4,-4,0,0,3,-3,&
     0,0,0,4,-4,0,0,3,-3,0,1,1,5,5,-1,-1,-5,-5,0,0,3,3,4,4,5,-3,-3,-4,&
     -4,-5,1,1,4,4,-1,-1,-4,-4,1,1,3,3,-1,-1,-3,-3,0,0,1,-1,0,0,3,3,4,4,&
     5,-3,-3,-4,-4,-5,1,1,4,4,-1,-1,-4,-4,1,1,3,3,-1,-1,-3,-3,0,0,1,-1,1,1,&
     5,5,-1,-1,-5,-5,3,3,-3,-3,1,1,-1,-1,1,1,5,5,-1,-1,-5,-5,3,3,-3,-3,1,1,&
     -1,-1,2,2,5,5,-2,-2,-5,-5,0,0,3,3,4,4,5,-3,-3,-4,-4,-5,2,2,4,4,-2,-2,&
     -4,-4,2,2,3,3,-2,-2,-3,-3,0,0,2,-2,0,0,3,3,4,4,5,-3,-3,-4,-4,-5,2,2,&
     4,4,-2,-2,-4,-4,2,2,3,3,-2,-2,-3,-3,0,0,2,-2,2,2,5,5,-2,-2,-5,-5,1,1,&
     5,5,-1,-1,-5,-5,1,1,2,2,-1,-1,-2,-2,2,2,5,5,-2,-2,-5,-5,1,1,5,5,-1,-1,&
     -5,-5,1,1,2,2,-1,-1,-2,-2,4,4,-4,-4,0,0,4,-4,0,0,4,-4,4,4,-4,-4,2,2,&
     5,5,-2,-2,-5,-5,1,1,4,4,-1,-1,-4,-4,2,2,-2,-2,4,4,-4,-4,2,2,5,5,-2,-2,&
     -5,-5,1,1,4,4,-1,-1,-4,-4,2,2,-2,-2,3,3,5,5,-3,-3,-5,-5,0,0,3,3,4,4,&
     5,-3,-3,-4,-4,-5,3,3,-3,-3,0,0,3,-3,0,0,3,3,4,4,5,-3,-3,-4,-4,-5,3,3,&
     -3,-3,0,0,3,-3,3,3,5,5,-3,-3,-5,-5,1,1,5,5,-1,-1,-5,-5,1,1,3,3,-1,-1,&
     -3,-3,3,3,5,5,-3,-3,-5,-5,1,1,5,5,-1,-1,-5,-5,1,1,3,3,-1,-1,-3,-3,4,4,&
     -4,-4,2,2,4,4,-2,-2,-4,-4,4,4,-4,-4,2,2,4,4,-2,-2,-4,-4,3,3,5,5,-3,-3,&
     -5,-5,2,2,5,5,-2,-2,-5,-5,2,2,3,3,-2,-2,-3,-3,3,3,5,5,-3,-3,-5,-5,2,2,&
     5,5,-2,-2,-5,-5,2,2,3,3,-2,-2,-3,-3,4,4,5,5,-4,-4,-5,-5,4,4,-4,-4,0,0,&
     3,3,4,4,5,-3,-3,-4,-4,-5,0,0,4,-4,4,4,-4,-4,0,0,3,3,4,4,5,-3,-3,-4,&
     -4,-5,0,0,4,-4,4,4,5,5,-4,-4,-5,-5,1,1,5,5,-1,-1,-5,-5,1,1,4,4,-1,-1,&
     -4,-4,4,4,5,5,-4,-4,-5,-5,1,1,5,5,-1,-1,-5,-5,1,1,4,4,-1,-1,-4,-4,3,3,&
     5,5,-3,-3,-5,-5,3,3,-3,-3,3,3,5,5,-3,-3,-5,-5,3,3,-3,-3,4,4,5,5,-4,-4,&
     -5,-5,2,2,5,5,-2,-2,-5,-5,2,2,4,4,-2,-2,-4,-4,4,4,5,5,-4,-4,-5,-5,2,2,&
     5,5,-2,-2,-5,-5,2,2,4,4,-2,-2,-4,-4,4,4,-4,-4,4,4,-4,-4,5,5,-5,-5,4,4,&
     5,5,-4,-4,-5,-5,3,3,5,5,-3,-3,-5,-5,0,0,3,3,4,4,5,-3,-3,-4,-4,-5,4,4,&
     5,5,-4,-4,-5,-5,3,3,5,5,-3,-3,-5,-5,0,0,3,3,4,4,5,-3,-3,-4,-4,-5,5,5,&
     -5,-5,1,1,5,5,-1,-1,-5,-5,5,5,-5,-5,1,1,5,5,-1,-1,-5,-5,5,5,-5,-5,2,2,&
     5,5,-2,-2,-5,-5,5,5,-5,-5,2,2,5,5,-2,-2,-5,-5,4,4,5,5,-4,-4,-5,-5,4,4,&
     -4,-4,4,4,5,5,-4,-4,-5,-5,4,4,-4,-4,5,5,-5,-5,3,3,5,5,-3,-3,-5,-5,5,5,&
     -5,-5,3,3,5,5,-3,-3,-5,-5,5,5,-5,-5,4,4,5,5,-4,-4,-5,-5,5,5,-5,-5,4,4,&
     5,5,-4,-4,-5,-5,5,5,-5,-5,5,5,-5,-5/
data shiftz/0,0,0,0,1,-1,0,0,0,0,1,1,1,1,-1,-1,-1,-1,&
     1,1,1,1,-1,-1,-1,-1,0,0,&
     0,0,2,-2,0,0,0,0,0,0,0,0,1,1,1,1,2,2,2,2,-1,-1,-1,-1,-2,-2,-2,-2,&
     1,1,1,1,1,1,1,1,2,2,2,2,-1,-1,-1,-1,-1,-1,-1,-1,-2,-2,-2,-2,0,0,0,0,&
     2,2,2,2,-2,-2,-2,-2,0,0,0,0,1,1,1,1,2,2,2,2,2,2,2,2,3,-1,-1,-1,&
     -1,-2,-2,-2,-2,-2,-2,-2,-2,-3,0,0,0,0,0,0,0,0,1,1,1,1,3,3,3,3,-1,-1,&
     -1,-1,-3,-3,-3,-3,1,1,1,1,1,1,1,1,3,3,3,3,-1,-1,-1,-1,-1,-1,-1,-1,-3,-3,&
     -3,-3,2,2,2,2,-2,-2,-2,-2,0,0,0,0,0,0,0,0,2,2,2,2,3,3,3,3,-2,-2,&
     -2,-2,-3,-3,-3,-3,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3,3,&
     3,3,-1,-1,-1,-1,-1,-1,-1,-1,-2,-2,-2,-2,-2,-2,&
     -2,-2,-3,-3,-3,-3,-3,-3,-3,-3,0,0,&
     0,0,4,-4,0,0,0,0,0,0,0,0,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,&
     4,4,4,4,-1,-1,-1,-1,-2,-2,-2,-2,-2,-2,-2,&
     -2,-3,-3,-3,-3,-4,-4,-4,-4,0,0,0,0,&
     1,1,1,1,1,1,1,1,3,3,3,3,4,4,4,4,-1,-1,-1,-1,-1,-1,-1,-1,-3,-3,-3,-3,&
     -4,-4,-4,-4,1,1,1,1,3,3,3,3,3,3,3,3,-1,-1,-1,-1,-3,-3,-3,-3,-3,-3,-3,-3,&
     0,0,0,0,0,0,0,0,2,2,2,2,4,4,4,4,-2,-2,-2,-2,-4,-4,-4,-4,1,1,1,1,&
     1,1,1,1,2,2,2,2,2,2,2,2,4,4,4,4,4,4,4,4,-1,-1,-1,-1,-1,-1,-1,-1,&
     -2,-2,-2,-2,-2,-2,-2,-2,-4,-4,-4,-4,-4,-4,-4,-4,2,2,2,2,3,3,3,3,3,3,3,3,&
     -2,-2,-2,-2,-3,-3,-3,-3,-3,-3,-3,-3,2,2,2,2,2,2,2,2,4,4,4,4,-2,-2,-2,-2,&
     -2,-2,-2,-2,-4,-4,-4,-4,0,0,0,0,0,0,0,0,0,0,0,0,3,3,3,3,4,4,4,4,&
     5,-3,-3,-3,-3,-4,-4,-4,-4,-5,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,&
     1,1,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,5,5,5,5,-1,-1,-1,-1,-1,-1,&
     -1,-1,-1,-1,-1,-1,-3,-3,-3,-3,-3,-3,-3,-3,-4,&
     -4,-4,-4,-4,-4,-4,-4,-5,-5,-5,-5,1,1,&
     1,1,1,1,1,1,3,3,3,3,5,5,5,5,-1,-1,-1,-1,-1,-1,-1,-1,-3,-3,-3,-3,-5,-5,&
     -5,-5,0,0,0,0,0,0,0,0,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,&
     3,3,4,4,4,4,4,4,4,4,5,5,5,5,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-3,-3,&
     -3,-3,-3,-3,-3,-3,-4,-4,-4,-4,-4,-4,-4,-4,&
     -5,-5,-5,-5,1,1,1,1,1,1,1,1,2,2,&
     2,2,2,2,2,2,5,5,5,5,5,5,5,5,-1,-1,-1,-1,-1,-1,-1,-1,-2,-2,-2,-2,-2,-2,&
     -2,-2,-5,-5,-5,-5,-5,-5,-5,-5,0,0,0,0,4,4,4,4,-4,-4,-4,-4,1,1,1,1,2,2,&
     2,2,2,2,2,2,4,4,4,4,4,4,4,4,5,5,5,5,-1,-1,-1,-1,-2,-2,-2,-2,-2,-2,&
     -2,-2,-4,-4,-4,-4,-4,-4,-4,-4,-5,-5,-5,-5,0,0,0,0,0,0,0,0,3,3,3,3,3,3,&
     3,3,3,3,3,3,4,4,4,4,5,5,5,5,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-4,-4,&
     -4,-4,-5,-5,-5,-5,1,1,1,1,1,1,1,1,3,3,3,3,3,3,3,3,5,5,5,5,5,5,&
     5,5,-1,-1,-1,-1,-1,-1,-1,-1,-3,-3,-3,-3,-3,-3,-3,-3,-5,&
     -5,-5,-5,-5,-5,-5,-5,2,2,&
     2,2,4,4,4,4,4,4,4,4,-2,-2,-2,-2,-4,-4,-4,-4,-4,-4,-4,-4,2,2,2,2,2,2,&
     2,2,3,3,3,3,3,3,3,3,5,5,5,5,5,5,5,5,-2,-2,-2,-2,-2,-2,-2,-2,-3,-3,&
     -3,-3,-3,-3,-3,-3,-5,-5,-5,-5,-5,-5,-5,-5,0,0,0,0,0,0,0,0,3,3,3,3,4,4,&
     4,4,4,4,4,4,4,4,4,4,5,5,5,5,-3,-3,-3,-3,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,&
     -4,-4,-5,-5,-5,-5,1,1,1,1,1,1,1,1,4,4,4,4,4,4,4,4,5,5,5,5,5,5,&
     5,5,-1,-1,-1,-1,-1,-1,-1,-1,-4,-4,-4,-4,-4,-4,-4,&
     -4,-5,-5,-5,-5,-5,-5,-5,-5,3,3,&
     3,3,3,3,3,3,5,5,5,5,-3,-3,-3,-3,-3,-3,-3,-3,-5,-5,-5,-5,2,2,2,2,2,2,&
     2,2,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,-2,-2,-2,-2,-2,-2,-2,-2,-4,-4,&
     -4,-4,-4,-4,-4,-4,-5,-5,-5,-5,-5,-5,-5,&
     -5,4,4,4,4,-4,-4,-4,-4,0,0,0,0,3,3,&
     3,3,3,3,3,3,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,-3,-3,&
     -3,-3,-3,-3,-3,-3,-4,-4,-4,-4,-4,-4,-4,-4,-5,&
     -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,1,1,&
     1,1,5,5,5,5,5,5,5,5,-1,-1,-1,-1,-5,-5,-5,-5,-5,-5,-5,-5,2,2,2,2,5,5,&
     5,5,5,5,5,5,-2,-2,-2,-2,-5,-5,-5,-5,-5,-5,-5,-5,4,4,4,4,4,4,4,4,5,5,&
     5,5,-4,-4,-4,-4,-4,-4,-4,-4,-5,-5,-5,-5,3,3,3,3,5,5,5,5,5,5,5,5,-3,-3,&
     -3,-3,-5,-5,-5,-5,-5,-5,-5,-5,4,4,4,4,5,5,5,5,5,5,5,5,-4,-4,-4,-4,-5,-5,&
     -5,-5,-5,-5,-5,-5,5,5,5,5,-5,-5,-5,-5/
  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, dimension(:)    :: shift_b
  integer :: num_int_points  ! the total # of points to shift
   real*8 :: EPSRND,dxeps(3),small
   parameter (EPSRND = 1.0d6,small=1.0d-08)

  !
  integer :: i,j,k,test,n,l
  integer, dimension(3) :: pnt,ind,shift
  real*8, dimension(3) :: hm1,dxs
  real*8 :: INVALID,ONE

!----------------------------------------------------------
  real*8, dimension(:) :: params
  real*8, dimension(1025) :: x,y,z    !temporary storage
  integer li,lj,lk,bias
  real*8    :: lx,ly,lz,dx1,dy1,dz1
!----------------------------------------------------------

  parameter(INVALID = 0.d0,ONE = 1.d0)
  integer :: NOPOINT,NOSHIFT
  parameter(NOPOINT = -1,NOSHIFT = 0)
  !
  hm1 = ONE/h
  num_int_points = 0
  do li = 1 , int(params(P_NX))
    x(li) = params(P_XMIN) + h(1)*(li-1)
  end do
  do lj = 1 , int(params(P_NY))
    y(lj) = params(P_YMIN) + h(2)*(lj-1)
  end do
  do lk = 1 , int(params(P_NZ))
    z(lk) = params(P_ZMIN) + h(3)*(lk-1)
  end do


  ! First set up the interior of the computational domain.
  do k=1,ex(3)
    do j=1,ex(2)
      pointloop: do i=1,ex(1)
   if(dstMask(i,j,k) /= INVALID) then
     pnt = (/ i, j, k /)

     if (abs(dx(i,j,k)).le.small) dx(i,j,k) = 0.0d0
     if (abs(dy(i,j,k)).le.small) dy(i,j,k) = 0.0d0
     if (abs(dz(i,j,k)).le.small) dz(i,j,k) = 0.0d0

     dxs = hm1* (/ dx(i,j,k), dy(i,j,k), dz(i,j,k) /)
     dxeps = EPSRND*SPACING(dxs)
     ind = FLOOR(dxs+dxeps) + pnt - 1


     where(ind < 1) ind = 1        ! Find the starting point index for
     where(ind + 3 > ex) ind = ex - 3  ! 4 point interpolation
     !
     shift = 0
     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.or.( (abs(dx(i,j,k)).le.small).and.(abs(dy(i,j,k)) &
      .le.small).and.(abs(dz(i,j,k)).le.small) ) ) then
      test = 1 
        else
      test = 0
        end if
     end if

     !
     if(test == 1) then ! no shifting is required
!       shifts(i,j,k) = TRANSFER(NOSHIFT,ONE)
       shifts(i,j,k) = ONE * NOSHIFT 

     else
       ind = FLOOR(dxs+dxeps) + pnt - 1  !scott put dxeps here
!---------------------------------------------------------------------------
! here we will determine the bias in a particular direction, this says
! that we have a hole

!  shift..1 is biased in x, shift..2 is in y, shift..3 is in z
! this assumes the hole is centered on the grid

       lx = x(i) + dx(i,j,k) 
       ly = y(j) + dy(i,j,k) 
       lz = z(k) + dz(i,j,k) 
       dx1  = abs(lx-params(P_HOLE1_X))
       dy1  = abs(ly-params(P_HOLE1_Y))
       dz1  = abs(lz-params(P_HOLE1_Z))
       if (abs(abs(dx1)-abs(dy1)).le.small) dy1 = dx1
       if (abs(abs(dx1)-abs(dz1)).le.small) dz1 = dx1
       if (abs(abs(dz1)-abs(dy1)).le.small) dz1 = dy1


		 if (dx1.eq.dy1.and.dx1.eq.dz1) then
	       bias = 12
       else if (dx1.gt.dy1.and.dx1.gt.dz1) then
          if (dy1.gt.dz1) then
             bias = 12
          else  if (dz1.gt.dy1) then
             bias = 13
			 else  if (lz.lt.ly) then  !dy=dz
             bias = 12
			 else
             bias = 13
          end if
       else if (dy1.gt.dz1) then
          if (dx1.gt.dz1) then
             bias = 21
          else  if (dz1.gt.dx1) then
             bias = 23
			 else if (lx.lt.lz) then !dx=dz
             bias = 23
          else
             bias = 21
          end if
      else
          if (dx1.gt.dy1) then
             bias = 31
          else  if (dy1.gt.dx1) then
             bias = 32
			 else if (ly.lt.lx) then !dy=dx
             bias = 31
          else
             bias = 32
          end if
       end if

!---------------------------------------------------------------------------
       do n=1,Nshifts

       select case(bias)
     case(12)
         shift(1) = shiftx(n)
         shift(2) = shifty(n)
         shift(3) = shiftz(n)
     case(13)
         shift(1) = shiftx(n)
         shift(2) = shiftz(n)
         shift(3) = shifty(n)
     case(21)
         shift(1) = shifty(n)
         shift(2) = shiftx(n)
         shift(3) = shiftz(n)
     case(23)
         shift(1) = shiftz(n)
         shift(2) = shiftx(n)
         shift(3) = shifty(n)
     case(31)
         shift(1) = shifty(n)
         shift(2) = shiftz(n)
         shift(3) = shiftx(n)
     case(32)
         shift(1) = shiftz(n)
         shift(2) = shifty(n)
         shift(3) = shiftx(n)
       end select

       
     if (ind(1)+shift(1)<1.or.ind(2)+shift(2)<1.or.ind(3)+shift(3)<1.or.   &
      ind(1)+3+shift(1)>ex(1).or.ind(2)+3+shift(2)>ex(2).or.   &
      ind(3)+3+shift(3)>ex(3)) then
        test = 1
     else
        test = int(sum(srcmask(ind(1)+shift(1):ind(1)+shift(1)+3,   &
        ind(2)+shift(2):ind(2)+shift(2)+3,ind(3)+shift(3):ind(3)+shift(3)+3)))

        if (test.eq.64) then
          test = 1 
        else
          test = 0
        end if
     end if

     if(test == 1) then
         shifts(i,j,k) = n*ONE
         num_int_points = num_int_points + 1
         shift_b(num_int_points) = bias 
         shift_p1(num_int_points) = i
         shift_p2(num_int_points) = j 
         shift_p3(num_int_points) = k 
         cycle pointloop
     end if
  end do
22 format(7i5)


!---------------------------------------------------------------------------

       write(*,*) "WARNING!! : PRESHIFT - Cannot interpolate to x grid!"
       shifts(i,j,k) = NOPOINT*ONE
     end if
   else
!     shifts(i,j,k) = TRANSFER(NOPOINT,ONE)
     shifts(i,j,k) = NOPOINT*ONE
   end if
      end do pointloop
    end do
  end do
  !
return
end subroutine preshift


!--------------------------------------------------------------------------
!
! Interpolate the data to the uniform x grid as defined by the change in
! the x-tilde coordinates (dx,dy,dz)
!
! INPUT:
!         h : Physical spatial discretization.
!         U : Function on a uniform grid in x-tilde.
!      dx,dy,dz : The change in the x-tilde coordinate value for each point.
!     Mask : The Mask for the uniform grid in x-tilde.
!        shifts : Kind of point interpolation to be used.
!        NOPOINT -> the point should not be filled.
!        NOSHIFT -> the point can be interpolated without shifting
!        positive -> index into shift arrays to get necessary shift
!           value
!        ex : Extents of the arrays
!
! OUTPUT:
!     dest : Function on a uniform grid in x.
!
!--------------------------------------------------------------------------
subroutine make_xgrid(h,U,dx,dy,dz,Mask,shifts,dest,ex,shift_p1,shift_p2,&
             shift_p3,shift_b,num_int_points,params)
use input_params
  implicit none
  integer :: Nshifts
  parameter (Nshifts = 1330)
  integer, dimension(Nshifts),save :: shiftx,shifty,shiftz
  real*8,dimension(:) :: params
  
!--------------------------
  real*8, dimension(1025) :: xx,yy,zz    !temporary storage
  real*8 :: dxt,dyt,dzt,lx,ly,lz
  integer li,lj,lk

!--------------------------------------------------------------------------
!

data shiftx/1,-1,0,0,0,0,1,-1,1,-1,1,-1,0,0,1,-1,0,0,1,-1,1,-1,1,-1,1,-1,2,-2,&
     0,0,0,0,2,-2,1,-1,2,-2,1,-1,2,-2,0,0,1,-1,0,0,2,-2,0,0,1,-1,0,0,&
     2,-2,1,-1,2,-2,1,-1,1,-1,1,-1,2,-2,1,-1,2,-2,1,-1,1,-1,1,-1,2,-2,2,-2,&
     2,-2,0,0,2,-2,0,0,3,-3,0,0,2,-2,2,-2,2,-2,1,-1,2,-2,1,-1,0,2,-2,2,&
     -2,2,-2,1,-1,2,-2,1,-1,0,3,-3,1,-1,3,-3,1,-1,3,-3,0,0,1,-1,0,0,3,-3,&
     0,0,1,-1,0,0,3,-3,1,-1,3,-3,1,-1,1,-1,1,-1,3,-3,1,-1,3,-3,1,-1,1,-1,&
     1,-1,2,-2,2,-2,2,-2,2,-2,3,-3,2,-2,3,-3,2,-2,3,-3,0,0,2,-2,0,0,3,-3,&
     0,0,2,-2,0,0,3,-3,2,-2,3,-3,2,-2,3,-3,1,-1,3,-3,1,-1,2,-2,1,-1,2,-2,&
     1,-1,3,-3,2,-2,3,-3,2,-2,3,-3,1,-1,3,-3,1,-1,2,-2,1,-1,2,-2,1,-1,4,-4,&
     0,0,0,0,4,-4,1,-1,4,-4,1,-1,4,-4,0,0,3,-3,2,-2,3,-3,2,-2,2,-2,2,-2,&
     1,-1,0,0,4,-4,0,0,3,-3,2,-2,3,-3,2,-2,2,-2,2,-2,1,-1,0,0,3,-3,3,-3,&
     4,-4,1,-1,4,-4,1,-1,3,-3,0,0,1,-1,1,-1,4,-4,1,-1,4,-4,1,-1,3,-3,0,0,&
     1,-1,1,-1,3,-3,3,-3,3,-3,1,-1,3,-3,1,-1,3,-3,3,-3,3,-3,1,-1,3,-3,1,-1,&
     4,-4,2,-2,4,-4,2,-2,4,-4,0,0,2,-2,0,0,4,-4,0,0,2,-2,0,0,4,-4,2,-2,&
     4,-4,2,-2,4,-4,1,-1,4,-4,1,-1,2,-2,1,-1,2,-2,1,-1,4,-4,2,-2,4,-4,2,-2,&
     4,-4,1,-1,4,-4,1,-1,2,-2,1,-1,2,-2,1,-1,3,-3,3,-3,3,-3,2,-2,3,-3,2,-2,&
     3,-3,3,-3,3,-3,2,-2,3,-3,2,-2,4,-4,2,-2,4,-4,2,-2,2,-2,2,-2,4,-4,2,-2,&
     4,-4,2,-2,2,-2,2,-2,5,-5,4,-4,3,-3,0,4,-4,3,-3,0,4,-4,0,0,3,-3,0,0,&
     0,4,-4,0,0,3,-3,0,0,0,5,-5,1,-1,5,-5,1,-1,5,-5,4,-4,3,-3,0,4,-4,3,&
     -3,0,4,-4,1,-1,4,-4,1,-1,3,-3,1,-1,3,-3,1,-1,1,-1,0,0,5,-5,4,-4,3,-3,&
     0,4,-4,3,-3,0,4,-4,1,-1,4,-4,1,-1,3,-3,1,-1,3,-3,1,-1,1,-1,0,0,5,-5,&
     1,-1,5,-5,1,-1,3,-3,3,-3,1,-1,1,-1,5,-5,1,-1,5,-5,1,-1,3,-3,3,-3,1,-1,&
     1,-1,5,-5,2,-2,5,-5,2,-2,5,-5,4,-4,3,-3,0,4,-4,3,-3,0,4,-4,2,-2,4,-4,&
     2,-2,3,-3,2,-2,3,-3,2,-2,2,-2,0,0,5,-5,4,-4,3,-3,0,4,-4,3,-3,0,4,-4,&
     2,-2,4,-4,2,-2,3,-3,2,-2,3,-3,2,-2,2,-2,0,0,5,-5,2,-2,5,-5,2,-2,5,-5,&
     1,-1,5,-5,1,-1,2,-2,1,-1,2,-2,1,-1,5,-5,2,-2,5,-5,2,-2,5,-5,1,-1,5,-5,&
     1,-1,2,-2,1,-1,2,-2,1,-1,4,-4,4,-4,4,-4,0,0,4,-4,0,0,4,-4,4,-4,5,-5,&
     2,-2,5,-5,2,-2,4,-4,1,-1,4,-4,1,-1,2,-2,2,-2,4,-4,4,-4,5,-5,2,-2,5,-5,&
     2,-2,4,-4,1,-1,4,-4,1,-1,2,-2,2,-2,5,-5,3,-3,5,-5,3,-3,5,-5,4,-4,3,-3,&
     0,4,-4,3,-3,0,3,-3,3,-3,3,-3,0,0,5,-5,4,-4,3,-3,0,4,-4,3,-3,0,3,-3,&
     3,-3,3,-3,0,0,5,-5,3,-3,5,-5,3,-3,5,-5,1,-1,5,-5,1,-1,3,-3,1,-1,3,-3,&
     1,-1,5,-5,3,-3,5,-5,3,-3,5,-5,1,-1,5,-5,1,-1,3,-3,1,-1,3,-3,1,-1,4,-4,&
     4,-4,4,-4,2,-2,4,-4,2,-2,4,-4,4,-4,4,-4,2,-2,4,-4,2,-2,5,-5,3,-3,5,-5,&
     3,-3,5,-5,2,-2,5,-5,2,-2,3,-3,2,-2,3,-3,2,-2,5,-5,3,-3,5,-5,3,-3,5,-5,&
     2,-2,5,-5,2,-2,3,-3,2,-2,3,-3,2,-2,5,-5,4,-4,5,-5,4,-4,4,-4,4,-4,5,-5,&
     4,-4,3,-3,0,4,-4,3,-3,0,4,-4,0,0,4,-4,4,-4,5,-5,4,-4,3,-3,0,4,-4,3,&
     -3,0,4,-4,0,0,5,-5,4,-4,5,-5,4,-4,5,-5,1,-1,5,-5,1,-1,4,-4,1,-1,4,-4,&
     1,-1,5,-5,4,-4,5,-5,4,-4,5,-5,1,-1,5,-5,1,-1,4,-4,1,-1,4,-4,1,-1,5,-5,&
     3,-3,5,-5,3,-3,3,-3,3,-3,5,-5,3,-3,5,-5,3,-3,3,-3,3,-3,5,-5,4,-4,5,-5,&
     4,-4,5,-5,2,-2,5,-5,2,-2,4,-4,2,-2,4,-4,2,-2,5,-5,4,-4,5,-5,4,-4,5,-5,&
     2,-2,5,-5,2,-2,4,-4,2,-2,4,-4,2,-2,4,-4,4,-4,4,-4,4,-4,5,-5,5,-5,5,-5,&
     4,-4,5,-5,4,-4,5,-5,3,-3,5,-5,3,-3,5,-5,4,-4,3,-3,0,4,-4,3,-3,0,5,-5,&
     4,-4,5,-5,4,-4,5,-5,3,-3,5,-5,3,-3,5,-5,4,-4,3,-3,0,4,-4,3,-3,0,5,-5,&
     5,-5,5,-5,1,-1,5,-5,1,-1,5,-5,5,-5,5,-5,1,-1,5,-5,1,-1,5,-5,5,-5,5,-5,&
     2,-2,5,-5,2,-2,5,-5,5,-5,5,-5,2,-2,5,-5,2,-2,5,-5,4,-4,5,-5,4,-4,4,-4,&
     4,-4,5,-5,4,-4,5,-5,4,-4,4,-4,4,-4,5,-5,5,-5,5,-5,3,-3,5,-5,3,-3,5,-5,&
     5,-5,5,-5,3,-3,5,-5,3,-3,5,-5,5,-5,5,-5,4,-4,5,-5,4,-4,5,-5,5,-5,5,-5,&
     4,-4,5,-5,4,-4,5,-5,5,-5,5,-5,5,-5/
data shifty/0,0,1,-1,0,0,1,1,-1,-1,0,0,1,-1,0,0,1,-1,1,1,-1,-1,1,1,-1,-1,0,0,&
     2,-2,0,0,1,1,2,2,-1,-1,-2,-2,0,0,2,-2,0,0,1,-1,0,0,2,-2,0,0,1,-1,&
     1,1,2,2,-1,-1,-2,-2,1,1,-1,-1,1,1,2,2,-1,-1,-2,-2,1,1,-1,-1,2,2,-2,-2,&
     0,0,2,-2,0,0,2,-2,0,0,3,-3,2,2,-2,-2,1,1,2,2,-1,-1,-2,-2,0,2,2,-2,&
     -2,1,1,2,2,-1,-1,-2,-2,0,1,1,3,3,-1,-1,-3,-3,0,0,3,-3,0,0,1,-1,0,0,&
     3,-3,0,0,1,-1,1,1,3,3,-1,-1,-3,-3,1,1,-1,-1,1,1,3,3,-1,-1,-3,-3,1,1,&
     -1,-1,2,2,-2,-2,2,2,-2,-2,2,2,3,3,-2,-2,-3,-3,0,0,3,-3,0,0,2,-2,0,0,&
     3,-3,0,0,2,-2,2,2,3,3,-2,-2,-3,-3,1,1,3,3,-1,-1,-3,-3,1,1,2,2,-1,-1,&
     -2,-2,2,2,3,3,-2,-2,-3,-3,1,1,3,3,-1,-1,-3,-3,1,1,2,2,-1,-1,-2,-2,0,0,&
     4,-4,0,0,1,1,4,4,-1,-1,-4,-4,0,0,4,-4,2,2,3,3,-2,-2,-3,-3,2,2,-2,-2,&
     0,0,1,-1,0,0,4,-4,2,2,3,3,-2,-2,-3,-3,2,2,-2,-2,0,0,1,-1,3,3,-3,-3,&
     1,1,4,4,-1,-1,-4,-4,0,0,3,-3,1,1,-1,-1,1,1,4,4,-1,-1,-4,-4,0,0,3,-3,&
     1,1,-1,-1,3,3,-3,-3,1,1,3,3,-1,-1,-3,-3,3,3,-3,-3,1,1,3,3,-1,-1,-3,-3,&
     2,2,4,4,-2,-2,-4,-4,0,0,4,-4,0,0,2,-2,0,0,4,-4,0,0,2,-2,2,2,4,4,&
     -2,-2,-4,-4,1,1,4,4,-1,-1,-4,-4,1,1,2,2,-1,-1,-2,-2,2,2,4,4,-2,-2,-4,-4,&
     1,1,4,4,-1,-1,-4,-4,1,1,2,2,-1,-1,-2,-2,3,3,-3,-3,2,2,3,3,-2,-2,-3,-3,&
     3,3,-3,-3,2,2,3,3,-2,-2,-3,-3,2,2,4,4,-2,-2,-4,-4,2,2,-2,-2,2,2,4,4,&
     -2,-2,-4,-4,2,2,-2,-2,0,0,3,3,4,4,5,-3,-3,-4,-4,-5,0,0,4,-4,0,0,3,-3,&
     0,0,0,4,-4,0,0,3,-3,0,1,1,5,5,-1,-1,-5,-5,0,0,3,3,4,4,5,-3,-3,-4,&
     -4,-5,1,1,4,4,-1,-1,-4,-4,1,1,3,3,-1,-1,-3,-3,0,0,1,-1,0,0,3,3,4,4,&
     5,-3,-3,-4,-4,-5,1,1,4,4,-1,-1,-4,-4,1,1,3,3,-1,-1,-3,-3,0,0,1,-1,1,1,&
     5,5,-1,-1,-5,-5,3,3,-3,-3,1,1,-1,-1,1,1,5,5,-1,-1,-5,-5,3,3,-3,-3,1,1,&
     -1,-1,2,2,5,5,-2,-2,-5,-5,0,0,3,3,4,4,5,-3,-3,-4,-4,-5,2,2,4,4,-2,-2,&
     -4,-4,2,2,3,3,-2,-2,-3,-3,0,0,2,-2,0,0,3,3,4,4,5,-3,-3,-4,-4,-5,2,2,&
     4,4,-2,-2,-4,-4,2,2,3,3,-2,-2,-3,-3,0,0,2,-2,2,2,5,5,-2,-2,-5,-5,1,1,&
     5,5,-1,-1,-5,-5,1,1,2,2,-1,-1,-2,-2,2,2,5,5,-2,-2,-5,-5,1,1,5,5,-1,-1,&
     -5,-5,1,1,2,2,-1,-1,-2,-2,4,4,-4,-4,0,0,4,-4,0,0,4,-4,4,4,-4,-4,2,2,&
     5,5,-2,-2,-5,-5,1,1,4,4,-1,-1,-4,-4,2,2,-2,-2,4,4,-4,-4,2,2,5,5,-2,-2,&
     -5,-5,1,1,4,4,-1,-1,-4,-4,2,2,-2,-2,3,3,5,5,-3,-3,-5,-5,0,0,3,3,4,4,&
     5,-3,-3,-4,-4,-5,3,3,-3,-3,0,0,3,-3,0,0,3,3,4,4,5,-3,-3,-4,-4,-5,3,3,&
     -3,-3,0,0,3,-3,3,3,5,5,-3,-3,-5,-5,1,1,5,5,-1,-1,-5,-5,1,1,3,3,-1,-1,&
     -3,-3,3,3,5,5,-3,-3,-5,-5,1,1,5,5,-1,-1,-5,-5,1,1,3,3,-1,-1,-3,-3,4,4,&
     -4,-4,2,2,4,4,-2,-2,-4,-4,4,4,-4,-4,2,2,4,4,-2,-2,-4,-4,3,3,5,5,-3,-3,&
     -5,-5,2,2,5,5,-2,-2,-5,-5,2,2,3,3,-2,-2,-3,-3,3,3,5,5,-3,-3,-5,-5,2,2,&
     5,5,-2,-2,-5,-5,2,2,3,3,-2,-2,-3,-3,4,4,5,5,-4,-4,-5,-5,4,4,-4,-4,0,0,&
     3,3,4,4,5,-3,-3,-4,-4,-5,0,0,4,-4,4,4,-4,-4,0,0,3,3,4,4,5,-3,-3,-4,&
     -4,-5,0,0,4,-4,4,4,5,5,-4,-4,-5,-5,1,1,5,5,-1,-1,-5,-5,1,1,4,4,-1,-1,&
     -4,-4,4,4,5,5,-4,-4,-5,-5,1,1,5,5,-1,-1,-5,-5,1,1,4,4,-1,-1,-4,-4,3,3,&
     5,5,-3,-3,-5,-5,3,3,-3,-3,3,3,5,5,-3,-3,-5,-5,3,3,-3,-3,4,4,5,5,-4,-4,&
     -5,-5,2,2,5,5,-2,-2,-5,-5,2,2,4,4,-2,-2,-4,-4,4,4,5,5,-4,-4,-5,-5,2,2,&
     5,5,-2,-2,-5,-5,2,2,4,4,-2,-2,-4,-4,4,4,-4,-4,4,4,-4,-4,5,5,-5,-5,4,4,&
     5,5,-4,-4,-5,-5,3,3,5,5,-3,-3,-5,-5,0,0,3,3,4,4,5,-3,-3,-4,-4,-5,4,4,&
     5,5,-4,-4,-5,-5,3,3,5,5,-3,-3,-5,-5,0,0,3,3,4,4,5,-3,-3,-4,-4,-5,5,5,&
     -5,-5,1,1,5,5,-1,-1,-5,-5,5,5,-5,-5,1,1,5,5,-1,-1,-5,-5,5,5,-5,-5,2,2,&
     5,5,-2,-2,-5,-5,5,5,-5,-5,2,2,5,5,-2,-2,-5,-5,4,4,5,5,-4,-4,-5,-5,4,4,&
     -4,-4,4,4,5,5,-4,-4,-5,-5,4,4,-4,-4,5,5,-5,-5,3,3,5,5,-3,-3,-5,-5,5,5,&
     -5,-5,3,3,5,5,-3,-3,-5,-5,5,5,-5,-5,4,4,5,5,-4,-4,-5,-5,5,5,-5,-5,4,4,&
     5,5,-4,-4,-5,-5,5,5,-5,-5,5,5,-5,-5/
data shiftz/0,0,0,0,1,-1,0,0,0,0,1,1,1,1,-1,-1,-1,-1,&
     1,1,1,1,-1,-1,-1,-1,0,0,&
     0,0,2,-2,0,0,0,0,0,0,0,0,1,1,1,1,2,2,2,2,-1,-1,-1,-1,-2,-2,-2,-2,&
     1,1,1,1,1,1,1,1,2,2,2,2,-1,-1,-1,-1,-1,-1,-1,-1,-2,-2,-2,-2,0,0,0,0,&
     2,2,2,2,-2,-2,-2,-2,0,0,0,0,1,1,1,1,2,2,2,2,2,2,2,2,3,-1,-1,-1,&
     -1,-2,-2,-2,-2,-2,-2,-2,-2,-3,0,0,0,0,0,0,0,0,1,1,1,1,3,3,3,3,-1,-1,&
     -1,-1,-3,-3,-3,-3,1,1,1,1,1,1,1,1,3,3,3,3,-1,-1,-1,-1,-1,-1,-1,-1,-3,-3,&
     -3,-3,2,2,2,2,-2,-2,-2,-2,0,0,0,0,0,0,0,0,2,2,2,2,3,3,3,3,-2,-2,&
     -2,-2,-3,-3,-3,-3,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3,3,&
     3,3,-1,-1,-1,-1,-1,-1,-1,-1,-2,-2,-2,-2,-2,-2,&
     -2,-2,-3,-3,-3,-3,-3,-3,-3,-3,0,0,&
     0,0,4,-4,0,0,0,0,0,0,0,0,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,&
     4,4,4,4,-1,-1,-1,-1,-2,-2,-2,-2,-2,-2,-2,&
     -2,-3,-3,-3,-3,-4,-4,-4,-4,0,0,0,0,&
     1,1,1,1,1,1,1,1,3,3,3,3,4,4,4,4,-1,-1,-1,-1,-1,-1,-1,-1,-3,-3,-3,-3,&
     -4,-4,-4,-4,1,1,1,1,3,3,3,3,3,3,3,3,-1,-1,-1,-1,-3,-3,-3,-3,-3,-3,-3,-3,&
     0,0,0,0,0,0,0,0,2,2,2,2,4,4,4,4,-2,-2,-2,-2,-4,-4,-4,-4,1,1,1,1,&
     1,1,1,1,2,2,2,2,2,2,2,2,4,4,4,4,4,4,4,4,-1,-1,-1,-1,-1,-1,-1,-1,&
     -2,-2,-2,-2,-2,-2,-2,-2,-4,-4,-4,-4,-4,-4,-4,-4,2,2,2,2,3,3,3,3,3,3,3,3,&
     -2,-2,-2,-2,-3,-3,-3,-3,-3,-3,-3,-3,2,2,2,2,2,2,2,2,4,4,4,4,-2,-2,-2,-2,&
     -2,-2,-2,-2,-4,-4,-4,-4,0,0,0,0,0,0,0,0,0,0,0,0,3,3,3,3,4,4,4,4,&
     5,-3,-3,-3,-3,-4,-4,-4,-4,-5,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,&
     1,1,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,5,5,5,5,-1,-1,-1,-1,-1,-1,&
     -1,-1,-1,-1,-1,-1,-3,-3,-3,-3,-3,-3,-3,-3,-4,&
     -4,-4,-4,-4,-4,-4,-4,-5,-5,-5,-5,1,1,&
     1,1,1,1,1,1,3,3,3,3,5,5,5,5,-1,-1,-1,-1,-1,-1,-1,-1,-3,-3,-3,-3,-5,-5,&
     -5,-5,0,0,0,0,0,0,0,0,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,&
     3,3,4,4,4,4,4,4,4,4,5,5,5,5,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-3,-3,&
     -3,-3,-3,-3,-3,-3,-4,-4,-4,-4,-4,-4,-4,-4,&
     -5,-5,-5,-5,1,1,1,1,1,1,1,1,2,2,&
     2,2,2,2,2,2,5,5,5,5,5,5,5,5,-1,-1,-1,-1,-1,-1,-1,-1,-2,-2,-2,-2,-2,-2,&
     -2,-2,-5,-5,-5,-5,-5,-5,-5,-5,0,0,0,0,4,4,4,4,-4,-4,-4,-4,1,1,1,1,2,2,&
     2,2,2,2,2,2,4,4,4,4,4,4,4,4,5,5,5,5,-1,-1,-1,-1,-2,-2,-2,-2,-2,-2,&
     -2,-2,-4,-4,-4,-4,-4,-4,-4,-4,-5,-5,-5,-5,0,0,0,0,0,0,0,0,3,3,3,3,3,3,&
     3,3,3,3,3,3,4,4,4,4,5,5,5,5,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-4,-4,&
     -4,-4,-5,-5,-5,-5,1,1,1,1,1,1,1,1,3,3,3,3,3,3,3,3,5,5,5,5,5,5,&
     5,5,-1,-1,-1,-1,-1,-1,-1,-1,-3,-3,-3,-3,-3,-3,-3,-3,-5,&
     -5,-5,-5,-5,-5,-5,-5,2,2,&
     2,2,4,4,4,4,4,4,4,4,-2,-2,-2,-2,-4,-4,-4,-4,-4,-4,-4,-4,2,2,2,2,2,2,&
     2,2,3,3,3,3,3,3,3,3,5,5,5,5,5,5,5,5,-2,-2,-2,-2,-2,-2,-2,-2,-3,-3,&
     -3,-3,-3,-3,-3,-3,-5,-5,-5,-5,-5,-5,-5,-5,0,0,0,0,0,0,0,0,3,3,3,3,4,4,&
     4,4,4,4,4,4,4,4,4,4,5,5,5,5,-3,-3,-3,-3,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,&
     -4,-4,-5,-5,-5,-5,1,1,1,1,1,1,1,1,4,4,4,4,4,4,4,4,5,5,5,5,5,5,&
     5,5,-1,-1,-1,-1,-1,-1,-1,-1,-4,-4,-4,-4,-4,-4,-4,&
     -4,-5,-5,-5,-5,-5,-5,-5,-5,3,3,&
     3,3,3,3,3,3,5,5,5,5,-3,-3,-3,-3,-3,-3,-3,-3,-5,-5,-5,-5,2,2,2,2,2,2,&
     2,2,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,-2,-2,-2,-2,-2,-2,-2,-2,-4,-4,&
     -4,-4,-4,-4,-4,-4,-5,-5,-5,-5,-5,-5,-5,&
     -5,4,4,4,4,-4,-4,-4,-4,0,0,0,0,3,3,&
     3,3,3,3,3,3,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,-3,-3,&
     -3,-3,-3,-3,-3,-3,-4,-4,-4,-4,-4,-4,-4,-4,-5,&
     -5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,1,1,&
     1,1,5,5,5,5,5,5,5,5,-1,-1,-1,-1,-5,-5,-5,-5,-5,-5,-5,-5,2,2,2,2,5,5,&
     5,5,5,5,5,5,-2,-2,-2,-2,-5,-5,-5,-5,-5,-5,-5,-5,4,4,4,4,4,4,4,4,5,5,&
     5,5,-4,-4,-4,-4,-4,-4,-4,-4,-5,-5,-5,-5,3,3,3,3,5,5,5,5,5,5,5,5,-3,-3,&
     -3,-3,-5,-5,-5,-5,-5,-5,-5,-5,4,4,4,4,5,5,5,5,5,5,5,5,-4,-4,-4,-4,-5,-5,&
     -5,-5,-5,-5,-5,-5,5,5,5,5,-5,-5,-5,-5/
  !
  real*8,  dimension(3)        :: h
  integer, dimension(3)        :: ex
  real*8,  dimension(:,:,:) :: U,dx,dy,dz,Mask,shifts,dest
  !
  integer :: i,j,k,l,n,maxshift,inc,offset,limit,check,si,sj,sk
  integer, dimension(1) :: maxdir
  integer, dimension(3) :: pnt,ind,shift
  real*8 :: W,Wm1,Wm2,Wm3,Wm4
  real*8, dimension(3) :: hm1,dxs
  real*8, dimension(3) :: X,Xm1,Xm2,Xm3
  real*8, dimension(4,4) :: U1
  real*8, dimension(4) :: U2
  real*8 :: INVALID,ZERO,F1o6,HALF,ONE
  parameter(INVALID = 0.d0,ZERO = 0.d0, F1o6 = 1.d0/6.d0, HALF = 0.5d0, &
       & ONE = 1.d0)
  integer :: NOPOINT,NOSHIFT
  parameter(NOPOINT = -1,NOSHIFT = 0)

  integer        num_int_points
  integer,dimension(:) :: shift_p1,shift_p2,shift_p3,shift_b


  integer :: pn_i, pn_j, pn_k, ind_i, ind_j, ind_k
  real*8 :: dxx, dxy, dxz
  real*8 :: X1, X2, X3
  real*8 :: Xm11, Xm12, Xm13
  real*8 :: XM21, Xm22, Xm23
  real*8 :: Xm31, Xm32, Xm33
  real*8 :: u11, u21, u31, u41
  real*8 :: u12, u22, u32, u42
  real*8 :: u13, u23, u33, u43
  real*8 :: u14, u24, u34, u44
  real*8 :: uu1, uu2, uu3, uu4

  real*8 :: dxepsx,dxepsy,dxepsz,dxeps(3)
   real*8 :: EPSRND,small
   parameter (EPSRND = 1.d06,small=1.0d-10) ! Should prevent most vertex centered rounding
   !            ! anomalies.
  integer i1,i2,i3


  !
  hm1 = ONE/h
  !

  do li = 1 , int(params(P_NX))
    xx(li) = params(P_XMIN) + h(1)*(li-1)
  end do
  do lj = 1 , int(params(P_NY))
    yy(lj) = params(P_YMIN) + h(2)*(lj-1)
  end do
  do lk = 1 , int(params(P_NZ))
    zz(lk) = params(P_ZMIN) + h(3)*(lk-1)
  end do


  do k=1,ex(3)
    do j=1,ex(2)
      do i=1,ex(1)
     pn_i=i
     pn_j=j
     pn_k=k

     dxx=hm1(1) * dx(i,j,k)   
     dxy=hm1(2) * dy(i,j,k)
     dxz=hm1(3) * dz(i,j,k)
     dxepsx= epsrnd*spacing(dxx)
     dxepsy= epsrnd*spacing(dxy)
     dxepsz= epsrnd*spacing(dxz)


     ind_i = max(floor(dxx+dxepsx) + pn_i -1,1)
     ind_j = max(floor(dxy+dxepsy) + pn_j -1,1)
     ind_k = max(floor(dxz+dxepsz) + pn_k -1,1)
!     dest(i,j,k) = u(ind_i,ind_j,ind_k)
!	  goto 800

!     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_i + 3 > ex(1) ) ind_i = ex(1) - 3
     if(ind_j + 3 > ex(2) ) ind_j = ex(2) - 3
     if(ind_k + 3 > ex(3) ) ind_k = ex(3) - 3

     X1 = dxx + pn_i - ind_i  ! Get normalized distance from starting point.
     X2 = dxy + pn_j - ind_j
     X3 = dxz + pn_k - ind_k
     Xm11 = X1 - ONE
     Xm12 = X2 - ONE
     Xm13 = X3 - ONE
     Xm21 = Xm11 - ONE
     Xm22 = Xm12 - ONE
     Xm23 = Xm13 - ONE
     Xm31 = Xm21 - ONE
     Xm32 = Xm22 - ONE
     Xm33 = Xm23 - ONE
     !
      u11 = F1o6*Xm11*Xm21*(     X1*U(ind_i+3,ind_j  ,ind_k  )    &
     &          -Xm31*U(ind_i  ,ind_j  ,ind_k  ) )    &
     & + HALF*X1*Xm31*(   Xm21*U(ind_i+1,ind_j  ,ind_k  )    &
     &          -Xm11*U(ind_i+2,ind_j  ,ind_k  ) )

      u21 = F1o6*Xm11*Xm21*(     X1*U(ind_i+3,ind_j+1,ind_k  )    &
     &          -Xm31*U(ind_i  ,ind_j+1,ind_k  ) )    &
     & + HALF*X1*Xm31*(   Xm21*U(ind_i+1,ind_j+1,ind_k  )    &
     &          -Xm11*U(ind_i+2,ind_j+1,ind_k  ) )

      u31 = F1o6*Xm11*Xm21*(     X1*U(ind_i+3,ind_j+2,ind_k  )    &
     &          -Xm31*U(ind_i  ,ind_j+2,ind_k  ) )    &
     & + HALF*X1*Xm31*(   Xm21*U(ind_i+1,ind_j+2,ind_k  )    &
     &          -Xm11*U(ind_i+2,ind_j+2,ind_k  ) )

      u41 = F1o6*Xm11*Xm21*(     X1*U(ind_i+3,ind_j+3,ind_k  )    &
     &          -Xm31*U(ind_i  ,ind_j+3,ind_k  ) )    &
     & + HALF*X1*Xm31*(   Xm21*U(ind_i+1,ind_j+3,ind_k  )    &
     &          -Xm11*U(ind_i+2,ind_j+3,ind_k  ) )

      u12 = F1o6*Xm11*Xm21*(     X1*U(ind_i+3,ind_j  ,ind_k+1)    &
     &          -Xm31*U(ind_i  ,ind_j  ,ind_k+1) )    &
     & + HALF*X1*Xm31*(   Xm21*U(ind_i+1,ind_j  ,ind_k+1)    &
     &          -Xm11*U(ind_i+2,ind_j  ,ind_k+1) )

      u22 = F1o6*Xm11*Xm21*(     X1*U(ind_i+3,ind_j+1,ind_k+1)    &
     &          -Xm31*U(ind_i  ,ind_j+1,ind_k+1) )    &
     & + HALF*X1*Xm31*(   Xm21*U(ind_i+1,ind_j+1,ind_k+1)    &
     &          -Xm11*U(ind_i+2,ind_j+1,ind_k+1) )

      u32 = F1o6*Xm11*Xm21*(     X1*U(ind_i+3,ind_j+2,ind_k+1)    &
     &          -Xm31*U(ind_i  ,ind_j+2,ind_k+1) )    &
     & + HALF*X1*Xm31*(   Xm21*U(ind_i+1,ind_j+2,ind_k+1)    &
     &          -Xm11*U(ind_i+2,ind_j+2,ind_k+1) )

      u42 = F1o6*Xm11*Xm21*(     X1*U(ind_i+3,ind_j+3,ind_k+1)    &
     &          -Xm31*U(ind_i  ,ind_j+3,ind_k+1) )    &
     & + HALF*X1*Xm31*(   Xm21*U(ind_i+1,ind_j+3,ind_k+1)    &
     &          -Xm11*U(ind_i+2,ind_j+3,ind_k+1) )

      u13 = F1o6*Xm11*Xm21*(     X1*U(ind_i+3,ind_j  ,ind_k+2)    &
     &          -Xm31*U(ind_i  ,ind_j  ,ind_k+2) )    &
     & + HALF*X1*Xm31*(   Xm21*U(ind_i+1,ind_j  ,ind_k+2)    &
     &          -Xm11*U(ind_i+2,ind_j  ,ind_k+2) )

      u23 = F1o6*Xm11*Xm21*(     X1*U(ind_i+3,ind_j+1,ind_k+2)    &
     &          -Xm31*U(ind_i  ,ind_j+1,ind_k+2) )    &
     & + HALF*X1*Xm31*(   Xm21*U(ind_i+1,ind_j+1,ind_k+2)    &
     &          -Xm11*U(ind_i+2,ind_j+1,ind_k+2) )

      u33 = F1o6*Xm11*Xm21*(     X1*U(ind_i+3,ind_j+2,ind_k+2)    &
     &          -Xm31*U(ind_i  ,ind_j+2,ind_k+2) )    &
     & + HALF*X1*Xm31*(   Xm21*U(ind_i+1,ind_j+2,ind_k+2)    &
     &          -Xm11*U(ind_i+2,ind_j+2,ind_k+2) )

      u43 = F1o6*Xm11*Xm21*(     X1*U(ind_i+3,ind_j+3,ind_k+2)    &
     &          -Xm31*U(ind_i  ,ind_j+3,ind_k+2) )    &
     & + HALF*X1*Xm31*(   Xm21*U(ind_i+1,ind_j+3,ind_k+2)    &
     &          -Xm11*U(ind_i+2,ind_j+3,ind_k+2) )

      u14 = F1o6*Xm11*Xm21*(     X1*U(ind_i+3,ind_j  ,ind_k+3)    &
     &          -Xm31*U(ind_i  ,ind_j  ,ind_k+3) )    &
     & + HALF*X1*Xm31*(   Xm21*U(ind_i+1,ind_j  ,ind_k+3)    &
     &          -Xm11*U(ind_i+2,ind_j  ,ind_k+3) )

      u24 = F1o6*Xm11*Xm21*(     X1*U(ind_i+3,ind_j+1,ind_k+3)    &
     &          -Xm31*U(ind_i  ,ind_j+1,ind_k+3) )    &
     & + HALF*X1*Xm31*(   Xm21*U(ind_i+1,ind_j+1,ind_k+3)    &
     &          -Xm11*U(ind_i+2,ind_j+1,ind_k+3) )

      u34 = F1o6*Xm11*Xm21*(     X1*U(ind_i+3,ind_j+2,ind_k+3)    &
     &          -Xm31*U(ind_i  ,ind_j+2,ind_k+3) )    &
     & + HALF*X1*Xm31*(   Xm21*U(ind_i+1,ind_j+2,ind_k+3)    &
     &          -Xm11*U(ind_i+2,ind_j+2,ind_k+3) )

      u44 = F1o6*Xm11*Xm21*(     X1*U(ind_i+3,ind_j+3,ind_k+3)    &
     &          -Xm31*U(ind_i  ,ind_j+3,ind_k+3) )    &
     & + HALF*X1*Xm31*(   Xm21*U(ind_i+1,ind_j+3,ind_k+3)    &
     &          -Xm11*U(ind_i+2,ind_j+3,ind_k+3) )

!  U1 =   F1o6*Xm11*Xm21*(                  &
!       &      X1*U(ind_i+3,ind_j:ind_j+3,ind_k:ind_k+3)  &
!       &       - Xm31*U(ind_i  ,ind_j:ind_j+3,ind_k:ind_k+3) )&
!       & + HALF*X1*Xm31*(                  &
!       &         Xm21*U(ind_i+1,ind_j:ind_j+3,ind_k:ind_k+3)  &
!       &       - Xm11*U(ind_i+2,ind_j:ind_j+3,ind_k:ind_k+3) )

!  U2 =   F1o6*Xm12*Xm22*(  X2*U1(4,:) - Xm32*U1(1,:))        &
!       & + HALF*X2*Xm32*(Xm22*U1(2,:) - Xm12*U1(3,:))

      uu1 = F1o6*Xm12*Xm22*(   X2*u41 - Xm32*u11 )    &
     & + HALF*X2*Xm32*( Xm22*u21 - Xm12*u31 )
      uu2 = F1o6*Xm12*Xm22*(   X2*u42 - Xm32*u12 )    &
     & + HALF*X2*Xm32*( Xm22*u22 - Xm12*u32 )
      uu3 = F1o6*Xm12*Xm22*(   X2*u43 - Xm32*u13 )    &
     & + HALF*X2*Xm32*( Xm22*u23 - Xm12*u33 )
      uu4 = F1o6*Xm12*Xm22*(   X2*u44 - Xm32*u14 )    &
     & + HALF*X2*Xm32*( Xm22*u24 - Xm12*u34 )

!  dest(i,j,k) = F1o6*Xm13*Xm23*(  X3*U2(4) - Xm33*U2(1))     &
!       &      + HALF*X3*Xm33*(Xm23*U2(2) - Xm13*U2(3))

      dest(i,j,k) = F1o6*Xm13*Xm23*(   X3*uu4 - Xm33*uu1 )    &
         & + HALF*X3*Xm33*( Xm23*uu2 - Xm13*uu3 )
!        end if
800   continue
      end do
    end do
  end do

  do l = 1, num_int_points
   i = shift_p1(l)
   j = shift_p2(l)
   k = shift_p3(l)
   n = int(shifts(i,j,k))
     pnt(1) = i
     pnt(2) = j
     pnt(3) = k
     dxs = hm1* (/ dx(i,j,k), dy(i,j,k), dz(i,j,k) /)
     dxeps = EPSRND*SPACING(dxs)
       select case(shift_b(l))
     case(12)
         shift(1) = shiftx(n)
         shift(2) = shifty(n)
         shift(3) = shiftz(n)
     case(13)
         shift(1) = shiftx(n)
         shift(2) = shiftz(n)
         shift(3) = shifty(n)
     case(21)
         shift(1) = shifty(n)
         shift(2) = shiftx(n)
         shift(3) = shiftz(n)
     case(23)
         shift(1) = shiftz(n)
         shift(2) = shiftx(n)
         shift(3) = shifty(n)
     case(31)
         shift(1) = shifty(n)
         shift(2) = shiftz(n)
         shift(3) = shiftx(n)
     case(32)
         shift(1) = shiftz(n)
         shift(2) = shifty(n)
         shift(3) = shiftx(n)
       end select

!make short vectors into scalars
     ind = FLOOR(dxs+dxeps) + pnt - 1 + shift
     X = dxs + pnt - ind  ! Get normalized distance from starting point.
     Xm1 = X - ONE
     Xm2 = Xm1 - ONE
     Xm3 = Xm2 - ONE
     !
     ! Prepare for optimizing largest shift:
     ! In the direction of the largest shift, we will "slide" the 1-D
     ! lines to them minimum value which is still valid

!     maxdir = MAXLOC(max(x-3.0d0,-x)) ! index of maximum shift
!-------------------------------------------------------

       lx = xx(i) + dx(i,j,k)
       ly = yy(j) + dy(i,j,k)
       lz = zz(k) + dz(i,j,k)
       dxt  = abs(lx-params(P_HOLE1_X))
       dyt  = abs(ly-params(P_HOLE1_Y))
       dzt  = abs(lz-params(P_HOLE1_Z))


       if (abs(dxt-dyt).le.small) dyt = dxt
       if (abs(dxt-dzt).le.small) dzt = dxt
       if (abs(dzt-dyt).le.small) dzt = dyt


       if (dxt.lt.dyt.and.dxt.lt.dzt) then
          maxdir = 1
       else if (dyt.lt.dzt.and.dyt.lt.dxt)  then
          maxdir = 2
       else if (dzt.lt.dyt.and.dzt.lt.dxt) then
          maxdir = 3
		 else if (dxt.eq.dyt.and.dxt.eq.dzt) then
			 maxdir = 1
       else if (dyt.eq.dzt) then   !two terms are equal
         if (ly.lt.lz) then
            maxdir = 2
         else   if (lz.lt.ly) then
            maxdir = 3
         else
            maxdir = 2
         end if
       else if (dyt.eq.dxt) then   !two terms are equal
         if (ly.lt.lx) then
            maxdir = 2
         else   if (lx.lt.ly) then
            maxdir = 1
			else
				maxdir = 1
         end if
       else if (dxt.eq.dzt) then   !two terms are equal
         if (lx.lt.lz) then
            maxdir = 1
         else   if (lz.lt.lx) then
            maxdir = 3
			else
				maxdir = 3
         end if
       end if

!-------------------------------------------------------

     maxshift = shift(maxdir(1)) ! value of maximum shift
     !
     if(maxshift > 0) then  ! If we shifted up, then we will test the low
       inc = -1         ! end as we slide back down.
       offset = 0

       limit = MIN(ind(maxdir(1))-1,maxshift)
     else          ! If we shifted down, then we will test the
       inc = 1     ! high end as we slide back up.
       offset = 3
       limit = MIN(ex(maxdir(1))-ind(maxdir(1))-3,-maxshift)
     end if
    select case (maxdir(1)) ! Minimize biggest shift for line interp
    case (1) ! Minimize shift for each X line interp/extrap
!unroll sk and sj
       do sk=ind(3),ind(3)+3
         do sj=ind(2),ind(2)+3
        W = X(1)
        check = ind(1) + offset + inc
!line shift

if (n.lt.999) then
        do si=1,limit ! Unshift until the data is invalid
      if(Mask(check,sj,sk) == INVALID) then
         exit
      else
         W = W - inc
         check = check + inc
      end if
        end do
end if



        si = check - offset - inc ! then interpolate from this index 
        Wm1 = W - ONE
        Wm2 = Wm1 - ONE
        Wm3 = Wm2 - ONE
        U1(sj-ind(2)+1,sk-ind(3)+1) =               &
       F1o6*Wm1*Wm2*(W  *U(si+3,sj,sk) - Wm3*U(si  ,sj,sk))      &
     + HALF*W  *Wm3*(Wm2*U(si+1,sj,sk) - Wm1*U(si+2,sj,sk))
         end do
      end do

      U2 =   F1o6*Xm1(2)*Xm2(2)*(  X(2)*U1(4,:) - Xm3(2)*U1(1,:))       &
        + HALF*X(2)*Xm3(2)*(Xm2(2)*U1(2,:) - Xm1(2)*U1(3,:))

      dest(i,j,k) = F1o6*Xm1(3)*Xm2(3)*(  X(3)*U2(4) - Xm3(3)*U2(1))    &
         + HALF*X(3)*Xm3(3)*(Xm2(3)*U2(2) - Xm1(3)*U2(3))


    case (2) ! Minimize shift for each Y line interp/extrap
       do sk=ind(3),ind(3)+3
         do si=ind(1),ind(1)+3
        W = X(2)
        check = ind(2) + offset + inc
! line shift
if (n.lt.999) then
        do sj=1,limit ! Unshift until the data is invalid
      if(Mask(si,check,sk) == INVALID) then
         exit
      else
         W = W - inc
         check = check + inc
      end if
        end do
end if
        sj = check - offset - inc ! then interpolate from this index 
        Wm1 = W - ONE
        Wm2 = Wm1 - ONE
        Wm3 = Wm2 - ONE
        U1(si-ind(1)+1,sk-ind(3)+1) =               &
       F1o6*Wm1*Wm2*(W  *U(si,sj+3,sk) - Wm3*U(si,sj  ,sk)) &
     + HALF*W  *Wm3*(Wm2*U(si,sj+1,sk) - Wm1*U(si,sj+2,sk))
         end do
       end do
       U2 =   F1o6*Xm1(1)*Xm2(1)*(  X(1)*U1(4,:) - Xm3(1)*U1(1,:))       &
        + HALF*X(1)*Xm3(1)*(Xm2(1)*U1(2,:) - Xm1(1)*U1(3,:))

       dest(i,j,k) = F1o6*Xm1(3)*Xm2(3)*(  X(3)*U2(4) - Xm3(3)*U2(1))    &
         + HALF*X(3)*Xm3(3)*(Xm2(3)*U2(2) - Xm1(3)*U2(3))

    case (3) ! Minimize shift for each Z line interp/extrap
       do sj=ind(2),ind(2)+3
     do si=ind(1),ind(1)+3
         W = X(3)
         check = ind(3) + offset + inc
!line shifts
if (n.lt.999) then
         do sk=1,limit ! Unshift until the data is invalid
       if(Mask(si,sj,check) == INVALID) then
          exit
       else
          W = W - inc
          check = check + inc
       end if
         end do
end if
         sk = check - offset - inc ! then interpolate from this index 
         Wm1 = W - ONE
         Wm2 = Wm1 - ONE
         Wm3 = Wm2 - ONE
         U1(si-ind(1)+1,sj-ind(2)+1) =               &
        F1o6*Wm1*Wm2*(W  *U(si,sj,sk+3) - Wm3*U(si,sj,sk  )) &
      + HALF*W  *Wm3*(Wm2*U(si,sj,sk+1) - Wm1*U(si,sj,sk+2))
     end do
       end do
       U2 =   F1o6*Xm1(1)*Xm2(1)*(  X(1)*U1(4,:) - Xm3(1)*U1(1,:))       &
     + HALF*X(1)*Xm3(1)*(Xm2(1)*U1(2,:) - Xm1(1)*U1(3,:))

       dest(i,j,k) = F1o6*Xm1(2)*Xm2(2)*(  X(2)*U2(4) - Xm3(2)*U2(1))    &
         + HALF*X(2)*Xm3(2)*(Xm2(2)*U2(2) - Xm1(2)*U2(3))


   end select

! dump this out .............................
!		write(88,222) sqrt(lx*lx+ly*ly+lz*lz),dest(i,j,k),u(i,j,k),i,j,k,shift(1),shift(2),shift(3),maxdir

end do

222 format(3e14.4,7i5)
  !
end subroutine make_xgrid
!-----------------------------------------------------------------------
