%  function [Fs] = calc_Fs(Te,Ta,To)
%
    function [Fs,Feva] = calc_Fs(Te,Ta,Ts)
 
 %Read parameters
 Param
 
 %Compute emissivity of atmospheric column
  EPSA = calc_EPSA(Ta);

 %Radiative surface flux
 Frad = SIGMA*(Ts^4-Te^4-EPSA*Ta^4);
  
% Evaporation
 if Ts-Ta>DTcrit_conv
   Feva = BIGONE * (Ts-Ta-DTcrit_conv); 
 else
   Feva = 0;
 end

%Net surface heat flux
 Fs = Frad + Feva;