anchor.diff:  !       write(6,*) 'ux lx: ',minval(lx1),maxval(ux1),minval(lx2),maxval(ux2), minval(lx3), maxval(ux3)
   ! Deal with poles
   !        write(6,*) 'poles'
  -        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
  +        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=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)
  
positive.code: !        write(6,*) 'uniform phi'
  ! uniform in phi
          pd=(zphi-minph)/(uniqph(2)-uniqph(1))
          td=abs(theta-calcthmksbl3(dble(uniqx2(lx2)),dble(zr)))/dth
          rd=(zr-uniqr(lx1))/(uniqr(ux1)-uniqr(lx1))
  ! for now just do nearest neighbor
          rd(:)=1.; td(:)=1.; pd(:)=1.
  
negative.code:           uks%data(1)=fum%data(1)
          end function umks2uksbl3
  
          function umksh2uks(fum,r,x2) result(uks)
            ! Converts a Kerr-Schild spherical 4-velocity to a Boyer-Lindquist one.
            ! JAD 11/10/2008
            ! Do simple transformation first:
  