function [licorval,h2ommol] = licorWebb(co2air,h2oair,pair,tair,RH,tvirt)
% pair=PairL';tair=TairL';RH=RHL';
% Calculate the concentration of CO2 using the raw licor data
% from the 7500.  Uses an approach that is similar to the closed
% path sensors.  
%
% The inpu is:
%
% co2air: the number density of CO2 in moist air (mmol/m^3); this is rho_c in
%         Webb et al (1980) and c in McGillis et al. (2001)
% h2oair: the number density of H2O in moist air (mmol/m^3)
% tair: the air temperature in C from the vaisala or other reference T
% RH: the relative humidity in % from the vaisala or other reference RH
% pair: the air pressure in mb from the vaisala or fast response pressure
% tsonic: the sonic virtual temperature
%
% fs is the sonic sampling frequency
%
%  NOTE: This program assumes that all variables have a common time base,
%  i.e., they are all "sampled" at fs Hz using, e.g., interp1.m
%
% This version provides the option to correct for hydrostatic 
% pressure fluctuations on moving platforms if a fast response pressure 
% signal is available.  For example, the SO GasEx analysis used 
% mean Vaisala P + dp = P - rho g dz where dz is taken from the heave
% measurement.

S1=length(tvirt);
t1=(1:S1)/S1;
V1=length(tair);
tV=(1:V1)./V1;
Tair=interp1(tV,tair,t1,'spline','extrap');
%Pair=interp1(tV,pair,t1,'spline','extrap');
Pair=pair;
Exx=6.1121*exp(17.502*Tair./(240.97+Tair));
Exx=Exx.*(1.0007+Pair*3.46E-6);
Qsat=.622*Exx./(Pair-.375*Exx)*1000;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Determine the median values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rh = median(RH)/100;                       %Relative humidity as a fraction 
P0 = mean(pair);                           %Atmospheric pressure (mb)
Tair = median(tair);                       %Air temperature (C)     
meanP=P0/10;                               %Atmospheric pressure (kPa)
flucP=pair/10;                             %if Fast response pressure (kPa)     
Press=P0*100;                              %Atmospheric pressure (Pa)
Tairk=Tair+273.16;                         %Air temperature (K)
Exx=6.1121*exp(17.502*Tair./(240.97+Tair));
Esat=Exx.*(1.0007+P0*3.46E-6);             %Saturation vapor pressure
Ems=Esat*rh;                               %Vapor pressure
Qms=.622*Ems/(P0-.378*Ems)*1000;           %Specific humidity (g/Kg)
Tvirt=Tairk*(1+0.61*Qms/1000);             %Virtual temperature
rhoair = (Press-Ems*100)/287/Tairk;        %Density of dry air
h2ommol=29000*meanP*Qms/18/8.314/Tvirt;    %mmol/m3
Rms = Qms./(1-Qms/1000);                   %Mixing ratio g/Kg

% The following equation computes the mole fractions (or mass concentrations) using 
% expressions given by the Licor 7500 manual pages 2-9 and 2-10.  The problem with 
% these expressions is that they are appropriate for number density computed in dry 
% air, while the number densities output by the Licor are the moles per unit volume
% of moist air.  
%
%For future reference, note that the following steps are equivalent
%
%q = 18*h2oair*8.314/29*Tvirt/meanP/1000;      %Specific Humidity g/Kg
%
%or
%
%wcell = h2oair*8.314*Tairk/meanP/1000;     %H2O mass concentration mmol/mol
%M = 29 - 11/1000*wcell;                    %Molecular weight of moist air
%q = 18*wcell./M;                           %Specific humidity g/Kg
%
%M*q = 18*wcell;
%29*q - 11/1000*wcell*q = 18*wcell
%29*q = 11/1000*wcell*q + 18*wcell
%29*q = 18*wcell*(11/18*q/1000 + 1)
%29*q = 18*wcell*(1 + 0.61*q/1000)
%29*q = 18*h2oair*8.314*Tairk/meanP/1000*(1 + 0.61*q/1000)
%29*q = 18*h2oair*8.314/meanP/1000*Tvirt
%q = 18*h2oair*8.314/meanP/1000*Tvirt/29
%
%
%Now let's check
%
%co2 = (44*8.314/29/meanP)*co2air.*Tvirt;  % Specific content of CO2
%
%versus
%
%ccell = co2air*8.314*Tairk/meanP;     %CO2 mass concentration mmol/mol
%co2 = 44*ccell./M;                         % Specific content of CO2
%co2 = 44*co2air*8.314*Tairk/meanP/M;
%co2 = 44*co2air*8.314*Tairk/meanP/29/(1 - 11/29/1000*wcell);
%co2 = 44*co2air*8.314*Tairk/meanP/29/(1 - 11/29/1000*q*M/18);
%co2 = 44*co2air*8.314*Tairk/meanP/29/(1 - 11/29/1000*q/18*(29-11*wcell/1000));
%co2 = 44*co2air*8.314*Tairk/meanP/29/(1 - 11/18*q/1000+11/29/1000*q/18*11*wcell/1000));
%co2 ~ 44*co2air*8.314*Tairk/meanP/29/(1 - 0.61*q/1000));
%co2 ~ 44*co2air*8.314*Tairk/meanP/29*(1 + 0.61*q/1000));
%co2 ~ (44*8.314/29/meanP)*co2air.*Tvirt;  % Specific content of CO2
%
% Convert from number to mass density
h2oair=18*h2oair; %g/m3
co2air=44*co2air; %mg/m3

% Get specific Humidity g/Kg from Licor using Tvirt from sonic+Vaisala

qair = (8.314/29/1000)*h2oair.*tvirt./flucP;   %g/Kg
tair = tvirt./(1+0.61*qair./1000);             %Fluctuating air temperature
qco2 = (8.314/29)*co2air.*tvirt./flucP;        %mg/Kg

%***********************************************
% Webb correct 
%***********************************************

rair = qair./(1-qair/1000);                    %Mixing ratio g/Kg
rco2 = qco2./(1-qair/1000);                    %Mixing ratio mg/Kg
% rair = qair./(1-1.6*qair/1000);                    %Mixing ratio g/Kg
% rco2 = qco2./(1-1.6*qair/1000);                    %Mixing ratio mg/Kg

rhfast=qair./Qsat;                             %fast rh from Licor and Vaisala
% NOTE:  this version outputs the mixing ratio

licorval = [rco2; rair; tvirt; tair; qair; rhfast];
