function  rhoa=air_dens(Ta,RH,Pa)
% AIR_DENS: computes computes the density of moist air.
% rhoa=AIR_DENS(Ta,RH,Pa) computes the density of moist air.
% Air pressure is optional.
%
% INPUT:   Ta  - air temperature Ta  [C]
%          RH  - relative humidity  [%]
%          Pa  - air pressure (optional) [mb] 
%
% OUTPUT:  rhoa - air density  [kg/m^3]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4/7/99: version 1.2 (contributed by AA)
% 8/5/99: version 2.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% load constants
as_consts;

if nargin == 2
  Pa = P_default; 
end

o61  = 1/eps_air-1;                % 0.61 (moisture correction for temp.)
Q    = (0.01.*RH).*qsat(Ta,Pa);    % specific humidity of air [kg/kg]
T    = Ta+CtoK;                    % convert to K
Tv   = T.*(1 + o61*Q);             % air virtual temperature
rhoa = (100*Pa)./(gas_const_R*Tv); % air density [kg/m^3]