subroutine causual_cn_pre(x_shift_,srcMask_, dstMask_,shifts_,   &
                          params,num_int_points, &
                          shift_p1,shift_p2,shift_p3, shift_b, &
                          gm1_,km1_,gp1_,kp1_,m,nl)

  use input_params
!  use parameters
  use tensors
  implicit none

  real*8, dimension(100)                         :: params
  type(GridFunction), dimension(:)               :: x_shift_
  type(GridFunction)                             :: srcMask_, dstMask_
  type(GridFunction)                             :: shifts_
  type(GridFunction), dimension(3,3)             :: gm1_,km1_,gp1_,kp1_
  integer, dimension(:)                          :: m
  integer                                        :: nl, num_int_points
  integer, dimension(:)                          :: shift_p1,shift_p2,shift_p3 
  integer, dimension(:)                          :: shift_b


  real*8,  dimension(3)                          :: h3d
  integer, dimension(3)                          :: ext
  integer                                        :: i,j


  
  interface
     subroutine preshift(h,dx,dy,dz,srcMask,dstMask,shifts,ex,shift_p1, &
          shift_p2,shift_p3,shift_b,params,num_int_points, &
          mask_list,num_mask_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 
       integer, dimension(:)    :: shift_b
       integer :: num_int_points  ! the total # of points to shift
       integer, dimension(:)                      :: mask_list
       integer                  ::  num_mask_points
       real*8 :: EPSRND,dxeps(3),small
       real*8,dimension(:) :: params
       parameter (EPSRND = 1.0d06,small=1.0e-10)
     end subroutine preshift

     subroutine preshift1o(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,parameter :: Nshifts=728
       integer, dimension(Nshifts),save :: shiftx,shifty,shiftz
       integer, dimension(:) :: m
       integer               :: nl,num_int_points
       integer,dimension(:)  :: shift_p1,shift_p2,shift_p3,shift_b
       integer,dimension(:)  :: ex
       real*8,dimension(:)   :: params
       real*8,dimension(:,:,:) :: dx,dy,dz
       real*8,dimension(:,:,:) :: srcmask,dstmask,shifts
       real*8,dimension(3)   :: h,hm1
     end subroutine

     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
       real*8,  dimension(3)                 :: h
       integer, dimension(3)                 :: ex
       real*8,  dimension(:,:,:) :: U,dx,dy,dz,Mask,shifts,dest
       integer num_int_points
       integer,dimension(:) :: shift_p1,shift_p2,shift_p3,shift_b
                 real*8,dimension(:) :: params
     end subroutine make_xgrid

  end interface


  h3d=params(P_GRIDSPACING)
  ext=shape(srcMask_%d)

  call preshift(h3d,x_shift_(1)%d,x_shift_(2)%d,x_shift_(3)%d,      &
       & srcMask_%d,dstMask_%d,shifts_%d,ext,shift_p1,shift_p2,shift_p3, &
       & shift_b,params,num_int_points,m,nl)
 
     !---------------------------\
     ! Now interpolate variables  \
     ! to their final locations    >
     ! Use kavg, gavg as temps    /
     !---------------------------/
  do i=1,3
    do j=i,3
      call make_xgrid(h3d,gp1_(i,j)%d,                         &
     &          x_shift_(1)%d,x_shift_(2)%d,x_shift_(3)%d, &
     &          srcMask_%d,shifts_%d,gm1_(i,j)%d,ext,     &
     &          shift_p1,shift_p2,shift_p3,                &
     &          shift_b,num_int_points,params)
      gp1_(i,j)%d = gm1_(i,j)%d
      call make_xgrid(h3d,kp1_(i,j)%d,                         &
     &          x_shift_(1)%d,x_shift_(2)%d,x_shift_(3)%d, &
     &          srcMask_%d,shifts_%d,km1_(i,j)%d,ext,        &
     &          shift_p1,shift_p2,shift_p3,                &
     &          shift_b,num_int_points,params)
      kp1_(i,j)%d = km1_(i,j)%d
    end do
  end do


return
end
