anchor.diff:            endif
             r_arr=exp(x1_arr)
             ph_arr=x3_arr
  -! HARD-CODED MKS or MMKS
  -          th_arr=calcthmmks(x2_arr,x1_arr)
  -!          th_arr=real(calcthmksh(dble(x2_arr)))
  +          if(IS_MMKS.eq.1) then
  +             th_arr=calcthmmks(x2_arr,x1_arr)
  +          else
  +             th_arr=real(calcthmksh(dble(x2_arr)))
  +          endif
             allocate(p(nx1*nx2*nx3))
             allocate(rho(nx1*nx2*nx3)); allocate(pmin(nx1*nx2*nx3))
             allocate(u(nx1*nx2*nx3)); allocate(b(nx1*nx2*nx3))
  
positive.code:           r_arr=exp(x1_arr)
  ! the 2*pi factor is from my conversion choice from HDF5
            ph_arr=x3_arr*2.*pi
            th_arr=calcthmmks(x2_arr,x1_arr)
  ! CHANGING FOR REGULAR MKS
  !          th_arr=real(calcthmksh(dble(x2_arr)))
            write(6,*) 'iharm grid: ',maxval(x1_arr),maxval(x2_arr),minval(th_arr)
          end subroutine read_iharm_grid_file
  
  
negative.code:           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
            if(IS_MMKS.eq.1) then
               uks = ummks2uks(u,dble(x1_arr),dble(x2_arr))
               bks = ummks2uks(b,dble(x1_arr),dble(x2_arr))
            else
               uks = umksh2uks(u,r_arr,x2_arr)
               bks = umksh2uks(b,r_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)),&
  