!/*@@
!
! $Id: intp_levels.f90,v 1.3 1997/08/26 07:42:51 haupt Exp $
!
!@@*/
!----------------------------------------------------------------------
! Given points in a list, this routine will interpolate to determine the
! values
! We must interface this with Greg's interpolators
! this routine determines the shifts and the bias to be used later
!-----------------------------------------------------------------------
subroutine pre_interp_levels(mask,params,si,sj,sk,np,shift_t,shift_b)
use input_params
implicit none
real*8,dimension(:,:,:) :: mask
real*8,dimension(:)     :: params
integer                 :: np
integer,dimension(:)    :: si,sj,sk
integer,dimension(:)    :: shift_t,shift_b 
integer,parameter       :: Nshifts=729
integer, dimension(Nshifts),save :: shiftx,shifty,shiftz
integer                 :: i,j,k,bias
integer                 :: l,li,lj,lk,test,n
real*8                  :: lx,ly,lz,dx1,dy1,dz1
integer                 :: shifti,shiftj,shiftk,nx,ny,nz,shp(3)
real*8                  :: h
real*8, dimension(1025) :: x,y,z    !temporary storage



data shiftx /0,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,4,-4,3,-3,4,-4,3,-3,4,-4,&
     0,0,3,-3,0,0,4,-4,0,0,3,-3,0,0,4,-4,3,-3,4,-4,3,-3,4,-4,1,-1,&
     4,-4,1,-1,3,-3,1,-1,3,-3,1,-1,4,-4,3,-3,4,-4,3,-3,4,-4,1,-1,4,-4,&
     1,-1,3,-3,1,-1,3,-3,1,-1,3,-3,3,-3,3,-3,3,-3,4,-4,3,-3,4,-4,3,-3,&
     4,-4,2,-2,4,-4,2,-2,3,-3,2,-2,3,-3,2,-2,4,-4,3,-3,4,-4,3,-3,4,-4,&
     2,-2,4,-4,2,-2,3,-3,2,-2,3,-3,2,-2,4,-4,4,-4,4,-4,0,0,4,-4,0,0,&
     4,-4,4,-4,4,-4,1,-1,4,-4,1,-1,4,-4,4,-4,4,-4,1,-1,4,-4,1,-1,4,-4,&
     3,-3,4,-4,3,-3,3,-3,3,-3,4,-4,3,-3,4,-4,3,-3,3,-3,3,-3,4,-4,4,-4,&
     4,-4,2,-2,4,-4,2,-2,4,-4,4,-4,4,-4,2,-2,4,-4,2,-2,4,-4,4,-4,4,-4,&
     3,-3,4,-4,3,-3,4,-4,4,-4,4,-4,3,-3,4,-4,3,-3,4,-4,4,-4,4,-4,4,-4/
data shifty /0,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,3,3,4,4,-3,-3,-4,-4,0,0,&
         4,-4,0,0,3,-3,0,0,4,-4,0,0,3,-3,3,3,4,4,-3,-3,-4,-4,1,1,4,4,&
         -1,-1,-4,-4,1,1,3,3,-1,-1,-3,-3,3,3,4,4,-3,-3,-4,-4,1,1,4,4,-1,-1,&
         -4,-4,1,1,3,3,-1,-1,-3,-3,3,3,-3,-3,3,3,-3,-3,3,3,4,4,-3,-3,-4,-4,&
         2,2,4,4,-2,-2,-4,-4,2,2,3,3,-2,-2,-3,-3,3,3,4,4,-3,-3,-4,-4,2,2,&
         4,4,-2,-2,-4,-4,2,2,3,3,-2,-2,-3,-3,4,4,-4,-4,0,0,4,-4,0,0,4,-4,&
         4,4,-4,-4,1,1,4,4,-1,-1,-4,-4,4,4,-4,-4,1,1,4,4,-1,-1,-4,-4,3,3,&
         4,4,-3,-3,-4,-4,3,3,-3,-3,3,3,4,4,-3,-3,-4,-4,3,3,-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,4,4,-4,-4,3,3,&
         4,4,-3,-3,-4,-4,4,4,-4,-4,3,3,4,4,-3,-3,-4,-4,4,4,-4,-4,4,4,-4,-4/
