anchor.diff:  ! 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.
  +           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)
  
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:         u%data(4)=u%data(1)*dble(vph0)
          call assign_metric(u,transpose(kerr_metric(zr,real(x0%data(3)) &
          ,a)))
          write(6,*) 'min vals u', minval(u%data(1)), maxval(abs(u*u+1))
  !        write(6,*) 'harm vals r: ',zr
  !        write(6,*) 'harm vals rho: ',rho
  !        write(6,*) 'harm vals rd: ',rd
  !        write(6,*) 'harm vals ux: ',nx1,ux1
  !        write(6,*) 'harm vals udotu: ',(abs(u*u+1))(minloc(zr))
          write(6,*) 'harm vals minloc: ',minloc(zr),rd(minloc(zr)), &
               td(minloc(zr)),x0(minloc(zr))%data(3),uniqth(lx2(minloc(zr))), &
               lx2(minloc(zr)),pd(minloc(zr)),lx1(minloc(zr)),ux1(minloc(zr)), &
               uniqx1(1),x1(minloc(zr))
          write(6,*) 'harm vals minloc b0i: ',b0i(minloc(zr),:)
          write(6,*) 'harm vals minloc bri: ',bri(minloc(zr),:)
          write(6,*) 'harm vals minloc bthi: ',bthi(minloc(zr),:)
          write(6,*) 'harm vals minloc bphi: ',bphi(minloc(zr),:)
          write(6,*) 'harm vals interp bi: ',b(minloc(zr))%data(1),b(minloc(zr))%data(2), &
               b(minloc(zr))%data(3),b(minloc(zr))%data(4)
          write(6,*) 'harm vals minloc x0: ', ri(minloc(zr),:),thi(minloc(zr),:),&
               phii(minloc(zr),:)
  !        write(6,*) 'harm coords: ',zr(minloc(p)),theta(minloc(p)),zphi(minloc(p))
  !        write(6,*) 'harm vals rhoi: ',rhoi(minloc(p),:)
          end subroutine harm3d_vals
  