%Codes a simple high latitude climate model (energy balance)that includes 
%the following physical processes:
%
%  (1) SW Radiation
%  (2) LW Radiation
%  (3) Ocean Heat Transport
%  (4) Atm. Heat Transport
%
% Among the processes NOT included in the model are:
%
%  (1) Convection
%  (2) Seasonal Cycle
%
% (dsa, 3/13/07)
function doall=energy_balance_high_lat()

%default climate parameters

%latitude range of interest
lat1=60;
lat2=90;

%see Trenberth and Caron (2001)
TRo0=0.5;
TRa0=3.5;

%see Harrison et al. (1990)
%remember, albedo is weighted toward summer value (when insolation is higher)
alpha0=0.4;

%choose epsilon to give a reasonable average temperature (~-5 to -10)
epsilon0=0.6;
  

%plotting options
MarkSize=8;
FontSize=14;
LineWidth=2;
TX=0.015;
TY=1.09;



[Ts0,Ta0]=energy_balance(TRo0,TRa0,alpha0,epsilon0,lat1,lat2);


epsilon=0.4:0.025:1;
[Ts,Ta]=energy_balance(TRo0,TRa0,alpha0,epsilon,lat1,lat2);
subplot(2,2,1)
plot(epsilon,Ts,'b',epsilon0,Ts0,'ro','MarkerSize',MarkSize,'LineWidth',LineWidth)
xlabel('Atm. Emissivity []','FontSize',FontSize)
set(gca,'FontSize',FontSize)
%text(TX,TY,'(a) Atm. Emissivity [ ]','Units','Normalized','FontSize',14)
ylabel('Surf. Temp. [^\circC]')

alpha=0.2:0.05:0.6;
[Ts,Ta]=energy_balance(TRo0,TRa0,alpha,epsilon0,lat1,lat2);
subplot(2,2,2)
plot(alpha,Ts,'b',alpha0,Ts0,'ro','MarkerSize',MarkSize,'LineWidth',LineWidth)
xlabel('Albedo []','FontSize',FontSize)
set(gca,'FontSize',FontSize)
%text(TX,TY,'(b) Albedo [ ]','Units','Normalized','FontSize',14)


TRo=0:0.25:2;
[Ts,Ta]=energy_balance(TRo,TRa0,alpha0,epsilon0,lat1,lat2);
subplot(2,2,3)
plot(TRo,Ts,'b',TRo0,Ts0,'ro','MarkerSize',MarkSize,'LineWidth',LineWidth)
xlabel('Ocean Heat Transport [PW]','FontSize',FontSize)
set(gca,'FontSize',FontSize)
%text(TX,TY,'(c) Ocean Heat Transport [PW]','Units','Normalized','FontSize',14)
ylabel('Surf. Temp. [^\circC]')



TRa=2:0.25:7;
[Ts,Ta]=energy_balance(TRo0,TRa,alpha0,epsilon0,lat1,lat2);
subplot(2,2,4)
plot(TRa,Ts,'b',TRa0,Ts0,'ro','MarkerSize',MarkSize,'LineWidth',LineWidth)
xlabel('Atm. Heat Transport [PW]','FontSize',FontSize)
set(gca,'FontSize',FontSize)
%text(TX,TY,'(c) Atm. Heat Transport [PW]','Units','Normalized','FontSize',14)


set(1,'Position',[560 200 600 600],'PaperPositionMode','auto')
print -djpeg -f1 energy_balance.jpg


%solve simple high latitude energy balance model
function [Ts,Ta]=energy_balance(TRo,TRa,alpha,epsilon,lat1,lat2)

R=6.37e6;
A=2*pi*R^2*(sin(lat2*pi/180)-sin(lat1*pi/180));
sigma=5.67e-8;
S0=avg_solar(lat1,lat2);

Fa=TRa*1e15/A;
Fo=TRo*1e15/A;

Ts=((Fo+0.5*Fa+(1-alpha)*S0)./((1-epsilon/2)*sigma)).^0.25;
Ta=((epsilon*Fo+Fa+(1-alpha)*S0*epsilon)./(2*(1-epsilon/2).*epsilon*sigma)).^0.25;

Ts=Ts-273.15;
Ta=Ta-273.15;

%North (1975) - approximate annual average incoming radiation to within
%2 percent at all latitudes
function S=solar(lat)
  
lat=lat*pi/180;  
  
Q0=1365;  
Q2=-0.482;  

S=(Q0/4)*(1+(Q2/2)*(3*sin(lat).^2-1));  


%analytical area-weighted average of the North (1975) function
function S_bar=avg_solar(lat1,lat2)
  
lat1=lat1*pi/180;
lat2=lat2*pi/180;
  
Q0=1365;  
Q2=-0.482;  

S_bar=(Q0/4)*(1-Q2/2+Q2/2*((sin(lat2).^3-sin(lat1).^3)./(sin(lat2)-sin(lat1))));
