function nice_length_unit, length, lengthlabel=lengthlabel

  factor = 3.086d27
  lengthlabel = 'Gpc'
  if length lt 3.086d27 then begin
     factor = 3.086d24
     lengthlabel = 'Mpc'
  end
  if length lt 3.086d24 then begin
     factor = 3.086d21
     lengthlabel = 'kpc'
  end
  if length lt 3.086d21 then begin
     factor = 3.086d18
     lengthlabel = 'pc'
  end
  if length lt 3.086d16 then begin
     factor = 1.49598d13
     lengthlabel = 'AU'
  end
  if length lt 1.49598d12 then begin
     factor = 6.69d10
     lengthlabel = 'Rsun'
  end
  if length lt 6.69d10 then begin
     factor = 1d5
     lengthlabel = 'km'
  end
  if length lt 1d5 then begin
     factor = 1d2
     lengthlabel = 'm'
  end
  if length lt 1d2 then begin
     factor = 1.D
     lengthlabel = 'cm'
  end

return, length/factor
end

function nice_time_unit, time, timelabel=timelabel
  factor =   3.15576e+16
  timelabel = 'Gyr'
  if time lt 3.15576e+16 then begin
     factor = 3.15576e+13
     timelabel = 'Myr'
  end
  if time lt 3.15576e+13 then begin
     factor = 3.15576e+10
     timelabel = 'kyr'
  end
  if time lt 3.15576e+10 then begin
     factor = 3.15576e+7
     timelabel = 'yr'
  end
  if time lt 3.15576e7 then begin
     factor = 86400.
     timelabel = 'day'
  end
  if time lt 86400. then begin
     factor = 1.
     timelabel = 's'
  end
  if time lt 0.1 then begin
     factor = 1.e-3
     timelabel = 'ms'
  end
return, time/factor
end


function get_unit, label, ulabel=ulabel, force_use_units=force_use_units, value=valuein
@common_blocks.inc
;value: use this code value to pick nice time and length units

  zz = 0.
  if N_elements(valuein) ne 1 then value = 1. else value = abs(valuein); default
  if (N_elements(redshifts) eq N_elements(times) ) then $
     if redshifts[time_index] ne -1 then zz=redshifts[time_index]
  protonmass = 1.67e-24
  tlabel = label
  if (strmatch(label, '*adius'))[0] ge 1 then tlabel = 'Length'
  if (where(label eq ['x','y','z','dx','dy','dz']))[0] ge 0 then tlabel = 'Length'
  if (strmatch(label, '*ensity'))[0] ge 1 then tlabel = 'Density'
  if (strmatch(label, '*mass*'))[0] ge 1 then tlabel = 'mass'
  if (strmatch(label, '*_time'))[0] ge 1 then tlabel = 'Time'
  if (strmatch(label, '* age'))[0] ge 1 then tlabel = 'Time'
  if (strmatch(label, '*_position*'))[0] ge 1 then tlabel = 'Length'
  if (((strmatch(label, '*elocity*'))[0] ge 1) or strmatch(label,'*speed')) and not (strmatch(label, 'Velocity_Div'))[0] then tlabel = 'Velocity'
  if (strmatch(label, '*Energy'))[0] ge 1 then tlabel = 'Energy'
  if (strmatch(label, 'Tb*'))[0] ge 1 then tlabel = 'Tb'
  if (strmatch(label, '*Balmer*'))[0] ge 1 then tlabel = 'Balmer'
  if (strmatch(label, '*Intensity*'))[0] ge 1 then tlabel = 'Intensity'
  if (strmatch(label, '*_kph'))[0] ge 1 then tlabel = 'RadiativeRate'
  if label eq 'Bx' or label eq 'By' or label eq 'Bz' or label eq 'abs(B)' then $
     tlabel = 'Bfield'
;  print, tlabel
  case 1 of 
     (tlabel eq 'Time'): Begin
;        unitfactor = timeunit/3.156e13
;        unitlabel = 'Myr'
        unitfactor = nice_time_unit(value*timeunit, timelabel=unitlabel)/value
     end
     (tlabel eq "Cooling_Time"): BEGIN
        print, value
        unitfactor = nice_time_unit(value, timelabel=unitlabel)/value
        print, unitfactor, " :: ", unitlabel
     END
     (tlabel eq 'Density'): BEGIN
        unitfactor = densityunit/protonmass*(1.+zz)^3
        unitlabel  = textoidl('amu. cm^{-3}')
        if (unitfactor gt 1.d20) then begin  ; instead of showing huge number density values use g/cm^3
           unitfactor = densityunit*(1.+zz)^3
           unitlabel  = 'g cm^{-3}'
        endif
     END
     (tlabel eq 'Length'): Begin