data shiftz /0,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,3,3,&
     3,3,4,4,4,4,-3,-3,-3,-3,-4,-4,-4,-4,1,1,1,1,1,1,1,1,3,3,3,3,&
         3,3,3,3,4,4,4,4,4,4,4,4,-1,-1,-1,-1,-1,-1,-1,-1,-3,-3,-3,-3,-3,-3,&
         -3,-3,-4,-4,-4,-4,-4,-4,-4,-4,3,3,3,3,-3,-3,-3,-3,2,2,2,2,2,2,2,2,&
    3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,-2,-2,-2,-2,-2,-2,-2,-2,-3,-3,&
    -3,-3,-3,-3,-3,-3,-4,-4,-4,-4,-4,-4,-4,-4,0,0,0,0,4,4,4,4,-4,-4,-4,-4,&
      1,1,1,1,4,4,4,4,4,4,4,4,-1,-1,-1,-1,-4,-4,-4,-4,-4,-4,-4,-4,3,3,&
      3,3,3,3,3,3,4,4,4,4,-3,-3,-3,-3,-3,-3,-3,-3,-4,-4,-4,-4,2,2,2,2,&
      4,4,4,4,4,4,4,4,-2,-2,-2,-2,-4,-4,-4,-4,-4,-4,-4,-4,3,3,3,3,4,4,&
      4,4,4,4,4,4,-3,-3,-3,-3,-4,-4,-4,-4,-4,-4,-4,-4,4,4,4,4,-4,-4,-4,-4/

  h = params(P_GRIDSPACING)

  shp = shape(mask)
  nx = shp(1) 
  ny = shp(2) 
  nz = shp(3) 


  do li = 1 , nx 
    x(li) = params(P_XMIN) + h*(li-1+int(params(P_LOCAL_LOWER_X)))
  end do
  do lj = 1 , ny 
    y(lj) = params(P_YMIN) + h*(lj-1+int(params(P_LOCAL_LOWER_Y)))
  end do
  do lk = 1 , nz 
    z(lk) = params(P_ZMIN) + h*(lk-1+int(params(P_LOCAL_LOWER_Z)))
  end do


pointloop:  do l = 1 , np

  i = si(l)
  j = sj(l)
  k = sk(l)
  lx = x(i) - params(P_HOLE1_X)
  ly = y(j) - params(P_HOLE1_Y)
  lz = z(k) - params(P_HOLE1_Z)

  dx1 = abs(lx)
  dy1 = abs(ly)
  dz1 = abs(lz)

  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)
            shifti = shiftx(n)
            shiftj = shifty(n)
            shiftk = shiftz(n)
         case(13)
            shifti = shiftx(n)
            shiftj = shiftz(n)
            shiftk = shifty(n)
         case(21)
            shifti = shifty(n)
            shiftj = shiftx(n)
            shiftk = shiftz(n)
         case(23)
            shifti = shiftz(n)
            shiftj = shiftx(n)
            shiftk = shifty(n)
         case(31)
            shifti = shifty(n)
            shiftj = shiftz(n)
            shiftk = shiftx(n)
         case(32)
            shifti = shiftz(n)
            shiftj = shifty(n)
            shiftk = shiftx(n)
      end select
! change this from i:i+3  , to i-1:i+2
      test = int( sum(mask(i+shifti-1:i+shifti+2,j+shiftj-1:j+shiftj+2, &
                           k+shiftk-1:k+shiftk+2)) )
      if (test.eq.64) then
         shift_t(l) = n
         shift_b(l) = bias 
         cycle pointloop
      end if
  end do

end do pointloop

return
end

