function [State,StateAll]=EBM(State,Var,I);

%An EBM coupled designed to be coupled to an icesheet
%function [State,StateAll]=EBM(State,Var,I);

if nargout==2,
  %Initialize daily totals
  Daily.Ti=zeros(1,Var.num); 
  Daily.Ts=zeros(1,Var.num); 
  Daily.Ta=zeros(1,Var.num); 
  Daily.Da=zeros(1,Var.num); 
  Daily.Is=zeros(1,Var.num); 
  Daily.Ia=zeros(1,Var.num); 
  Daily.alpha=zeros(1,Var.num); 
  Daily.Acc=zeros(1,Var.num); 
  Daily.Abl=zeros(1,Var.num); 
  Daily.count=0;
end;

%==================================================
%Forward integation
for t=0:Var.dt:365*24*3600-Var.dt;

  %------------------------------
  %Eddy heat flux
  small_steps=4*24; %15 minute time-step
  Da=zeros(1,Var.num); 
  for iter=1:small_steps,
    dTa=diff(State.Ta);
    Fmerid=1000*abs(dTa).*dTa;
    temp=diff([0 Fmerid 0]); 
    State.Ta=State.Ta+Var.dt/small_steps*temp/Var.Ca;
    Da=Da+temp; 
  end;
  Da=Da/small_steps;
  
  %--------------------------------
  %albedo 
  alpha(State.hh>0)=Var.alpha_i;  %ice
  alpha(State.hh<=0)=Var.alpha_l; %land 

  %-------------------------------
  %Shortwave flux  
  S=I(:,floor(t/(3600*24))+1)'; 
  Ss=S.*Var.T.*(1-alpha)./(1-alpha.*Var.R); 
  Sa=S.*Var.A+S.*Var.T.*alpha*Var.A./(1-alpha.*Var.R); 
  
  %------------------------------
  %Long-wave radiation
   %Hm-hh is distance from surface to middle atmosphere
   Tbl=State.Ta+Var.lr*(Var.Hm-(State.hh+State.H)/1000);
   Tupper=State.Ta-Var.lr*Var.Htm;
   %fluxes from atmosphere and surface
   Ia=Var.sigma_r.*( State.Ts.^4 - Var.eps_a*Tbl.^4 - Var.eps_a*Tupper.^4);
   Is=Var.sigma_r.*(-State.Ts.^4 + Var.eps_a*Tbl.^4);
    
  %------------------------------
  %Sensible heat flux
   Fs=Var.Ksensible*(State.Ts-Tbl); 
   
  %Flux heat from surface into interior
   Fg=Var.Ki*(State.Ts-State.Ti); 
   State.Ti=State.Ti+Fg*Var.dt./Var.Ci; 

  %---------------------------------
  %Accumulation 
  
  %Snow accumulation where below freezing and over land
   Flp=zeros(1,Var.num);
   Acc=zeros(1,Var.num);
   pl=find(Tbl<Var.K);
   Acc(pl)=Var.P*Var.dt;   %meters
   State.hh=State.hh+Acc;  %meters
   Flp(pl)=(Var.Cp*(Var.K-Tbl(pl)))*Var.rho_i*Var.P; %W/m^2

  %glacial transport
   State.hh=State.hh+State.Flow*Var.dt/(365*24*3600);

  %--------------------------------
  %Heating and ablation

   %Melt all ice at ocean grid boxes
   State.hh(end)=0; 
  
   %Atmospheric heating
   State.Ta=State.Ta+Var.dt*(Sa+Ia+Fs+Flp)/Var.Ca; 

   %Energy available for heating and melting
   B=(Ss+Is-Fs-Fg)*Var.dt;

   %Heating of surface 
   pl=find(State.hh<=0); 
   State.Ts(pl)=State.Ts(pl)+B(pl)./Var.Cs; 

   %Heating and cooling of ice 
   pl=find(State.hh>0 & (State.Ts<Var.K | B<0)); 
   State.Ts(pl)=State.Ts(pl)+B(pl)./Var.Cs;  

   %Ablation: only when there is ice, temperature is at freezing, and energy is added.
   pl=find(State.hh>0 & State.Ts>=Var.K & B>0); 
   Abl=zeros(1,Var.num); 
   Abl(pl)=B(pl)/(Var.Lm*Var.rho_i); 

   %Extra heating converted to melting
   pl=find(State.hh>0 & State.Ts>Var.K); 
   Abl(pl)=Abl(pl)+(State.Ts(pl)-Var.K).*Var.Cs/(Var.Lm*Var.rho_i); 
   State.Ts(pl)=Var.K; 
   
   %Extra melting converted to heating
   pl=find(Abl>State.hh);  
   State.Ts(pl)=State.Ts(pl)+(Abl(pl)-State.hh(pl))*Var.Lm*Var.rho_i./Var.Cs; 
   Abl(pl)=State.hh(pl); 
   State.hh=State.hh-Abl;
   
   if nargout==2,
     %Daily running totals
     Daily.Ti=Daily.Ti+State.Ti;
     Daily.Ts=Daily.Ts+State.Ts;
     Daily.Ta=Daily.Ta+State.Ta;
     Daily.Da=Daily.Da+Da;
     Daily.Is=Daily.Is+Is;
     Daily.Ia=Daily.Ia+Ia;
     Daily.alpha=Daily.alpha+alpha;
     Daily.Acc=Daily.Acc+Acc;
     Daily.Abl=Daily.Abl+Abl;
     Daily.count=Daily.count+1; 
     if rem(t,24*3600)==0,
       %Store daily totals or averages
       day=t/(24*3600)+1;
       StateAll.Ti(day,:)=Daily.Ti/Daily.count; 
       StateAll.Ts(day,:)=Daily.Ts/Daily.count; 
       StateAll.Ta(day,:)=Daily.Ta/Daily.count; 
       StateAll.Da(day,:)=Daily.Da/Daily.count; 
       StateAll.Is(day,:)=Daily.Is/Daily.count; 
       StateAll.Ia(day,:)=Daily.Ia/Daily.count; 
       StateAll.alpha(day,:)=Daily.alpha/Daily.count; 
       StateAll.Acc(day,:)=Daily.Acc; 
       StateAll.Abl(day,:)=Daily.Abl; 
       StateAll.S(day,:)=S; 
     
       %Re-initialize daily totals
       Daily.Ti=zeros(1,Var.num); 
       Daily.Ts=zeros(1,Var.num); 
       Daily.Ta=zeros(1,Var.num); 
       Daily.Da=zeros(1,Var.num); 
       Daily.Is=zeros(1,Var.num); 
       Daily.Ia=zeros(1,Var.num); 
       Daily.alpha=zeros(1,Var.num); 
       Daily.Acc=zeros(1,Var.num); 
       Daily.Abl=zeros(1,Var.num); 
       Daily.count=0;
     end; 
  end;

end;