%
%  [ESAT,QSAT] = Clausius_Clapeyron(TS,PS)
%
% Computes saturation water vapor pressure (ESAT in mb)
% and saturation specific humidity (QSAT in kg/kg)  
% at a given temperature TS(in K) and pressure PS(in 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