!----------------------------------------------------------------------
! Given points in a list, this routine will interpolate to determine the
! values
! We must interface this with Greg's interpolators
!-----------------------------------------------------------------------
subroutine interp_levels(u,mask,ssi,ssj,ssk,np,shift_t,shift_b,dest,ii,jj)
use input_params
implicit none
real*8,dimension(:,:,:) :: mask
integer                 :: np
integer,dimension(:)    :: ssi,ssj,ssk
integer,dimension(:)    :: shift_t,shift_b 
integer,parameter       :: Nshifts=729
integer, dimension(Nshifts),save :: shiftx,shifty,shiftz
integer                 :: i,j,k,bias
real*8                  :: lx,ly,lz,dx1,dy1,dz1,shifti,shiftj,shiftk,dx,dy,dz
real*8                  :: dxt,dyt,dzt
real*8                  :: h
real*8, dimension(1025) :: xx,yy,zz    !temporary storage
integer                 :: si,sj,sk,offset
real*8,dimension(:,:,:) :: u,dest
integer                 :: ii,jj
integer,dimension(3)    :: pnt,shift,ind
real*8,dimension(3)     :: dxs,x,xm1,xm2,xm3
real*8                  :: w,wm1,wm2,wm3
real*8,parameter        :: one = 1.0d0, half=0.5d0,F1o6=1.d0/6.0d0,zero=0.0d0
real*8,parameter        :: small = 1.0e-10, invalid = 0.0d0
real*8, dimension(4)    :: u2 
real*8, dimension(4,4)  :: U1
integer,dimension(3)     :: ex

integer                 :: li,lj,lk,limit,maxshift,check,l,n,inc
integer,dimension(1)    :: maxdir


data shiftx /0,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,4,-4,3,-3,4,-4,3,-3,4,-4,&
     0,0,3,-3,0,0,4,-4,0,0,3,-3,0,0,4,-4,3,-3,4,-4,3,-3,4,-4,1,-1,&
     4,-4,1,-1,3,-3,1,-1,3,-3,1,-1,4,-4,3,-3,4,-4,3,-3,4,-4,1,-1,4,-4,&
     1,-1,3,-3,1,-1,3,-3,1,-1,3,-3,3,-3,3,-3,3,-3,4,-4,3,-3,4,-4,3,-3,&
     4,-4,2,-2,4,-4,2,-2,3,-3,2,-2,3,-3,2,-2,4,-4,3,-3,4,-4,3,-3,4,-4,&
     2,-2,4,-4,2,-2,3,-3,2,-2,3,-3,2,-2,4,-4,4,-4,4,-4,0,0,4,-4,0,0,&
     4,-4,4,-4,4,-4,1,-1,4,-4,1,-1,4,-4,4,-4,4,-4,1,-1,4,-4,1,-1,4,-4,&
     3,-3,4,-4,3,-3,3,-3,3,-3,4,-4,3,-3,4,-4,3,-3,3,-3,3,-3,4,-4,4,-4,&
     4,-4,2,-2,4,-4,2,-2,4,-4,4,-4,4,-4,2,-2,4,-4,2,-2,4,-4,4,-4,4,-4,&
     3,-3,4,-4,3,-3,4,-4,4,-4,4,-4,3,-3,4,-4,3,-3,4,-4,4,-4,4,-4,4,-4/
data shifty /0,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,3,3,4,4,-3,-3,-4,-4,0,0,&
         4,-4,0,0,3,-3,0,0,4,-4,0,0,3,-3,3,3,4,4,-3,-3,-4,-4,1,1,4,4,&
         -1,-1,-4,-4,1,1,3,3,-1,-1,-3,-3,3,3,4,4,-3,-3,-4,-4,1,1,4,4,-1,-1,&
         -4,-4,1,1,3,3,-1,-1,-3,-3,3,3,-3,-3,3,3,-3,-3,3,3,4,4,-3,-3,-4,-4,&
         2,2,4,4,-2,-2,-4,-4,2,2,3,3,-2,-2,-3,-3,3,3,4,4,-3,-3,-4,-4,2,2,&
         4,4,-2,-2,-4,-4,2,2,3,3,-2,-2,-3,-3,4,4,-4,-4,0,0,4,-4,0,0,4,-4,&
         4,4,-4,-4,1,1,4,4,-1,-1,-4,-4,4,4,-4,-4,1,1,4,4,-1,-1,-4,-4,3,3,&
         4,4,-3,-3,-4,-4,3,3,-3,-3,3,3,4,4,-3,-3,-4,-4,3,3,-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,4,4,-4,-4,3,3,&
         4,4,-3,-3,-4,-4,4,4,-4,-4,3,3,4,4,-3,-3,-4,-4,4,4,-4,-4,4,4,-4,-4/
