function ice_density,t_in
   
;-------------------------------------------------------------------------------------
; Routine returns the density of ice for the input temperature, for 93K < T < 273K.
; 
; Input:
;  t_in.......temperature (K)
;
; Output:
;  den_out....ice density (g/cm3)
;
; Source: Mark Hervig, GATS Inc.
;-------------------------------------------------------------------------------------

;- ice density vs. T from http://www.engineeringtoolbox.com/ice-thermal-properties-d_576.html

   t   = [    0,   -5,  -10,  -15,  -20,  -25,  -30,  -35,  -40,  -50,  -60,  -70,  -80,  -90, -100, -180] + 273.16 ; K
   den = [916.2,917.5,918.9,919.4,919.4,919.6,920.0,920.4,920.8,921.6,922.4,923.3,924.1,924.9,925.7,934.0] * 1e-3   ; g/cm3

;- interpolate on that data

   den_out = -999.
   
   if (t_in ge min(t) and t_in le max(t)) then begin
   
     den_out = interpol(den,t,t_in)
     
   endif
   
;- done
   
   return,den_out
   end