%  function [Fa] = calc_Fs(Te,Ta,To)
%
    function [Fa] = calc_Fa(Te,Ta,To)
 
 %Read parameters
 Param
 
%Radiative surface flux
 Frad = SIGMA*(To^4-Te^4-EPSA*Ta^4);
  
% Evaporation
 if To-Ta>DTcrit_conv
   Feva = BIGONE * (To-Ta-DTcrit_conv); 
 else
   Feva = 0;
 end

%Net surface heat flux
 Fs = Frad + Feva;