subroutine interpolate_oneblock(x,iw,out)
double precision, dimension(3),intent(in) :: x
integer, intent(in) :: iw
double precision, intent(out) :: out
! .. local ..
double precision, dimension(3) :: xloc
integer :: ic1,ic2,ic3, ic11,&
ic12,ic13, ic21,ic22,ic23
double precision :: xd1,xd2,xd3
double precision :: c0, c1, c00, c10,&
c01, c11
integer :: ipivot1,ipivot2,&
ipivot3, idir
xloc=x
!--------------------------------------------
! Hunt for the index closest to the point
! This is a bit slow but allows for stretched grids
! (still need to be orthogonal for interpolation though)
!--------------------------------------------
ipivot1=1;ipivot2=1;ipivot3=1;
do idir = 1, 3
select case (idir)
case (1)
ic1 = minloc(dabs(xloc(1)-xoneblock(:,ipivot2,ipivot3,idir)),1,&
mask=.true.)
case (2)
ic2 = minloc(dabs(xloc(2)-xoneblock(ipivot1,:,ipivot3,idir)),1,&
mask=.true.)
case (3)
ic3 = minloc(dabs(xloc(3)-xoneblock(ipivot1,ipivot2,:,idir)),1,&
mask=.true.)
case default
call mpistop("error1 in interpolate_oneblock")
end select
end do
! flat interpolation would simply be:
!out = woneblock(ic^D,iw)
!return
!-------------------------------------------
! Get the left and right indices
!-------------------------------------------
if (xoneblock(ic1,ic2,ic3,1) .lt. xloc(1)) then
ic11 = ic1
else
ic11 = ic1 -1
end if
ic21 = ic11 + 1
if (xoneblock(ic1,ic2,ic3,2) .lt. xloc(2)) then
ic12 = ic2
else
ic12 = ic2 -1
end if
ic22 = ic12 + 1
if (xoneblock(ic1,ic2,ic3,3) .lt. xloc(3)) then
ic13 = ic3
else
ic13 = ic3 -1
end if
ic23 = ic13 + 1
!--------------------------------------------
! apply flat interpolation if outside of range,
! change point-location to make this easy!
!--------------------------------------------
if (ic11 .lt. 1) then
ic11 = 1
ic21 = ic11 + 1
xloc(1) = xoneblock(ic1,ic2,ic3,1)
end if
if (ic12 .lt. 1) then
ic12 = 1
ic22 = ic12 + 1
xloc(2) = xoneblock(ic1,ic2,ic3,2)
end if
if (ic13 .lt. 1) then
ic13 = 1
ic23 = ic13 + 1
xloc(3) = xoneblock(ic1,ic2,ic3,3)
end if
if (ic21 .gt. nc1) then
ic21 = nc1
ic11 = ic21 - 1
xloc(1) = xoneblock(ic1,ic2,ic3,1)
end if
if (ic22 .gt. nc2) then
ic22 = nc2
ic12 = ic22 - 1
xloc(2) = xoneblock(ic1,ic2,ic3,2)
end if
if (ic23 .gt. nc3) then
ic23 = nc3
ic13 = ic23 - 1
xloc(3) = xoneblock(ic1,ic2,ic3,3)
end if
!-------------------------------------------
! linear, bi- and tri- linear interpolations
!-------------------------------------------
xd1 = (xloc(1)-xoneblock(ic11,ic12,ic13,1)) / (xoneblock(ic21,ic12,ic13,&
1) - xoneblock(ic11,ic12,ic13,1))
xd2 = (xloc(2)-xoneblock(ic11,ic12,ic13,2)) / (xoneblock(ic11,ic22,ic13,&
2) - xoneblock(ic11,ic12,ic13,2))
xd3 = (xloc(3)-xoneblock(ic11,ic12,ic13,3)) / (xoneblock(ic11,ic12,ic23,&
3) - xoneblock(ic11,ic12,ic13,3))
c00 = woneblock(ic11,ic12,ic13,iw) * (1.0d0 - xd1) + woneblock(ic21,ic12,&
ic13,iw) * xd1
c10 = woneblock(ic11,ic22,ic13,iw) * (1.0d0 - xd1) + woneblock(ic21,ic22,&
ic13,iw) * xd1
c01 = woneblock(ic11,ic12,ic23,iw) * (1.0d0 - xd1) + woneblock(ic21,ic12,&
ic23,iw) * xd1
c11 = woneblock(ic11,ic22,ic23,iw) * (1.0d0 - xd1) + woneblock(ic21,ic22,&
ic23,iw) * xd1
c0 = c00 * (1.0d0 - xd2) + c10 * xd2
c1 = c01 * (1.0d0 - xd2) + c11 * xd2
out = c0 * (1.0d0 - xd3) + c1 * xd3
end subroutine interpolate_oneblock