%
%  [ESAT,QSAT] = saturation_humidity(TS,PS)
%
% Computes saturation water vapor pressure (ESAT)
% and saturation specific humidity (QSAT in kg/kg)  
% at a given temperature TS(Kelvin) and (total) pressure PS(mb)
%
% refs: Physics package of the MIT GCM
%
% !!! WHEN TS IS A VECTOR ... BUG SO NO
%     VECTORISATION ALLOWED
%
        function [ESAT,QSAT] = Clausius_Clapeyron(TS,PS)
%
          TS=max(TS,35.0);
          TC=TS-273.15;
          DENOM=243.5+TC;
          EPS = 0.62197; %RD / RV

          if(TC >= 0)  
%
% Satu. vap pressure for liquid water (mb)
           ESAT = 6.112*exp(17.67*TC/DENOM);
%
          else
%
% Satu. vap. pressure for ice (mb)
           ESAT = exp(23.33086-6111.72784/TS+0.15215*log(TS));

          end
  
% Saturation specific humidity at TS, PS
%
          QSAT = EPS*ESAT/(PS-ESAT*(1-EPS));

return