function do_all=ebm_1dm
%  	1-dimensional Budyko-Sellers type Energy Balance Model
%	
%	This is the same program as EBM-1D, 
%	except that it allows for variations of the solar constant	
%	to explore the stability of the Budyko-Sellers EBM
% from http://www.pik-potsdam.de/~claussen/lectures/physikalische_klimatologie/scripts/ebm-1dm.sci, 
% 
% Eli Tziperman: changed language to matlab, added sellers option.
% -----------------------------------------------------------------------
close all; clear all

budyko=1; % 1 for budyko, 0 for sellers
ebm(budyko);

budyko=0; % 1 for budyko, 0 for sellers
ebm(budyko)


function ebm(budyko)
%  	physical constants and boundary conditions:
%
if budyko==1; 
  budyko_sellers=' budyko'; 
else 
  budyko_sellers=' sellers'; 
end

a=204; 		% infrared cooling (lin. app.) W/m2
b=2.17;		% R = A + B*T  (W/m2/oC)
c=3.81;		% diffusivity for budyko model (W/m2/oC)
so00=1365;	% solar constant (W/m2)
luma=0.5;	% range of luminosity variations
lume=1.5;
albi=0.62;	% albedo of an ice-covered grid cell
tc1=-10;	% critical temperature of complete ice cover (oC)
tc2=0;		% ... and of no ice
tcc=tc2-tc1;
kappa=400.0e0;     % diffusion coeff. for sellers model
dt=0.001;

%	numerical constants:
%
ny=30;	% number of meridional grid cells (in one hemisphere)
% number of iterations:
if budyko
  itmax=50;
  max_bal_converged=1.e-4;
else
  itmax=5000;
  max_bal_converged=1.e-4;
end

kl=20;	% number of luminosity increments

% albedo of ice-free grid cells
clear alb0;
alb0=[1:ny]*0+0.25;

%	
%	grid
%
clear phi phid ds cds;
dp=pi/(2*ny);
phi=(0.5*dp):dp:(ny*dp);
phid=phi*180/pi;
for j=1:ny
  ds(j)=sin(phi(j)+dp/2)-sin(phi(j)-dp/2);
end
%
%	To compute the meridional distribution of insolation one has to compute
%	theratio of areas of a circular element and a spherical element
%	between latitudes p2 and p1
for j=1:ny
  p2=phi(j)+dp/2;
  p1=phi(j)-dp/2;
  s2=sin(p2);
  c2=cos(p2);
  s1=sin(p1);
  c1=cos(p1);
  cds(j)=(s2*c2-s1*c1+p2-p1)/(2*pi*(s2-s1));
end
%
%	Because of the tilt of the Earth´s axis, a simple cosine distribution 
%	of solar radiation with latitude does not provide a realistic estimate
%	in the annual mean. Hence a more appropriate empirical wheighting scheme
%	has been suggested by, for example, Henderson-Sellers and McGuffie:
%
%	cds=[0.30475 0.29725 0.28 0.25525 0.223 0.1925 0.156 0.13275 0.125];
%
%	initialize matrices
%
clear so tj alb;
clear albmm tmm tjmm;
clear diff insol ircool bal;
%
%	initial value of mean albedo
%
albm=0.;
for j=1:ny
  albm=albm+alb0(j)*ds(j);
end
%
% -------- luminosity loop ---------------------------------
%	
%	luminosity variations
%
clear lum tm;
dlum=(lume-luma)/kl;
lum=[luma:dlum:lume lume-dlum:-dlum:luma];
%lum=[luma:dlum:lume];

for k=1:(2*kl+1)
%for k=1:(kl+1)

  so0=so00*lum(k);	% new solar constant
  sq=so0/4;		% annual, global mean insolation
  so=so0*cds;		% meridional distribution of insolation

  %% initial values of mean temperature, and temperature
  %% distribution.  Specified these only for first luminocity value.
  %% Taken from previous luminocity runs for other cases.

  if k==1
    tm(k)=(sq*(1-albm)-a)/b;
    for j=1:ny
      tj(j)=tm(k);
    end
  else
    tm(k)=tm(k-1);
  end
  %
  % Iteration loop, solving for steady state temperature and ice fields:
  %
  alb=alb0;
  %% for i=1:it
  i=0; max_bal=1;
  while i<itmax && max_bal>max_bal_converged
    i=i+1;
    %	ice distribution: max ice if T<tc1, no ice if T>tc2, and
    %	linear interpolation in between:
    for j=1:ny
      if (tj(j) < tc1)
        ice(j)=1;
      elseif (tj(j) > tc2)
        ice(j)=0;
      else
        ice(j)=(tc2-tj(j))/tcc;
      end
      alb(j)=ice(j)*albi+(1-ice(j))*alb0(j);
    end
    %	mean albedo
    albm=0;
    for j=1:ny
      albm=albm+alb(j)*ds(j);
    end
    albmm(i)=albm;
    %	temperature distribution, mean temperature, and energy budget
    tjm=0;
    for j=1:ny
      if budyko
        tj(j)=(so(j)*(1-alb(j))+c*tm(k)-a)/(b+c);
        tjm=tjm+tj(j)*ds(j);
      elseif i>1
        tj(j)=tj(j)+dt*bal(j);
      end
      if budyko
        %% budyko "diffusion":
        diff(j)=c*(tm(k)-tj(j));
      else
        %% sellers diffusion:
        if j==1
          diff(j)=kappa*(tj(j+1)-tj(j));
        elseif j==ny
          diff(j)=-kappa*(tj(j)-tj(j-1));
        else
          diff(j)=kappa*(tj(j+1)-2*tj(j)+tj(j-1));
        end
      end

      insol(j)=so(j)*(1-alb(j));
      ircool(j)=-a-b*tj(j);
      bal(j)=insol(j)+diff(j)+ircool(j);
      max_bal=max(abs(bal));
    end
    tjmm(i)=tjm;
    tm(k)=tjm; 
  end

  tplot(k,:)=tj;
  iceplot(k,:)=ice;
  albedoplot(k,:)=alb;
  diffplot(k,:)=diff;
  insolplot(k,:)=insol;
  ircoolplot(k,:)=ircool;
  balplot(k,:)=bal;
  
  %--------- end of iteration loop ----------------------------
  %
