%  function [Fs] = calc_Fs(Te,Ta,Ts,EPSa)
%
    function [Fs,Feva] = calc_Fs(Te,Ta,Ts,EPSa)
 
 %Read parameters
 TwoBox_Param
 
 %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;