 function [qsat, qsati] = t_to_qsat(t,p)
% t_to_qsat generates saturated specific humidity (g/kg) wrt water and ice
% given t (oC), and p (mb)
% [qsat, qsati] = t_to_qsat(t,p)

%convert to SI units
 tk = t + 273.15; %K
 p = p*100; %Pa

%calculate saturated specific humidities
 eps = 0.622; 	%Rd/Rv
 [es esi] = thermo_es(t);   	%es and es_wrt_ice in mb
 es = es*100; esi = esi*100;    % convert to Pa
 rsat = eps*es./p;		 	%saturation mixing ratio
 rsati = eps*esi./p;		 	%sat mixing ration wrt ice
 qsat = rsat./(rsat+1);			%sat specific humidity
 qsati = rsati./(rsati+1);		%sat specific humidity wrt ice

%convert to g/kg
 qsat = 1000*qsat;
 qsati =1000*qsati; 

