  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

