pro nat_equilib,q_hno3,q_h2o,p_mb,t,tnat,p_hno3_eq,vnat
   
;------------------------------------------------------------------
; purpose:
; 
;   Routine calculates the:
;
;      a) formation temperature of Nitric Acid Trihydrate (NAT)
;      b) equilibrium vapor pressure of hno3 over NAT
;      c) volume of condensed NAT
;
;   as a function of the h2o mixing ratio, hno3 mixing ratio
;   and ambient pressure and temperature.  The equation is from 
;   Hanson and Mauersberger,  GRL, vol 15, no8, pp855-858, 1988.
;
;   This routine will accept input parameters as arrays,  as long
;   as they all have the same dimensions.  
;
; parameters:
;
;   input:
;     q_hno3.......hno3 vapor mixing ratio (ppbv)
;     q_h2o........water vapor mixing ratio (ppmv)
;     p_mb.........pressure (mb)
;     t............temperature (K)
;
;   output:
;     tnat........NAT formation temperature (K)
;     vnat........volume of condensed NAT (um3/cm3)
;     p_hno3_eq...equilibrium hno3 vapor presure (mb)
;
; source:    Mark Hervig 
;-------------------------------------------------------------------

;- check inputs

   i = n_elements(q_h2o)
   j = n_elements(q_hno3)
   k = n_elements(p_mb)

   if i ne j then stop,'input arrays mismatched in tnat.pro'
   if i ne k then stop,'input arrays mismatched in tnat.pro'
   if k ne j then stop,'input arrays mismatched in tnat.pro'

;- setup some constants etc...

   tnat      = fltarr(i)  ; NAT temperature
   vnat      = fltarr(i)  ; NAT volume

   R       = 8.314       ; J/mol/K
   Mwn     = 63.0        ; molec wt. hno3, g/mol
   dn      = 1.60        ; density of NAT, g/cm3
   wfn     = 0.5385      ; wt. fraction hno3 in NAT

   torr_mb = 1.333224    ; converts pressure,  torr to mb
   atm_mb  = 1013.0      ; converts pressure,  atmospheres to mb

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

   c1 = -2.7836
   c2 = 0.00088
   c3 = 38.9855
   c4 = 11397.0
   c5 = 0.009179

;- solve for equilibrium NAT formation temperature 

   term1 = c3 - alog10(p_hno3) + c1 * alog10(p_h2o)
   term2 = c5 - c2 * alog10(p_h2o) 

   a  = 1.0
   b  = term1 / term2
   c  = -1.* c4 / term2
   x1 = (-1.*b + sqrt(b^2 - 4.*a*c)) / (2.*a)
   x2 = (-1.*b - sqrt(b^2 - 4.*a*c)) / (2.*a)

   k  = where(x1 gt 0.,n)
   if (n gt 0) then tnat(k) = x1(k)

   k  = where(x2 gt 0.,n)
   if (n gt 0) then tnat(k) = x2(k)

;- solve for equilibrium pressure of hno3 over NAT,  in mbar

   p_hno3_eq = torr_mb*10.0^((c1-c2*T)*alog10(p_h2o)+c3-c4/T+c5*T)  ; mbar

   q_hno3_eq = p_hno3_eq / (p_mb *1.e-9)      ; ppbv

;- for T < tnat,  compute NAT volume

   k = where(t le tnat, n)

   if (n gt 0) then begin

     q_hno3_xs = q_hno3(k) - q_hno3_eq(k)                  ; excess hno3, ppbv
     nn        = p_mb(k)*100.0 *q_hno3_xs*1.e-9 /(R*t(k))  ; mols condensed hno3 per m3 air
     vnat(k)   = 1.0e6 * (nn * Mwn) / (wfn * dn)           ; NAT volume, um3/cm3

   endif

;- if inputs were scalar,  trim down output

   if i eq 1 then begin
     tnat = tnat(0)
     vnat = vnat(0)
     p_hno3_eq = p_hno3_eq(0)
   endif

   return
   end