subroutine interpolate_var(igrid,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,&
ixImax3,ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,gf,x,xloc,gfloc)
use mod_global_parameters
integer, intent(in) :: igrid,ixImin1,ixImin2,ixImin3,&
ixImax1,ixImax2,ixImax3, ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,&
ixOmax3
double precision, intent(in) :: gf(ixImin1:ixImax1,&
ixImin2:ixImax2,ixImin3:ixImax3)
double precision, intent(in) :: x(ixImin1:ixImax1,ixImin2:ixImax2,&
ixImin3:ixImax3,1:ndim)
double precision, intent(in) :: xloc(1:3)
double precision, intent(out) :: gfloc
double precision :: xd1,xd2,xd3, myq, dx1
double precision :: c0, c1, c00, c10, c01, c11
integer :: ic1,ic2,ic3, ic11,ic12,ic13, ic21,&
ic22,ic23, idir
if(stretch_type(1)==stretch_uni) then
! uniform stretch from xprobmin
ic1 = ceiling(dlog((xloc(1)-xprobmin1)*(qstretch(ps(igrid)%level,&
1)-one)/dxfirst(ps(igrid)%level,1)+one)/dlog(qstretch(ps(igrid)%level,&
1)))-(node(pig1_,igrid)-1)*block_nx1+nghostcells
else if(stretch_type(1)==stretch_symm) then
! symmetric stretch about 0.5*(xprobmin+xprobmax)
if(xloc(1)<xprobmin1+xstretch1) then
! stretch to left from xprobmin+xstretch
ic1 = block_nx1-(int(dlog((xprobmin1+&
xstretch1-xloc(1))*(qstretch(ps(igrid)%level,&
1)-one)/dxfirst(ps(igrid)%level,&
1)+one)/dlog(qstretch(ps(igrid)%level,&
1)))-(nstretchedblocks(ps(igrid)%level,1)/2-node(pig1_,&
igrid))*block_nx1)+nghostcells
else if(xloc(1)>xprobmax1-xstretch1) then
! stretch to right from xprobmax-xstretch
ic1 = ceiling(dlog((xloc(1)-xprobmax1+&
xstretch1)*(qstretch(ps(igrid)%level,&
1)-one)/dxfirst(ps(igrid)%level,&
1)+one)/dlog(qstretch(ps(igrid)%level,1)))-((node(pig1_,&
igrid)-1+nstretchedblocks(ps(igrid)%level,&
1)/2)*block_nx1-domain_nx1*2**(ps(igrid)%level-1))+nghostcells
else
! possible non-stretched central part
ic1 = int((xloc(1)-rnode(rpxmin1_,igrid))/rnode(rpdx1_,&
igrid)) + 1 + nghostcells
end if
else
! no stretch
ic1 = int((xloc(1)-rnode(rpxmin1_,igrid))/rnode(rpdx1_,&
igrid)) + 1 + nghostcells
end if
if(stretch_type(2)==stretch_uni) then
! uniform stretch from xprobmin
ic2 = ceiling(dlog((xloc(2)-xprobmin2)*(qstretch(ps(igrid)%level,&
2)-one)/dxfirst(ps(igrid)%level,2)+one)/dlog(qstretch(ps(igrid)%level,&
2)))-(node(pig2_,igrid)-1)*block_nx2+nghostcells
else if(stretch_type(2)==stretch_symm) then
! symmetric stretch about 0.5*(xprobmin+xprobmax)
if(xloc(2)<xprobmin2+xstretch2) then
! stretch to left from xprobmin+xstretch
ic2 = block_nx2-(int(dlog((xprobmin2+&
xstretch2-xloc(2))*(qstretch(ps(igrid)%level,&
2)-one)/dxfirst(ps(igrid)%level,&
2)+one)/dlog(qstretch(ps(igrid)%level,&
2)))-(nstretchedblocks(ps(igrid)%level,2)/2-node(pig2_,&
igrid))*block_nx2)+nghostcells
else if(xloc(2)>xprobmax2-xstretch2) then
! stretch to right from xprobmax-xstretch
ic2 = ceiling(dlog((xloc(2)-xprobmax2+&
xstretch2)*(qstretch(ps(igrid)%level,&
2)-one)/dxfirst(ps(igrid)%level,&
2)+one)/dlog(qstretch(ps(igrid)%level,2)))-((node(pig2_,&
igrid)-1+nstretchedblocks(ps(igrid)%level,&
2)/2)*block_nx2-domain_nx2*2**(ps(igrid)%level-1))+nghostcells
else
! possible non-stretched central part
ic2 = int((xloc(2)-rnode(rpxmin2_,igrid))/rnode(rpdx2_,&
igrid)) + 1 + nghostcells
end if
else
! no stretch
ic2 = int((xloc(2)-rnode(rpxmin2_,igrid))/rnode(rpdx2_,&
igrid)) + 1 + nghostcells
end if
if(stretch_type(3)==stretch_uni) then
! uniform stretch from xprobmin
ic3 = ceiling(dlog((xloc(3)-xprobmin3)*(qstretch(ps(igrid)%level,&
3)-one)/dxfirst(ps(igrid)%level,3)+one)/dlog(qstretch(ps(igrid)%level,&
3)))-(node(pig3_,igrid)-1)*block_nx3+nghostcells
else if(stretch_type(3)==stretch_symm) then
! symmetric stretch about 0.5*(xprobmin+xprobmax)
if(xloc(3)<xprobmin3+xstretch3) then
! stretch to left from xprobmin+xstretch
ic3 = block_nx3-(int(dlog((xprobmin3+&
xstretch3-xloc(3))*(qstretch(ps(igrid)%level,&
3)-one)/dxfirst(ps(igrid)%level,&
3)+one)/dlog(qstretch(ps(igrid)%level,&
3)))-(nstretchedblocks(ps(igrid)%level,3)/2-node(pig3_,&
igrid))*block_nx3)+nghostcells
else if(xloc(3)>xprobmax3-xstretch3) then
! stretch to right from xprobmax-xstretch
ic3 = ceiling(dlog((xloc(3)-xprobmax3+&
xstretch3)*(qstretch(ps(igrid)%level,&
3)-one)/dxfirst(ps(igrid)%level,&
3)+one)/dlog(qstretch(ps(igrid)%level,3)))-((node(pig3_,&
igrid)-1+nstretchedblocks(ps(igrid)%level,&
3)/2)*block_nx3-domain_nx3*2**(ps(igrid)%level-1))+nghostcells
else
! possible non-stretched central part
ic3 = int((xloc(3)-rnode(rpxmin3_,igrid))/rnode(rpdx3_,&
igrid)) + 1 + nghostcells
end if
else
! no stretch
ic3 = int((xloc(3)-rnode(rpxmin3_,igrid))/rnode(rpdx3_,&
igrid)) + 1 + nghostcells
end if
! linear interpolation:
if (x(ic1,ic2,ic3,1) .lt. xloc(1)) then
ic11 = ic1
else
ic11 = ic1 -1
end if
ic21 = ic11 + 1
if (x(ic1,ic2,ic3,2) .lt. xloc(2)) then
ic12 = ic2
else
ic12 = ic2 -1
end if
ic22 = ic12 + 1
if (x(ic1,ic2,ic3,3) .lt. xloc(3)) then
ic13 = ic3
else
ic13 = ic3 -1
end if
ic23 = ic13 + 1
if(ic11.lt.ixGlo1+1 .or. ic21.gt.ixGhi1-1) then
print *, 'direction: ',1
print *, 'position: ',xloc(1:ndim)
print *, 'indices:', ic11,ic21
call mpistop('Trying to interpolate from out of grid!')
end if
if(ic12.lt.ixGlo2+1 .or. ic22.gt.ixGhi2-1) then
print *, 'direction: ',2
print *, 'position: ',xloc(1:ndim)
print *, 'indices:', ic12,ic22
call mpistop('Trying to interpolate from out of grid!')
end if
if(ic13.lt.ixGlo3+1 .or. ic23.gt.ixGhi3-1) then
print *, 'direction: ',3
print *, 'position: ',xloc(1:ndim)
print *, 'indices:', ic13,ic23
call mpistop('Trying to interpolate from out of grid!')
end if
xd1 = (xloc(1)-x(ic11,ic12,ic13,1)) / (x(ic21,ic12,ic13,1) - x(ic11,ic12,&
ic13,1))
xd2 = (xloc(2)-x(ic11,ic12,ic13,2)) / (x(ic11,ic22,ic13,2) - x(ic11,ic12,&
ic13,2))
xd3 = (xloc(3)-x(ic11,ic12,ic13,3)) / (x(ic11,ic12,ic23,3) - x(ic11,ic12,&
ic13,3))
c00 = gf(ic11,ic12,ic13) * (1.0d0 - xd1) + gf(ic21,ic12,ic13) * xd1
c10 = gf(ic11,ic22,ic13) * (1.0d0 - xd1) + gf(ic21,ic22,ic13) * xd1
c01 = gf(ic11,ic12,ic23) * (1.0d0 - xd1) + gf(ic21,ic12,ic23) * xd1
c11 = gf(ic11,ic22,ic23) * (1.0d0 - xd1) + gf(ic21,ic22,ic23) * xd1
c0 = c00 * (1.0d0 - xd2) + c10 * xd2
c1 = c01 * (1.0d0 - xd2) + c11 * xd2
gfloc = c0 * (1.0d0 - xd3) + c1 * xd3
end subroutine interpolate_var