anchor.diff:          ri=reshape(r_arr(indx),(/npts,2**(ndim+1)/))
           thi=reshape(th_arr(indx),(/npts,2**(ndim+1)/))
           phii=reshape(ph_arr(indx),(/npts,2**(ndim+1)/))
  -!        write(6,*) 'after reshape', minval(ttd), maxval(ttd)
  +                                                                                               !        write(6,*) 'after reshape', minval(ttd), maxval(ttd)
   !        rttd=ttd
           rttd=0.
  -        rho=merge(interp(rhoi,rttd,pd,rd,td),dzero,x1.gt.uniqx1(1))
  -        p=merge(interp(ppi,rttd,pd,rd,td),fone,x1.gt.uniqx1(1))
  +        rho=merge(interp(rhoi,rttd,pd,rd,td),dzero,x1.gt.uniqx1(1))*nfac
  +        p=merge(interp(ppi,rttd,pd,rd,td),fone,x1.gt.uniqx1(1))*pfac
   !        write(6,*) 'rho: ', rho, p
           vrl0=merge(interp(vrli,rttd,pd,rd,td),dzero,x1.gt.uniqx1(1))
           vtl0=merge(interp(vtli,rttd,pd,rd,td),dzero,x1.gt.uniqx1(1))
  
positive.code:         subroutine harm3d_vals(x0,a,rho,p,b,u,bmag)
          type (four_Vector), intent(in), dimension(:) :: x0
          real, intent(in) :: a
          real(kind=8), dimension(size(x0)) :: done
          real, dimension(size(x0)) :: x3,x2,x1,zm,theta,fone, &
           vpl0,vrl0,vtl0,rd,td, pd,rttd,zr,dzero, &
           vr0,vth0,vph0,bth,dummy,tt,ttd,zt,zphi,zpp
  
negative.code: ! 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.)))
             rd=1.
             pfac=1e-3
             nfac=1e-3
          elsewhere
             pfac=1.
             nfac=1.
          endwhere
  !        write(6,*) 'coords: ',
  !        write(6,*) 'rd td pd: ',minval(rd),maxval(rd),minval(td),maxval(td),minval(pd),maxval(pd)
  