% ebm_0_dim.m
% 
% Energy Balance Model based on Equation 3.8
% from 'A Climate Modeling Primer', 1st Edition, 1987, 
% A. Henderson-Sellers and K. McGuffie
%
% Written by Tom Huber, 14-JUL-97  for Envision-It Workshop 1997
%
% Requires ebm_o_dim_gui.m to set up values (the present routine does the work)

global SXInit SXLow FirstYear LastYear
global AlbedoVal SVal
global TStart
global AVal BVal
global EbmGuiFig
global AdvancedControls

   A = AVal;                 % Radiative Heat Loss Coefficient
   B = BVal;                 % Radiative Heat Loss Coefficient
   Cm = 2.08E08;              % Heat Capacity
   S = SVal;                 % Solar Constant
   Sigma = 5.6696E-8;        % Stefan-Boltzmann Constant (in W/(m^2K^4))
   nsteps = 500;
   dYears = .05;
   TMin = -15;
   TMax = 15;
   alb_ice = .60;
   alb_land = .32;
   TIce = -10;
   TLand = 10;

   T0 = 273.15;

   BlackBodyButton = get(findobj('Tag','BlackBodyButton'),'Value');
   AlbedoButton = get(findobj('Tag','AlbedoButton'),'Value');

   dt = dYears*365.25*24*60*60;
   CallBack = 0;
   Time = zeros(nsteps,1);

   Temp(1) = TStart;
   Time(1) = 0;
   for n = 2:nsteps;
      T = Temp(n-1);
      Time(n) = (n-1)*dYears;
      if Time(n)>=FirstYear & Time(n)<LastYear
        SX = SXLow;
      else
        SX = SXInit;
      end
      if AlbedoButton == 0
        albedo = AlbedoVal;
      else
        if T<TIce
          albedo = alb_ice;
        elseif T>TLand
          albedo = alb_land;
        else
          albedo = alb_ice + (alb_land-alb_ice)*(T-TIce)/(TLand-TIce);
        end
      end

     PGain = S*SX*(1-albedo)/4;
     if BlackBodyButton == 0
        PLoss = A + B*T;
     else
        PLoss = Sigma*(T+T0)^4;
     end
     DeltaTemp = (PGain - PLoss)*(dt/Cm);
     Temp(n) = T + DeltaTemp;
     AlbedoTime(n) = albedo;
   end

   hold off
   plot (Time,Temp);
   xlabel('Time (Years)');
   ylabel('Temperature (Degrees C)');

   disp(['Final Temperature = ' num2str(Temp(nsteps)) ' Degrees Centigrade']);
   
%  matl5
   HH = findobj(gcf,'Tag','TFinal');
   set(HH,'String',[num2str(floor(Temp(nsteps)*1000)/1000) ' C'])
   if BlackBodyButton == 0
      TSteady=(S*SX*(1-albedo)/4-A)/B;
   else
      TSteady=(S*SX*(1-albedo)/4/Sigma)^.25-T0;
   end
   disp(['Steady-State Temperature = ' num2str(TSteady) ' Degrees Centigrade']);
   
%  matl5
   HH = findobj(gcf,'Tag','TSteady');
   set(HH,'String',[num2str(floor(TSteady*1000)/1000) ' C']);