%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Relationship between H and tropical surface temperature:
%
% H = energy export from moist to dry tropical regions (in W)
% Ts = Average surface temperature of the Tropics (in K)
% DTs = Ts1-Ts2 = surface temperature diff. between moist and dry regions (in K)
%
% (c) A. Czaja, April 2006
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%-------------------------------------------------------------------------
% Main model parameter (denoted by a in the notes)
FRAC = input('Fraction of Tropics occupied by moist regions? ');
%-------------------------------------------------------------------------

% Planetary constants
RADIUS = 6371*1000;
PW = 1.e15;
SIGMA = 5.67*1.e-08;
So = 1367; 
ALPHA = 0.3;
ATOT = 2*pi*RADIUS^2; %area of 30S-30N
S = (1-ALPHA)*So*(RADIUS^2)*(pi/3+sqrt(3)/2); %energy exchanged between 30S-30N and Space

% Model parameters
A1 = FRAC*ATOT; %area of moist regions (in m2)
A2 = ATOT-A1; %area of dry regions (in m2)
GAMMA = 6.5; %lapse rate (in K / km)
Z1 = 10; %emission height of moist regions (in km)

% Solve for T1 and Ts2 as a function of H
H = [0:.01:20]*PW; %energy export from moist to dry regions
N = length(H);

for n = 1:N
   T1(n) = (S/(SIGMA*ATOT)-H(n)/(A1*SIGMA))^.25;
   Ts2(n)= (S/(SIGMA*ATOT)+H(n)/(A2*SIGMA))^.25;
end

% Compute DTS
Ts1 = T1 + GAMMA*Z1;
DTs = Ts1-Ts2;

% Mean surface temperature
Ts = (A1*Ts1+A2*Ts2)/ATOT;

% Approximate solution for Ts (valid if FRAC~=1/2 and H/S <<1)
TE=(S/(SIGMA*ATOT))^.25; %emission temperature (in K)
Ts_approx = TE+FRAC*GAMMA*Z1-3*TE*(H/S).^2/(32*FRAC*(1-FRAC));

% Plot Ts1-Ts2 as a function of H
figure(1), clf,
h=plot(H/PW,DTs,'b-'); grid on;
set(h,'Linewidth',2);
set(gca,'Fontsize',15);
xlabel('Energy export from moist to dry regions (PW)','Fontsize',15);
ylabel('Surface temperature difference (K)','Fontsize',15);
title('Ts1-Ts2 in K','Fontsize',20);
axis([0 20 -20 80]);

% Plot Ts, Ts1, Ts2 as a function of H
figure(2), clf,
h=plot(H/PW,Ts,'g-',H/PW,Ts_approx,'g--'); grid on; hold on
h=plot(H/PW,Ts1,'r-',H/PW,Ts2,'b-'); grid on; hold on
set(h,'Linewidth',2);
set(gca,'Fontsize',15);
xlabel('Energy export from moist to dry regions (PW)','Fontsize',15);
ylabel('Surface temperature (K)','Fontsize',15);
title('Ts (green), Ts1 (red) & Ts2 (blue) in K','Fontsize',20);
axis([0 20 280 320]);