;        unitfactor = lengthunit/(1.+zz)
;        unitlabel = 'cm'
        unitfactor = nice_length_unit(value*lengthunit/(1.+ zz), lengthlabel=unitlabel)/value
     end
     (tlabel eq 'volume'): Begin
        unitfactor = nice_length_unit(value^(1./3.)*lengthunit/(1.+ zz), lengthlabel=unitlabel)
        unitfactor = unitfactor^3
        unitlabel = textoidl(unitlabel+'^3')
     end

     (tlabel eq 'Velocity'): Begin
        unitfactor = lengthunit/timeunit/1e5/(1.+InitialRedshift)
        unitlabel = 'km/s'
     end
     (tlabel eq 'Energy'): Begin
        unitfactor = (lengthunit/timeunit)^2/((1.+InitialRedshift))^2 * protonmass *6.2415e11
        unitlabel = 'eV/amu.'
        if (densityunit/protonmass*(1.+zz)^3 gt 1.d20) then begin  ; instead of showing huge energy values use erg/g
           unitfactor = (lengthunit/timeunit)^2/((1.+InitialRedshift))^2 
           unitlabel = 'erg/g'
        endif

     end
     (tlabel eq 'Temperature'): Begin
        unitfactor = temperatureunit
        unitfactor = 1.
        unitlabel = 'K'

     end
     (tlabel eq 'Pressure'): Begin
        eunit = (lengthunit/timeunit)^2/(1.+InitialRedshift)^2/(1+zz)^2
        kb = 1.38D-16
        m_p = 1.67e-24
        unitfactor = densityunit *(1.+zz)^3 * eunit/kb
        unitlabel = textoidl('a.m.u. cm^{-3}  K')
     end
     (tlabel eq 'Bfield'): Begin
;        unitfactor = sqrt(densityunit*(lengthunit/timeunit)^2*4.*!pi)/1.e-6
;        unitlabel = 'micro Gauss'
        unitfactor = sqrt(densityunit*(1.+zz)^3* $
                          (lengthunit/(1.+zz)/timeunit)^2*4.*!pi)* $
                     ((1.+InitialRedshift)/(1.+zz))^(-1)
        unitlabel = 'Gauss'
     end
     (tlabel eq 'mass'): Begin
        unitfactor = densityunit*lengthunit^3/1.989d33
        unitlabel = textoidl('M')+'!D!9n!X!N'
     end
     (tlabel eq 'Dynamical time (gas)'): Begin
        unitfactor = 1./sqrt(densityunit*6.67259d-8*(1.+zz)^3)/3.156d7
        unitlabel = 'yrs'
     end
     (tlabel eq '|v_particles|^10'): begin
        unitfactor = (lengthunit/timeunit/1e5/(1.+InitialRedshift))^10
        unitlabel = textoidl('(km s^{-1})^{10}')
     end
     (tlabel eq 'TT12'): begin
        unitfactor = (lengthunit/timeunit/1e5/(1.+InitialRedshift))^4 * (densityunit*lengthunit^3/1.989d33)^2
        unitlabel = textoidl('(km s^{-1})^{4} ')+textoidl('M')+'!D!9n!X!N'+textoidl('^2')
     end

     (tlabel eq 'Tb'): begin
        unitfactor = 1.0
        unitlabel = 'mK'
     end
     (tlabel eq 'Balmer'): begin
        unitfactor = 1e-21
        unitlabel = textoidl('erg cm^{-3} s^{-1}')
     end
     (tlabel eq 'Intensity'): begin
        ; units are erg/cm2/s/sr, e.g.
        ; g / s^3 / sr
        ; = rho * (cm/s)^3 / (4pi)
        unitfactor = densityunit * (lengthunit/timeunit)^3 / (4.d*!pi)
        unitlabel  = textoidl('erg cm^{-2} s^{-1} sr^{-1}')
     end
     (tlabel eq 'RadiativeRate'): begin     
        unitfactor = 1./timeunit
        unitlabel  = textoidl('s^{-1}')
     end
     (tlabel eq 'radial dM/dt'): begin     
        unitfactor = densityunit * lengthunit^3/1.989d33/timeunit*3.156e7
        unitlabel  = textoidl('M')+'!D!9n!X!N'+textoidl('yr^{-1}')
     end
     (tlabel eq 'div V' or tlabel eq 'Velocity_Div'): begin     
        unitfactor = 1./timeunit*3.156e7
        unitlabel  = textoidl('yr^{-1}')
     end

     (tlabel eq 'Current'): begin     
        magunitfactor = sqrt(densityunit*(1.+zz)^3* $
                          (lengthunit/(1.+zz)/timeunit)^2*4.*!pi)* $
                     ((1.+InitialRedshift)/(1.+zz))^(-1)
        unitfactor = magunitfactor/lengthunit
        unitlabel  = 'Gauss/cm'
     end


     else: begin
        unitfactor = 1.
        unitlabel = ''
     end
  end

if not useunits and not keyword_set(force_use_units) then begin 
   unitlabel  = ''
   unitfactor = 1.
endif

if N_elements(ulabel) gt 0 then $ ; if you did  /ulabel  in the call 
   if ulabel eq 1 then return, unitlabel  ; then just return the label

ulabel = unitlabel
return, unitfactor

end