anchor.diff:             pfac=1.
              nfac=1.
           endwhere
  +!  nearest neighbor
  +        rd(:)=1.; td(:)=1.; pd(:)=1.
   !        write(6,*) 'coords: ',
   !        write(6,*) 'rd td pd: ',minval(rd),maxval(rd),minval(td),maxval(td),minval(pd),maxval(pd)
   !        write(6,*) 'ux lx: ',minval(lx1),maxval(ux1),minval(lx2),maxval(ux2),minval(lx3),maxval(ux3)
  
positive.code: !        write(6,*) 'nx1 nx2 nx3: ',nx1,nx2,nx3,minval(zphi),maxval(zphi),minval(zr),maxval(zr),minval(theta), &
  !             maxval(theta)
  !       write(6,*) 'ux lx: ',minval(lx1),maxval(ux1),minval(lx2),maxval(ux2), minval(lx3), maxval(ux3)
  ! Deal with poles
  !        write(6,*) 'poles'
          if(BL.eq.1) then
             where(ux2.ne.lx2)
                dth=uniqth(ux2)-uniqth(lx2)
             elsewhere
                dth=uniqth(1)
             endwhere
          else
  !           where(ux2.ne.lx2)
  !              dth=calcth_cylindrified(dble(uniqx2(ux2)),dble(zr))-calcth_cylindrified(dble(uniqx2(lx2)),dble(zr))
  !           elsewhere
  !              dth=calcth_cylindrified(dble(uniqx2(1))+0d0*dble(zr),dble(zr))
  !           endwhere
             where(ux2.ne.lx2)
                dth=calcthmksbl3(dble(uniqx2(ux2)),dble(zr))-calcthmksbl3(dble(uniqx2(lx2)),dble(zr))
             elsewhere
                dth=calcthmksbl3(dble(uniqx2(1))+0d0*dble(zr),dble(zr))
             endwhere
          end if
  ! periodic in phi
          minph=uniqph(lx3)
          where(ux3.gt.nx3)
              ux3=1
          endwhere
  
negative.code:         endwhere
  !        write(6,*) 'uniform phi'
  ! uniform in phi
  !        pd=(zphi-minph)/(uniqph(2)-uniqph(1))
  !        if(BL.eq.3) then
  !!           td=abs(theta-calcth_cylindrified(dble(uniqx2(lx2)),dble(zr)))/dth
  !           td=abs(theta-calcthmksbl3(dble(uniqx2(lx2)),dble(zr)))/dth
  !        else
  !           td=abs(theta-uniqth(lx2))/dth
   !       end if
   !       rd=(zr-uniqr(lx1))/(uniqr(ux1)-uniqr(lx1))
  ! When the geodesic is inside the innermost zone center located outside the horizon, 
  ! use nearest neighbor rather than interpolate
          where(uniqr(lx1).le.(1.+sqrt(1.-a**2.)))
  