clear;
%addpath /home/phuybers/Ice/Inso;

initialize;

if Var.qplot==1,
  figure(1); clf; hold on;
end;
ct=0;
while State.year<=Var.t_end;

  %-----------------------------------------------------------
  %EBM 
  I=get_inso(State.year,Var.phi2); 
  for iter=1:3,
    hh=State.hh; 
    [State]=EBM(State,Var,I);
    M(iter,:)=State.hh-hh-State.Flow;
  end;
  State.M=M(end,:); 
    
  %-------------------------------------------------------
  %Icesheet
  [State]=ice(State,Var);       

  %-------------------------------------------------------
  %Display and save output
  if Var.qplot,
    subplot(3,1,[2 3]);  cla;  hold on;
    h(1)=plot(Var.phi2,State.M*100,'r'); 
    h(1)=plot(Var.phi2,State.Flow*100,'b'); 
    h(2)=plot(Var.phi2,State.hh+State.H,'k'); 
    h(3)=plot(Var.phi2,State.H,'k');     
    axis([Var.phi2(1) Var.phi2(end) min([-1000 State.H]) max([2000 State.H+State.hh])]); 
    subplot(311); hold on; 
    plot(State.year/1000,sum(State.hh),'k.');
    drawnow;
  end;
  ct=ct+1;
  if rem(State.year,1000)==0, 
    yearall(ct)=State.year;
    hhall(ct,:)=State.hh;
    Hall(ct,:)=State.H; 
    display([State.year/1000 sum(hhall(ct,:))]);
    [dummy StateAll(ct)]=EBM(State,Var,I);
    if rem(ct,100)==0,
      eval(['save run_epsa',num2str(Var.eps_a),'_uo',num2str(Var.uo),'_hs',num2str(Var.hs),'_Oct15.mat;']);
    end;
  end;
end;