data shiftz /0,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,3,3,&
     3,3,4,4,4,4,-3,-3,-3,-3,-4,-4,-4,-4,1,1,1,1,1,1,1,1,3,3,3,3,&
         3,3,3,3,4,4,4,4,4,4,4,4,-1,-1,-1,-1,-1,-1,-1,-1,-3,-3,-3,-3,-3,-3,&
         -3,-3,-4,-4,-4,-4,-4,-4,-4,-4,3,3,3,3,-3,-3,-3,-3,2,2,2,2,2,2,2,2,&
    3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,-2,-2,-2,-2,-2,-2,-2,-2,-3,-3,&
    -3,-3,-3,-3,-3,-3,-4,-4,-4,-4,-4,-4,-4,-4,0,0,0,0,4,4,4,4,-4,-4,-4,-4,&
      1,1,1,1,4,4,4,4,4,4,4,4,-1,-1,-1,-1,-4,-4,-4,-4,-4,-4,-4,-4,3,3,&
      3,3,3,3,3,3,4,4,4,4,-3,-3,-3,-3,-3,-3,-3,-3,-4,-4,-4,-4,2,2,2,2,&
      4,4,4,4,4,4,4,4,-2,-2,-2,-2,-4,-4,-4,-4,-4,-4,-4,-4,3,3,3,3,4,4,&
      4,4,4,4,4,4,-3,-3,-3,-3,-4,-4,-4,-4,-4,-4,-4,-4,4,4,4,4,-4,-4,-4,-4/

	ex = shape(u)

   do l = 1, np
      i = ssi(l)
      j = ssj(l)
      k = ssk(l)
      n = shift_t(l) 
      pnt(1) = i
      pnt(2) = j
      pnt(3) = k
      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
      dxs(1) = 1.0d0*shift(1)
      dxs(2) = 1.0d0*shift(2)
      dxs(3) = 1.0d0*shift(3)

      X = - dxs +1
      ind(1) = i + int(shift(1)) -1
      ind(2) = j + int(shift(2)) -1
      ind(3) = k + int(shift(3)) -1
     maxdir = MAXLOC(max(x-3.0d0,-x)) ! index of maximum shift

!22 format(1P,6e10.2,10i4)
      Xm1 = X - ONE
      Xm2 = Xm1 - ONE
      Xm3 = Xm2 - ONE


      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
         do sk=ind(3),ind(3)+3
         do sj=ind(2),ind(2)+3
            W = X(1)
            check = ind(1) + offset + inc
            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))
!			if (ii.eq.1.and.jj.eq.1)write(39,34) i,j,k,U(si+3,sj,sk),U(si,sj,sk), &
!         U(si+1,sj,sk),U(si+2,sj,sk),&
!         mask(si+3,sj,sk),mask(si,sj,sk), &
!         mask(si+1,sj,sk),mask(si+2,sj,sk),si,sj,sk,1
      end do
      end do

34 format(3i4,1P,4e10.2,4e9.1,4i3)
      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
          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))
!			if (ii.eq.1.and.jj.eq.1)write(39,34) i,j,k,U(si,sj+3,sk),U(si,sj,sk), &
!                             U(si,sj+1,sk),U(si,sj+2,sk),  &
!         mask(si,sj+3,sk),mask(si,sj,sk), &
!         mask(si,sj+1,sk),mask(si,sj+2,sk),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(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) 
       do sj=ind(2),ind(2)+3
       do si=ind(1),ind(1)+3
          W = X(3)
          check = ind(3) + offset + inc
          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))
!			if (ii.eq.1.and.jj.eq.1)write(39,34) i,j,k,U(si,sj,sk+3),U(si,sj,sk), &
!                            U(si,sj,sk+1),U(si,sj,sk+2), &
!         mask(si,sj,sk+3 ),mask(si,sj,sk), &
!         mask(si,sj,sk+1),mask(si,sj,sk+2),si,sj,sk,3
       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
!	 if (ii.eq.1.and.jj.eq.1) write(22,'(3i5,1P,e14.4)') i,j,k,dest(i,j,k)
end do

!if (ii.eq.1.and.jj.eq.1) write(39,*) '------------'
!if (ii.eq.1.and.jj.eq.1) write(22,*) '------------'
end subroutine 
