function [CO2,H2O,PW,TOUT,RHO,varargout]=webb(co2,h2o,P,T,Tflag,inputunits)
% function [CO2,H2O,PW,TOUT,RHO]=webb(co2,h2o,P,T)
%
% March 2008: previously was airdry.m.
% modified to output the temperature and water vapor contributions to the Webb correction explicitly
% converts from LI7500 output
%
% INPUTS: co2 - mmol/m^3 wet air
%         h2o - mmol/m^3 wet air
%         P   - total pressure (KPa)
%         Ts  - sonic temperature (C or K)
%
% OUTPUTS:
%  CO2 - 3XN array
%      ROW 1: co2 in umol/mol dry air
%      ROW 2: co2 in umol/mol wet air
%      ROW 3: co2 in umol/m^3 wet air
%
%  H2O - 3XN array
%      ROW 1: h2o in mmol/mol dry air
%      ROW 2: h2o in mmol/mol wet air
%      ROW 3: H2O in mmol/m^3 wet air
%
%  PW - 1XN - partial pressure of water vapor (kPa)
%        1/15/2001 -  PW is output in KPa, not Pa !!!
%
%  RHO - 3XN array
%      ROW 1: mol dry air/m^3 wet air
%      ROW 2: mol wet air/m^3 wet air
%      ROW 3: Kg  moist air/m^3 moist air
%
%  TOUT - 2XN - temperatures
%      ROW 1: measured sonic temperature (C)
%      ROW 2: dried sonic temperature (K)
%
% Aug 26, 2006 - removed despiking - clean, continous signals
% should be sent to this script
%
% either dry or sonic temperature can be input - if it is dry
% temperature, need to have flag set to dry.  Sonic temperature is
% the default

if nargin<5
  Tflag='wet';
end

if nargin<6
  inputunits='density';
end

% - Convert temperature if sent in Celcius
if nanmedian(T)<100
  T=T+273.15;
end

% two scenarios considered
% 1. input co2/h2o units are density (mmol/m^3 moist air); ie an open path
% 2. input co2/h2o units are molar (mmol/mol moist air, umol/mol moist air)
%    and input temperature is dry temperature (ie a closed path)
%
% Need to generate PW, Td, total density, co2 molar and h2o molar.  Consider
% on the scenario 1 and 2 separately


if strcmp(inputunits,'density')
  % calculate the partial pressure of h2o vapor. Use the sonic
  % temperaure to approximate the dry air temperature.
  %
  % Pw (KPa) = (n/V) * R_u * T * 10e-6
  %
  %   - n/V   = LI7500 output in mmol/m^3 wet air
  %   - R_u   = 8.314 J/mol K
  %   - T     = dry air temp in K
  %   - 10e-6 = converts mmol to mol and Pa to KPa
  PW=1e-6*8.314*h2o.*T;

  if findstr(Tflag,'wet')
    % calculate dry air temperature from sonic temperature
    % using Gaynor eq:
    Td = T./(ones(size(T))+0.32*PW./P);

    % Make an iteration on the calculation of Pw, using the
    % dry air temperature
    PW = 1e-6*8.314*h2o.*Td;

    % recalculate dry temperature w/new pressure
    TD = T./(ones(size(T))+0.32*PW./P);
  elseif findstr(Tflag,'dry')
    TD=T;
  end

  % calculate wet air molar density (mol wet air / m^3 wet air)
  %  (n/V)_a = P/R_u/T = 1e3 / 8.314 * P /T
  rhomtotal = (1e3/8.314) * P./(TD);

  % calculate mol fraction of water vapor (mmol h2o/mol moist air)
  % in wet air
  h2owet = h2o./rhomtotal;

  % calculate mol fraction of co2 (umol co2/mol moist air)
  % in moist air
  co2wet = 1e3*co2./rhomtotal;

elseif strcmp(inputunits,'molar')

  h2owet=h2o;
  co2wet=co2;

  if findstr(Tflag,'wet')

    rhomtotal = (1e3/8.314) * P./(T);
    h2o=h2owet.*rhomtotal;

    PW=1e-6*8.314*h2o.*T;

    % calculate dry air temperature from sonic temperature
    % using Gaynor eq:
    Td = T./(ones(size(T))+0.32*PW./P);

    % Make an iteration on the calculation of Pw, using the
    % dry air temperature
    PW = 1e-6*8.314*h2o.*Td;

    % recalculate dry temperature w/new pressure
    TD = T./(ones(size(T))+0.32*PW./P);
  elseif findstr(Tflag,'dry')
    TD=T;
  end

  % calculate wet air molar density (mol wet air / m^3 wet air)
  %  (n/V)_a = P/R_u/T = 1e3 / 8.314 * P /T
  rhomtotal = (1e3/8.314) * P./(TD);
  h2o=h2owet.*rhomtotal;
  co2=co2wet.*rhomtotal/1e3;

  PW=1e-6*8.314*h2o.*TD;

