pro co2_equilib,p,t,qco2,t_co2,p_co2,v_co2

;-------------------------------------------------------------
;
; Purpose:
;
;  Describe the equilibrium over solid CO2,  the vapor pressures
;  and solid co2 density are from the CRC handbook
;
; Input:
;     p.......pressure in mb
;     t.......temp in K
;     qco2....co2 mixing ratio in ppmv
; Output:
;     t_co2...equilibrium formation temp of solid co2 (K)
;     p_co2...equilibrium vapor pressure over solid co2 (mb)
;     v_co2...volume of condensed solid co2 (um3/cm3)
;
; Source: Mark Hervig
;
;-------------------------------------------------------------

;- some constants

   R      = 8.314         ; J/mol/K
   Mwc    = 44.0          ; molec wt. CO2, g/mol
   dc     = 1.56          ; density of solid CO2, g/cm3
   t_melt = 273.16 - 57.5 ; co2 melting point, K

;- check the input conditions

   t_co2 = -999.
   p_co2 = -999.
   v_co2 = -999.

   if (t gt t_melt) then print, 'T > CO2 melting point in co2_equilib.pro'
   if (t gt t_melt) then return
 
;- psat vs. tsat describes the saturation vapor pressure 
;  of co2 gas over solid co2 vs. temperature.

   tsat = [-134.3,-119.5,-108.6,-100.2,-85.7,-78.2]  ; deg C
   psat = [1.0,10.0,40.0,100.0,400.0,760.0]          ; mm Hg

   tsat = tsat +273.16     ; convert to K
   psat = psat * 0.750062  ; convert to mbar

;- find t_co2 and p_co2, using the fact that
;  a linear relationship exists between 1/T and alog10(p_co2).
;  Do a linear interpolation between some hi and low temps

   np = n_elements(psat)-1

   p1 = alog10(psat(0))
   p2 = alog10(psat(np))
   pp = p * qco2 * 1.e-6   ; partial pressure of co2

   ta = 1./tsat(0)
   tb = 1./tsat(np)
   tc = 1./t

   dtdt  = (tb - tc ) / (tb - ta)
   p_co2 = 10.0 ^ (p2 - dtdt * (p2 - p1))

   dsds  = (p2 - alog10(pp) ) / (p2 - p1)
   t_co2 = 1. / ( tb - dsds * (tb - ta) )

;- if t is below t_co2,  find the volume of condensate

   if (t le t_co2) then begin

     q_xs  = qco2*1.e-6 - p_co2/p      ; excess co2 mixing ratio
     nn    = (p*100.0) * q_xs/(R*T)    ; mols condensed co2 per m3 air
     v_co2 = 1.0e6 * nn * Mwc / dc     ; solid co2 volume, um3/cm3

   endif

   return
   end