% Solve radiative tranfer model (Schawrtzwil eqation) for longwave range
clear
close all
Kabs=1e2;
emiss=1;
lambda=12;
mi=1;
% Read atmospheric profil
%x=load('jmls.lay');
x=load('prof.dat');
Z=x(:,1);
P=x(:,2);
T=x(:,3);
WV=x(:,4);

% 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);
%t(1:2)=[255,270];  % surface inversion
%t=t*0+290;
Tsurf=t(1);
dens=100*p./(287.*t);
densWV=dens.*q*0.622/1000;
densWV=densWV*0+0.001;
nrec=length(p);


% compute profile of optical deph
tau(nrec)=0;
for k=1:nrec-1
 tau(k)=trapz(z(k:end),Kabs*densWV(k:end));
end

% Find upward radiation
[B]=planck(Tsurf,lambda)/pi;
Bsurf=emiss.*B;
I(1)=Bsurf;
Isurf=Bsurf.*exp((tau-tau(1))/mi);
[B]=planck(t,lambda)/pi;  
for k=2:nrec
  term=B/mi.*exp((tau-tau(1))/mi);
  Ia(k)=trapz(-tau(1:k),term(1:k));
end  
Iup=Isurf+Ia;
disp(I(end))
temp_radiacyjna(Iup(end),lambda)

% compute the downward radiation
Idn(nrec)=0;
for k=nrec-1:-1:1
%  [B]=planck(t,lambda)/pi;  
  term=1/mi*B.*exp((tau-tau(1))/mi);
  Idn(k)=trapz(-tau(k:nrec),term(k:nrec));
end  
disp(Idn(1))
temp_radiacyjna(Idn(1),lambda)
[z',tau',Idn',Iup']

