function q=qsat(Ta,Pa)
% QSAT: computes specific humidity at saturation. 
% q=QSAT(Ta) computes the specific humidity (kg/kg) at satuation at
% air temperature Ta (deg C). Dependence on air pressure, Pa, is small,
% but is included as an optional input.
%
%    INPUT:   Ta - air temperature  [C]
%             Pa - (optional) pressure [mb]
%
%    OUTPUT:  q  - saturation specific humidity  [kg/kg]

% Version 1.0 used Tetens' formula for saturation vapor pressure 
% from Buck (1981), J. App. Meteor., 1527-1532.  This version 
% follows the saturation specific humidity computation in the COARE
% Fortran code v2.5b.  This results in an increase of ~5% in 
% latent heat flux compared to the calculation with version 1.0.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3/8/97: version 1.0
% 4/7/99: version 1.2 (revised as above by AA)
% 8/5/99: version 2.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin==1,
  as_consts;
  Pa=P_default; % pressure in mb
end;

% original code
% a=(1.004.*6.112*0.6220)./Pa;
% q=a.*exp((17.502.*Ta)./(240.97+Ta))

% as in Fortran code v2.5b for COARE
ew = 6.1121*(1.0007+3.46e-6*Pa).*exp((17.502*Ta)./(240.97+Ta)); % in mb
q  = 0.62197*(ew./(Pa-0.378*ew));                         % mb -> kg/kg