anchor.diff:          uniqx2=x2_arr(1:nx3*(nx2-1)+1:nx3)
           uniqx1=x1_arr(1:nx3*(nx2-1)*nx1+1:nx3*nx2)
           uniqr=exp(uniqx1)
  -!        uniqth=pi*uniqx2
  +        if(BL.eq.1) uniqth=pi*uniqx2+(1.-hslope)/2.*sin(2.*pi*uniqx2)
           uniqph=uniqx3
   !        write(6,*) 'uniqr: ',minval(uniqr), maxval(uniqr)
   !        write(6,*) 'uniqth: ',minval(uniqth), maxval(uniqth)
  
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:           write(6,*) 'read harm transform coords phi',minval(x3_arr), maxval(x3_arr)
            write(6,*) 'read harm transform coords u ',minval(u%data(1)),minval(u%data(2))
            write(6,*) 'read harm transform coords u size ',size(u),hslope
            if(BL.eq.3) then
               uks = umks2uksbl3(u,dble(x1_arr),dble(x2_arr))
               bks = umks2uksbl3(b,dble(x1_arr),dble(x2_arr))
            else
               uks = umksh2uks(u,x1_arr,x2_arr)
               bks = umksh2uks(b,x1_arr,x2_arr)
            end if
            write(6,*) 'after uks ',minval(uks%data(1))
            u = uks2ubl(uks,dble(r_arr),dble(asim))
            write(6,*) 'read harm transform coords b', minval(u%data(3)),maxval(u%data(2)),&
                 minval(b%data(1)),maxval(b%data(4))
            write(6,*) "min vals u", minval(u%data(1)), minval(u%data(2)), minval(u%data(3)), minval(u%data(4))
            write(6,*) "max vals u", maxval(u%data(1)), maxval(u%data(2)), maxval(u%data(3)), maxval(u%data(4))
            b   = uks2ubl(bks,dble(r_arr),dble(asim))
  ! test code...
            allocate(bdotu(n)); allocate(udotu(n)); allocate(bdotb(n))
  