else
  disp('webb.m scenario not recognized');
end
% calculate partial pressure of dry air - This is where the fortran
% code differs from mine.  That code assumes the open path irga
% outputs mols / m^3 of DRY air; however, i think (and Licor
% confirmed) that it is mols / m^3 of moist air, and so the partial
% pressure of dry air is the output of the irga minus the vapor pressure
Pa = P-PW;

% calculate dry air molar density (mol dry air / m^3 wet air)
%  (n/V)_a = Pa/R_u/T = 1e3 / 8.314 * Pa /T
rhomdry = (1e3/8.314) * Pa./(TD);
rhomwater = rhomtotal-rhomdry;
rhotot = rhomdry*29/1000+rhomwater*18/1000;

% calculate mol fraction of water vapor (mmol h2o/mol dry air) in dry air
h2odry = h2o./rhomdry;

% calculate mol fraction of co2 (umol co2/mol dry air) in dry air
co2dry = 1e3*co2./rhomdry;

% calculate the fluctuations caused by the webb terms (water vapor and
% temperature for co2 and h2o)

% calculate densities in grams/m^3 moist air to match Webb equations
Mc=44;
Ma=28.966;
Mv=18.016;

rhoa = rhomdry*Ma;
rhov = (rhomtotal-rhomdry)*Mv;
rhoc = co2*1000*44/10^6;

mu = Ma/Mv;
sigma =  mean(rhov)/mean(rhoa);

%co2
Xc_raw = detrend(rhoc,0)/mean(rhoa)*(10^6*Ma/Mc);

Xc_water_induced = mu*mean(rhoc)/mean(rhoa)^2*detrend(rhov,0)*(10^6*Ma/Mc);
Xc_T_induced = (1+mu*sigma)*mean(rhoc)/mean(rhoa)/mean(TD)*detrend(TD,0)*(10^6*Ma/Mc);
Xc_P_induced =-(1+mu*sigma)*mean(rhoc)/mean(rhoa)/mean(P)*detrend(P,0)*(10^6*Ma/Mc);

Xc_corr = Xc_raw+Xc_water_induced+Xc_T_induced+Xc_P_induced;

%water vapor
Xw_raw = detrend(rhov,0)/mean(rhoa)*(10^3*Ma/Mv);

Xw_water_induced = mu*sigma*detrend(rhov,0)/mean(rhoa)*(10^3*Ma/Mv);
Xw_T_induced = (1+mu*sigma)*mean(rhov)/mean(rhoa)/mean(TD)*detrend(TD,0)*(10^3*Ma/Mv);
Xw_P_induced =-(1+mu*sigma)*mean(rhov)/mean(rhoa)/mean(P)*detrend(P,0)*(10^3*Ma/Mv);

Xw_corr = Xw_raw+Xw_water_induced+Xw_T_induced+Xw_P_induced;

DWEBB = [Xc_raw(:) Xc_water_induced(:) Xc_T_induced(:) Xc_P_induced(:) Xw_raw(:) Xw_water_induced(:) Xw_T_induced(:) Xw_P_induced(:)]';
HWEBB = strvcat('xc_raw: raw co2 fluctuation (ppm) [xc_corr=xc_raw+xc_w+xc_T+xc_P]',...
  'xc_w: co2 fluctuation induced by h2o fluctuations (ppm)',...
  'xc_T: co2 fluctuation induced by T fluctuations (ppm)',...
  'xc_P: co2 fluctuation induced by P fluctuations (ppm)',...
  'xw_raw: raw h2o fluctuation (mmol/mol)',...
  'xw_w: h2o fluctuation induced by h2o fluctuations (mmol/mol)',...
  'xw_T: h2o fluctuation induced by T fluctuations (mmol/mol)',...
  'xw_P: h2o fluctuation induced by P fluctuations (mmol/mol)');
  
% output arrays
CO2 = [co2dry(:) co2wet(:) co2(:)*1000]';
H2O = [h2odry(:) h2owet(:) h2o(:)]';
RHO = [rhomdry(:) rhomtotal(:) rhotot(:)]';

varargout{1}=HWEBB;
varargout{2}=DWEBB;

if findstr(Tflag,'wet')
  TOUT = [T(:) TD(:)]';
elseif findstr(Tflag,'dry')
  TOUT = TD(:)';
end
return

