pro nad_equilib, q_hno3,q_h2o,p,t,t_nad,p_hno3_eq,v_nad
   
;------------------------------------------------------------------
;
; purpose:
; 
;   Routine calculates the:
;
;     a) equilibrium temperature of nitric acid dihydrate (NAD)
;     b) equilibrium vapor pressure of hno3 over NAD
;     c) volume of condensed NAD
;
;   as a function of the h2o mixing ratio, hno3 mixing ratio
;   and ambient pressure and temperature.  
;   The equilibrium over NAD is from Worsnop et al., Science, 259, 71-74, 1993.
;
;   and the data from Fig 3 were fit according to 
;      Kn = exp(a + bT) 
;      where b = 0.578, a = -144.15
;
; parameters:
;
;   input:
;     q_hno3.......hno3 vapor mixing ratio (ppbv)
;     q_h2o........water vapor mixing ratio (ppmv)
;     p............pressure (mb)
;     t............temperature (K)
;
;   output:
;     t_nad........NAD equilibrium temperature (K)
;     p_hno3_eq....equilibrium hno3 vapor presure (mb)
;     v_nad........volume of condensed NAT (um3/cm3)
;
; source:    Mark Hervig 
;
;-------------------------------------------------------------------

;- check inputs

   if (q_hno3 le 0.0 or q_h2o le 0.0 or p le 0.0 or t le 0.0) then begin
     t_nad=0.0
     p_hno3_eq=0.0
     v_nad=0.0
     return
   endif

;- setup some constants etc...

   p_hno3    = 0.0 
   p_h2o     = 0.0
   t_frost   = 0.0
   v_nad     = 0.0
   t_nad     = 0.0
   p_hno3_eq = 0.0

   if (q_hno3 le 0.0) then return 

;- some constants

   torr_mb = 1.333224    ; pressure,  torr to mb
   R       = 8.314       ; J/mol/K
   Mwd     = 99.0        ; molec wt. NAD, g/mol
   Mwn     = 63.0        ; molec wt. hno3, g/mol
   dn      = 1.67        ; density of NAD, g/cm3,  this number should be checked !!
   wfn     = 0.6364      ; wt. fraction hno3 in NAD
   a       = -144.15     ; constants describing Kn from Worsnop, fig 3
   b       = 0.578

;- partial pressures

   p_torr  = p / torr_mb
   p_hno3  = p_torr * q_hno3 * 1.0e-9
   p_h2o   = p_torr * q_h2o * 1.e-6

;- solve for equilibrium NAD temperature 

   t_nad = ( alog(p_hno3 * p_h2o^2) -a ) / b  ; K

;- solve for equilibrium pressure of hno3 over NAD

   Kn        = exp( a + b *T)              ; torr^3 
   p_hno3_eq = torr_mb * Kn * p_h2o ^(-2)  ; equilib p, mb
   q_hno3_eq = p_hno3_eq / (p *1.e-9)      ; equilib mix, ppbv

;- for T < t_nad,  compute NAD volume

   if (t le t_nad) then begin

     q_hno3_xs = q_hno3 - q_hno3_eq                   ; excess hno3
     nn        = (p*100.0) * (q_hno3_xs*1.e-9) /(R*T) ; mols condensed hno3 per m3 air
     v_nad     = 1.0e6 * (nn * Mwn) / (wfn * dn)      ; NAD volume, um3/cm3

   endif

   return
   end