real function trilinear_interp_var(var, cell, xgrid, ygrid, zgrid, x, y, z)
use modglobal, only : ib, ie, ih, jb, je, jh, kb, ke, kh, itot, jtot, ktot
use decomp_2d, only : zstart
implicit none
real, intent(in) :: var(ib-ih:ie+ih,jb-jh:je+jh,kb-kh:kb+kh)
integer, intent(in) :: cell(3) ! GLOBAL indices of cell containing the point
real, intent(in), dimension(ib:itot+ih) :: xgrid
real, intent(in), dimension(jb:jtot+jh) :: ygrid
real, intent(in), dimension(kb:ktot+kh) :: zgrid
real, intent(in) :: x, y, z ! location of point to interpolate at
real, dimension(8) :: corners(8)
real :: x0, y0, z0, x1, y1, z1
integer :: i, j, k
i = cell(1) - zstart(1) + 1
j = cell(2) - zstart(2) + 1
k = cell(3) - zstart(3) + 1
if ((i < ib-1) .or. (i > ie+1) .or. (j < jb-1) .or. (j > je+1)) then
write(*,*) "problem in trilinear_interp_var", i, j, k
stop 1
end if
corners = eval_corners(var, i, j, k)
x0 = xgrid(cell(1))
y0 = ygrid(cell(2))
z0 = zgrid(cell(3))
x1 = xgrid(cell(1)+1)
y1 = ygrid(cell(2)+1)
z1 = zgrid(cell(3)+1)
trilinear_interp_var = trilinear_interp(x, y, z, x0, y0, z0, x1, y1, z1, corners)
end function trilinear_interp_var