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:              u%data(2) = u%data(2)-gamma*beta%data(2)/alpha
               u%data(3) = u%data(3)-gamma*beta%data(3)/alpha
               u%data(4) = u%data(4)-gamma*beta%data(4)/alpha
               deallocate(alpha); deallocate(beta); deallocate(gamma)
               deallocate(zero)
               b%data(1) = u*b
               b%data(2) = (b%data(2)+b%data(1)*u%data(2))/u%data(1)
               b%data(3) = (b%data(3)+b%data(1)*u%data(3))/u%data(1)
               b%data(4) = (b%data(4)+b%data(1)*u%data(4))/u%data(1)
            endif
            deallocate(data)
  ! test code for these transformations
  !          allocate(bdotu(n)); allocate(udotu(n)); allocate(bdotb(n))
  !          udotu=abs(u*u+1.)
  !          write(6,*) 'udotu: ',maxval(udotu),udotu(30*nx1*nx2+50)
  !          bdotb = b*b; bdotu = b*u
  !          write(6,*) 'bdotu: ',maxval(abs(bdotu)),bdotu(30*nx1*nx2+50), &
  !               bdotb(30*nx1*nx2+50)
  !          deallocate(bdotu); deallocate(udotu); deallocate(bdotb)
            p=merge(p,pmin,p.gt.pmin)
            write(6,*) "min vals rho", minval(rho)
            write(6,*) "min vals p", minval(p)
  