% Solve radiative tranfer model (Schawrtzwil eqation) for longwave range in gray approximation
% Absorption coefficient is constant with wavelength
% TODO - improve the high TAU case - problem with "z" resolution
clear
close all
Kwv=0.5e2;
Kco2=1e2;
Ko3=1e1;
%Kwv=0;
%Kco2=0;
emiss=1;     % Earth's surface emissivity
sig=5.67e-8;  % Stefan Boltzmiann
% Read atmospheric profil
%x=load('jmls.lay');
x=load('prof.dat');
Z=x(:,1);
P=x(:,2);
T=x(:,3);
WV=x(:,4);
O3=x(:,5);

% Define a new vertical levels
z=[0:0.1:100];
z=Z';
z=[0:0.25:10,15:5:100];
t=interp1(Z,T,z);
p=interp1(Z,P,z);
q=interp1(Z,WV,z);
Xo3=interp1(Z,O3,z);
Xco2=410*1e-6;    %ppm
Mco2=44.01;
Mair=29.97;
Mo3=48.0;
%t(1:4)=[215,240,270,280];  % surface inversion
%q=mean(q)+q*0;
%t=t*0+290;
Tsurf=t(1);
dens=100*p./(287.*t);
densWV=dens.*q*0.622/1000;
densCO2=dens.*Xco2*Mco2/Mair;
densO3=dens.*Xo3*Mo3/Mair;

% Include cloud
ztop=10; 
dz=1;
COD=1;
ind=find(z>=ztop-dz&z<=ztop);
ABScloud=zeros(size(z));
ABScloud(ind)=COD/dz;
nrec=length(p);
ABScloud=ABScloud*0;

% compute profile of optical deph
tau(nrec)=0;
for k=1:nrec-1
 ABS=Kwv*densWV(k:end)+Kco2*densCO2(k:end)+Ko3*densO3(k:end)+ABScloud(k:end);
 tau(k)=trapz(z(k:end),ABS);
end

% Find downward flux
Fdn(nrec)=0;
for k=nrec-1:-1:1
  dtau=abs(tau(k)-tau(k:nrec));
  Eterm=Efun(dtau,2);
  Fdn(k)=2*trapz(-tau(k:nrec),Eterm.*sig.*t(k:nrec).^4);
end  

% Find upward flux
Fsurf(1)=emiss*sig*t(1).^4;
for k=2:1:nrec
  dtau=abs(tau(1)-tau(k));
  Eterm3=Efun(dtau,3);
  Fsurf(k)=2*Fsurf(1).*Eterm3;
  dtau=abs(tau(k)-tau(1:k));
  Eterm=Efun(dtau,2);
  Fa(k)=2*trapz(-tau(1:k),Eterm.*sig.*t(1:k).^4);
end  
Fup=Fsurf+Fa;

[z',t',tau',Fdn',Fup']
plot(Fdn,z,Fup,z)

NET=Fdn-Fup;



