anchor.diff:          real(kind=8), dimension(size(f%rho)), intent(out) :: ncgs,ncgsnth,bcgs,tempcgs
           real(kind=8), dimension(size(f%rho)) :: rhocgs,pcgs
           real(kind=8) :: lcgs,tcgs,mdot
  -        ! Converts Cosmos++ code units to standard cgs units. Follows Schnittman et al. (2006).
  -        ! JAD 11/26/2012 adapted from IDL code
  -        ! Black hole mass sets time and length scales:
  -!        write(6,*) 'convert mbh: ',sp%mbh
  -        lcgs=GC*sp%mbh*msun/c**2; tcgs=lcgs/c
  -!        write(6,*) 'lcgs: ',lcgs,tcgs
  -        ! Typical mb09 code mdot value.
           mdot=.003
  -        ! Now convert density using black hole to scale length/time, 
  -        ! accretion rate to scale torus mass:
  -!        write(6,*) 'convert mdot: ',mdot,sp%mdot
  -        rhocgs=sp%mdot/mdot/lcgs**3*tcgs*f%rho; ncgs=rhocgs/mp
  -!        write(6,*) 'n: ',sp%mdot/mdot/lcgs**3.*tcgs/mp
  -        ! Use this to convert pressure:
  -        pcgs=f%p*rhocgs/f%rho*c**2.
  -        ! Ideal gas temperature for single fluid (i.e., no separate e-/p):
  -        tempcgs=pcgs/ncgs/k
  -        ! And finally, bfield conversion is just square root of this:
  -        bcgs=f%bmag*sqrt(rhocgs/f%rho)*c
  -        ! Convert HL units to cgs:
  -        bcgs=bcgs*sqrt(4.*pi)
  -        ! non-thermal e- put in by hand
  +        call scale_sim_units(sp%mbh,sp%mdot,mdot,f%rho,f%p,f%bmag,ncgs, &
  +             bcgs,tempcgs)
           ncgsnth=ncgs
  -!        write(6,*) 'leaving convert', maxval(bcgs), maxval(ncgs), maxval(tempcgs)
  -!        write(6,*) 'convert b: ',bcgs
  -!        write(6,*) 'convert n: ',ncgs/1e7
  -!        write(6,*) 'convert temp: ',tempcgs/1e10
  -!        write(6,*) 'convert temp 2: ',f%p/f%rho*mp/k*c**2./1e10
  -!        write(6,*) 'convert mdot: ',sp%mdot, sp%mbh
  -!        write(6,*) 'convert bh: ',tcgs,lcgs
           end subroutine convert_fluidvars_harm
   
           subroutine convert_fluidvars_harm3d(f,ncgs,ncgsnth,bcgs,tempcgs,sp)
  
positive.code:            trat=Rhigh
             nrat=Nhigh
          endwhere
          
  !        Thetae_unit = (gam - 1d0)* mp / trat
  !        Thetae_unit = mp / trat
  
  ! changed to add 1+trat to have maximum T_e = T_tot / 2
          tempcgs=(f%p/f%rho)*mp*c*c/k/(1d0+trat)
  
          ! And finally, bfield conversion is just square root of this:
  
          !B_unit = CL * sqrt(4.*M_PI*RHO_unit);
          bcgs=f%bmag*sqrt(4.*pi*rhocgs/f%rho)*c
          ! Convert HL units to cgs:
  !        bcgs=bcgs*sqrt(4.*pi)
  ! scale ncgs to tune disk density manually
  !        ncgs=ncgs/nrat
          ! non-thermal e- put in by hand
  !        ncgsnth=ncgs
  !        if(any(isnan(tempcgs))) then
  !            write(6,*) 'temp: ',tempcgs
  !            write(6,*) 'trat: ',trat
  !            write(6,*) 'beta: ', beta
  !        endif
          where(f%bmag**2d0/f%rho.gt.1d0)
             ncgsnth=sp%jetalphaval*bcgs**2d0/8d0/pi/sp%gminval*(sp%p1-2d0)/(sp%p1-1d0)/8.2d-7
          elsewhere
             ncgsnth=0d0
          endwhere
          !        write(6,*) 'leaving convert', maxval(bcgs), maxval(ncgs), maxval(tempcgs)
          !        write(6,*) 'convert b: ',bcgs
          !        write(6,*) 'convert n: ',ncgs/1e7
          !        write(6,*) 'convert temp: ',tempcgs/1e10
          !        write(6,*) 'convert temp 2: ',f%p/f%rho*mp/k*c**2./1e10
          !        write(6,*) 'convert mdot: ',sp%mdot, sp%mbh
          !        write(6,*) 'convert bh: ',tcgs,lcgs
  !        write(6,*) 'harm3d convert params: ',sp%p1,sp%gminval,pi,sp%jetalphaval
  !        write(6,*) 'harm3d convert bcgs: ',minval(bcgs),maxval(bcgs)
  !        write(6,*) 'end convert fluidvars harm3d ncgsnth: ',minval(ncgsnth),maxval(ncgsnth)
          end subroutine convert_fluidvars_harm3d
  
          subroutine convert_fluidvars_ffjet(f,ncgs,ncgsnth,bcgs,tcgs,sp)
  
negative.code:         real(kind=8), dimension(size(f%rho)), intent(out) :: ncgs,ncgsnth,bcgs,tempcgs
          real(kind=8), dimension(size(f%rho)) :: rhocgs,pcgs
          real(kind=8) :: lcgs,tcgs,mdot
          mdot=.0013
          call scale_sim_units(sp%mbh,sp%mdot,mdot,f%rho,f%p,f%bmag,ncgs, &
               bcgs,tempcgs)
  ! non-thermal e- by hand
          ncgsnth=ncgs
          end subroutine convert_fluidvars_mb09
  
          subroutine convert_fluidvars_harm(f,ncgs,ncgsnth,bcgs,tempcgs,sp)
  