% Broadband radiative model: Gupta, NASA
% reguire sundist.m function
% input:
%       zen - solar zenith angle in deg
%       doy - Julian day
%       PW - Total water contant[g/cm^2]
%       DU - Total ozon in Dobson [DU]
%       AOT - Broadband Aerosol optical thickness
%       SSA - Broadnand Single scattering Albedo
%       g - Broadband asymetry parameter 
%       p - pressure in [hpa]
%       As - surface albedo
% output :
%         Ftot - total downward flux [Wm^-2] at the surface
%         Fdir - direct flux [Wm^-2]
%         Fdif - diffuse flux
clear
% input parameters
PW=2;     
DU=300;
p=1000;
aot=0.25;
omega=0.95;
g=0.75;
zen=30;
doy=200;
As=0.05;  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  Efective Attanuation  
aH2O=0.10*PW.^0.27;            
aO3=0.037*(DU/1000).^0.43;
aCO2=0.006*(p/1013*350/300).^0.29;
aO2=0.0075*(p/1013).^0.87;
aRay=0.035*(p/1013).^0.67;
aaer=aot.*(1-omega)+0.5*aot.*omega.*(1-g);
a=aH2O+aO3+aCO2+aO2+aRay+aaer;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tau=-log(1-a);
N=1.1-2*tau;
tauz=tau.*(1./cos(zen*pi/180)).^N;
B=0.065.*(p/1013).*As+2*As.*aot.*omega.*(1-g);

Ta=(1+B).*exp(-tauz);   
FTOA=1368*sundist(doy)*cos(zen*pi/180);   % TOA sun flux  
F=FTOA.*Ta;                                % total surface flux
TAU=-log(1-(aH2O+aO3+aCO2+aO2+aRay));
TAUZ=TAU.*(1./cos(zen*pi/180)).^N+aot./cos(zen*pi/180);
Tdir=exp(-TAUZ);
Fdir=1368*sundist(doy)*cos(zen*pi/180).*Tdir
Fdif=F-Fdir
Ftot=F