end
%--------- end of luminosity loop -----------------------------------
%
%  plot results:

klum=1:(2*kl+1);
nlabels=4;
klum_ticks=1:round(((2*kl+1)-1)/nlabels):(2*kl+1);
latitude=phi*180/pi;

figure
subplot(2,2,1)
plot(klum,lum)
V=axis; axis([klum(1) klum(2*kl+1) V(3) V(4)]); % ([XMIN XMAX YMIN YMAX])
xlabel('klum')
ylabel('luminocity')
set(gca,'fontsize',13)

subplot(2,2,2)
[c,h]=contourf(klum, latitude, tplot');
set(gca,'Xtick',klum_ticks); set(gca,'XtickLabel',lum(klum_ticks))
for i=1:length(h); set(h(i),'LineStyle','none');end
colorbar; set(gca,'fontsize',13)
xlabel('luminocity')
ylabel('latitude')
title(['temperature' budyko_sellers])
colorbar

subplot(2,2,3)
[c,h]=contourf(klum, latitude, iceplot');
set(gca,'Xtick',klum_ticks); set(gca,'XtickLabel',lum(klum_ticks))
for i=1:length(h); set(h(i),'LineStyle','none');end
colorbar; set(gca,'fontsize',13)
xlabel('luminocity')
ylabel('latitude')
title('ice')
colorbar

subplot(2,2,4)
[c,h]=contourf(klum, latitude, albedoplot');
set(gca,'Xtick',klum_ticks); set(gca,'XtickLabel',lum(klum_ticks))
for i=1:length(h); set(h(i),'LineStyle','none');end
colorbar; set(gca,'fontsize',13)
xlabel('luminocity')
ylabel('latitude')
title('albedo')
colorbar

% plot term in temperature equation:

figure

subplot(2,2,1)
[c,h]=contourf(klum, latitude, balplot');
set(gca,'Xtick',klum_ticks); set(gca,'XtickLabel',lum(klum_ticks))
for i=1:length(h); set(h(i),'LineStyle','none');end
colorbar; set(gca,'fontsize',13)
xlabel('luminocity')
ylabel('latitude')
title('bal')
colorbar

subplot(2,2,2)
[c,h]=contourf(klum, latitude, ircoolplot');
set(gca,'Xtick',klum_ticks); set(gca,'XtickLabel',lum(klum_ticks))
for i=1:length(h); set(h(i),'LineStyle','none');end
colorbar; set(gca,'fontsize',13)
xlabel('luminocity')
ylabel('latitude')
title('ircool')
colorbar

subplot(2,2,3)
[c,h]=contourf(klum, latitude, diffplot');
set(gca,'Xtick',klum_ticks); set(gca,'XtickLabel',lum(klum_ticks))
for i=1:length(h); set(h(i),'LineStyle','none');end
colorbar; set(gca,'fontsize',13)
xlabel('luminocity')
ylabel('latitude')
title('diff')
colorbar

subplot(2,2,4)
[c,h]=contourf(klum, latitude, insolplot');
set(gca,'Xtick',klum_ticks); set(gca,'XtickLabel',lum(klum_ticks))
for i=1:length(h); set(h(i),'LineStyle','none');end
colorbar; set(gca,'fontsize',13)
xlabel('luminocity')
ylabel('latitude')
title('insol')
colorbar

figure(3)
if budyko
  subplot(2,2,1)
  plot(latitude,tplot(35,:),'r');
  ylabel('temperature')
  title('budyko=r, sellers=b')
  hold on
  subplot(2,2,2)
  plot(latitude,iceplot(35,:),'r');  
  ylabel('ice')
  title('budyko=r, sellers=b')
  hold on
else
  subplot(2,2,1)
  plot(latitude,tplot(28,:),'b');  
  subplot(2,2,2)
  plot(latitude,iceplot(28,:),'b');  
end
xlabel('latitude')
