function index_air,tk,pmb,pvmb,lambda

;---------------------------------------------------------------------
; Calculate the refractive index of moist air
; adapted from the paper:
;
; Phillip E. Ciddor, "Refractive index of air: new equations for the 
; visible and near infrared," Appl. Optics 35, 1566-1573 (1996).
;
; see the website:  http://emtoolbox.nist.gov/Wavelength/Ciddor.asp
;
; Input:
;   pmb.....pressure mb
;   TK......temperature, K
;   pvmb....water vapor pressure, mb
;   lambda..wavelength, microns
;
; Output:
;   index...refractive index of air
;
; Source: Mark Hervig
;---------------------------------------------------------------------

   XCO2 = 450.0 ; co2 mixing ratio, ppmv, vol mix = molal mix

   p   = pmb * 100.   ; pressure, mb -> pascals
   pv  = pvmb * 100.  ; vapor pressure, mb -> Pa
   TC  = TK - 273.15  ; T in deg C

;- define constants

   w0 = 295.235d0 
   w1 = 2.6422d0 
   w2 = -0.03238d0 
   w3 = 0.004028d0
   k0 = 238.0185d0 
   k1 = 5792105.d0 
   k2 = 57.362d0 
   k3 = 167917.d0
   a0 = 1.58123 * 1d-6 
   a1 = -2.9331 * 1d-8
   a2 = 1.1043 * 1d-10
   b0 = 5.707 * 1d-6
   b1 = -2.051 * 1d-8
   c0 = 1.9898 * 1d-4 
   c1 = -2.376 * 1d-6
   d  = 1.83 * 1d-11 
   e  = -0.765 * 1d-8

   pR1   = 101325.d0  ; Pa
   TR1   = 288.15d0   ; K
   Za    = 0.9995922115d0
   Rhovs = 0.00985938d0
   R     = 8.314472d0 
   Mv    = 0.018015d0

;- express humidity as a mole fraction, calculate mole fraction xv,
;  for partial pressure pv (Pascals), air pressure p (Pascals), 
;  and air temperature t (Celsius)

   Alpha = 1.00062 
   Beta   = 3.14 * 1e-8 
   Gamma = 5.60 * 1e-7

   xv = (Alpha + Beta * p + Gamma * TC^2) * pv / p

;- for wavelength in micrometers,  find
 
   S = 1 / Lambda^2.

;- Calculate intermediate results that depend on S

   ras = 1e-8 * ( (k1 / (k0 - S)) + (k3 / (k2 - S)) )
   rvs = 1.022 * 1e-8 * (w0 + w1 * S + w2 * S^2 + w3 * S^3)
	
;- Given the CO2 concentration xCO2 in micro-mol/mol, calculate

   Ma   = 0.0289635 + 1.2011 * 1e-8 * ( XCO2 - 400.)
   raxs = ras * ( 1 + 5.34 * 1e-7 * (XCO2 - 450.) )

;- Find the compressibility, and density components

   Zm = 1-(p/tk) *(a0+a1*tc+a2*tc^2 +(b0+b1*tc)*xv +(c0+c1*tc)*xv^2) +(p/tk)^2 *(d+e*xv^2)
   Rhoaxs = pR1 * Ma / (Za * R * TR1)
   Rhov = xv * p * Mv / (Zm * R * TK)
   Rhoa = (1 - xv) * p * Ma / (Zm * R * TK)

;- Calculate the index of refraction

   index = 1 + (rhoa / rhoaxs) * raxs + (rhov / rhovs) * rvs

;- done

   return, index
   end