% Isothermal atmosphere in radiative equilibrium
% NB: assumed transparent to solar radiation

% ** Input:
% 
% TAUinf = optical depth of the atmosphere (no dim)
%   NB: this is proportional to the greenhouse gas concentration
%   times their absorption coeff 
%
% ** Output
% 
% Ta = temperature of atmosphere (K)
% Ts = surface temperature (K)
% Ze = emission level (m)
% Iup(z) = upward radiation (W m-2)
% Idown(z) = downward radiation (W m-2)
%
%(c) A. Czaja, Dec. 2005

%Planetary constants
GRAVI = 9.81;
ALPHA = 0.3, %planetary albedo (no dim)
So = 1367; % Solar constant (W m-2)
SIGMA = 5.67*1.e-08; %Stefan constant (W m-2 K-4)
TE=( (1-ALPHA)*So/(4*SIGMA) )^.25, %emission temperature (K)

%Greenhouse gas constituants: optical depth of the atmosphere
%TAUinf = 0.5; 
TAUinf = input('Optical depth?  ');

%Color for plot
RGBK = input('Plotting color?   '); %answer should be 'r' (red), 'g' (green), 'b' (blue) ...
%Vertical grid
km = 1000;
z = [0:.01:100]*km;

%Isothermal atmosphere
Rd = 287; %ideal gas constant for dry air (J kg-1 K-1)
Ta = TE/(1+exp(-TAUinf))^.25,
H = Rd * Ta / GRAVI, %scale height (m)

%Surface temperature
Ts = Ta*2^.25,

%Optical depth
TAU = TAUinf * (1-exp(-z/H));

%Upward radiation
Iup = SIGMA*exp(-TAU).*(Ts^4 + (exp(TAU)-1)*Ta^4);
Iup_inf = SIGMA*exp(-TAUinf).*(Ts^4 + (exp(TAUinf)-1)*Ta^4);

%Downward radiation
Idown = Iup-Iup_inf;

% Weighting function
W = abs((TAU-TAUinf).*exp(TAU-TAUinf));
W = W/max(W); %normalization

% Emission level
ZE = H*log(TAUinf);

% Figures
figure(1), hold on,
[h]=plot([0 1],[ZE/km ZE/km],[RGBK '--'],W,z/km,RGBK); hold on
set(h,'Linewidth',2);
set(gca,'Fontsize',15);
title('Weighting function as a function of height','Fontsize',30);
xlabel('Weighting function','Fontsize',25);
ylabel('Height (km)','Fontsize',25);
axis([0 1 0 50]);

figure(2), clf,
[h]=plot(Iup,z/km,'b-',[0 500],[ZE/km ZE/km],'b--'); hold on
set(h,'Linewidth',2);
[h]=plot(Idown,z/km,'r-',[0 500],[ZE/km ZE/km],'b--'); hold on
set(h,'Linewidth',2);
set(gca,'Fontsize',15);
xlabel('Radiation (W m-2)','Fontsize',25);
ylabel('Height (km)','Fontsize',25);
title('Upward (blue) and downward (red) radiation','Fontsize',20);
axis([0 500 0 50]);

figure(3), hold on, 
[h]=plot(Iup,z/km,RGBK,[0 500],[ZE/km ZE/km],[RGBK '--']); hold on
set(h,'Linewidth',2);
%[h]=plot([2*SIGMA*TE^4 2*SIGMA*TE^4],[0 500],'k--'); hold on
%set(h,'Linewidth',2);
%[h]=plot([SIGMA*TE^4 SIGMA*TE^4],[0 500],'g--'); hold on
%set(h,'Linewidth',2);
set(gca,'Fontsize',15);
xlabel('Radiation (W m-2)','Fontsize',25);
ylabel('Height (km)','Fontsize',25);
title('Upward radiation as a function of height','Fontsize',30);
axis([0 500 0 40]);


