% Run EBM (One-dimensional atmospheric energy balance model)
% This program is a one-dimensional energy balance model that calculates
% the surface air temperature as a function of latitude. It is based on the
% lecture notes by Stocker (2008) [problem 9]

function [tatmnew,tbndnew] = runmodel_ebm(tatmold,tbndold,EBM_time,Tocenew,thau,atmeffyr)

%% Load constants
persistent ebm_con1 ebm_con7

if EBM_time==1 % Load constant during first timestep
    load('constants.mat');
    ebm_con1=[jmt;stebol;dxu;ceff;dxt;ebm_res;yr2sec;dt_ebm];
    ebm_con7=[solin;apln;dyu;csu;diff_cnt;cst;dyt];
end

%% Initialize arrays
netswpln=zeros(1,ebm_con1(1));
lwocs=zeros(1,ebm_con1(1));
lwatm=zeros(1,ebm_con1(1));
grad_tatm = zeros(1,ebm_con1(1));
heddy=zeros(1,ebm_con1(1));
div_eddy=zeros(1,ebm_con1(1));
source=zeros(1,ebm_con1(1));
yp=zeros(1,ebm_con1(1));
oc2atm=zeros(1,ebm_con1(1));
diff_oa=zeros(1,ebm_con1(1));
tatmnew=zeros(1,ebm_con1(1));
tbndnew=zeros(1,ebm_con1(1));

%% Run EBM
if (EBM_time>1 || eq_length==0)
    % Shortwave energy budget
    netswpln(2:ebm_con1(1)-1) = ebm_con7(1,2:ebm_con1(1)-1).*(1.0 - ebm_con7(2,2:ebm_con1(1)-1));
    
    % Longwave energy budget
    % Compute outgoing longwave radiation of atmosphere
    lwatm(2:ebm_con1(1)-1) = ebm_con1(2).*tatmold(2:ebm_con1(1)-1).^4;
    % Compute outgoing longwave radiation at ocean surface
    lwocs(2:ebm_con1(1)-1) = netswpln(2:ebm_con1(1)-1)+lwatm(2:ebm_con1(1)-1);

    % Set source term for energy balance equation
    source(2:ebm_con1(1)-1)   = (2-2*atmeffyr(2:ebm_con1(1)-1)).*lwocs(2:ebm_con1(1)-1) - 2*lwatm(2:ebm_con1(1)-1);

    % Compute meridional heat transport
    % Compute meridional temperature gradient (temperature difference / distance)
    grad_tatm(3:ebm_con1(1)-1)   = (tatmold(3:ebm_con1(1)-1) - tatmold(2:ebm_con1(1)-2))./ebm_con7(3,3:ebm_con1(1)-1);
    % Compute meridional heat transport by transient eddies (cross-sectional area * heat capacity * diffusivity * temperature gradient)
    heddy(3:ebm_con1(1)-1)    = ebm_con1(3).*ebm_con7(4,3:ebm_con1(1)-1).*ebm_con1(4).*-ebm_con7(5,2:ebm_con1(1)-2).*grad_tatm(3:ebm_con1(1)-1);
    % Compute divergence of meridional heat transport after Sellers (1969) and North (1975)
    dareat = ebm_con1(5)*ebm_con7(6,2:ebm_con1(1)-1).*ebm_con7(7,2:ebm_con1(1)-1);
    div_eddy(2:ebm_con1(1)-1) = (heddy(3:ebm_con1(1)) - heddy(2:ebm_con1(1)-1))./dareat;

    % Ocean to Atmosphere heat diffusion
    diff_oa(3:ebm_con1(1)-2)   = Tocenew(1:ebm_con1(1)-4)-tbndold(3:ebm_con1(1)-2);
    oc2atm(3:ebm_con1(1)-2)    = diff_oa(3:ebm_con1(1)-2).*thau/(ebm_con1(6)*ebm_con1(7));
    
    % Total energy budget including forced temperature increase    
    % Compute "slope"
    yp(2:ebm_con1(1)-1) = (source(2:ebm_con1(1)-1)-div_eddy(2:ebm_con1(1)-1)+oc2atm(2:ebm_con1(1)-1))/ebm_con1(4);

    % Step atmospheric temperature forward in time (by "Euler forward" or "forward-in-time" method)
    tatmnew(2:ebm_con1(1)-1) = tatmold(2:ebm_con1(1)-1) + ebm_con1(8)*yp(2:ebm_con1(1)-1);
    tbndnew(2:ebm_con1(1)-1) = (lwocs(2:ebm_con1(1)-1)./max(0,ebm_con1(2))).^(1/4);
else % For initial timestep:
    netswpln(2:ebm_con1(1)-1) = ebm_con7(1,2:ebm_con1(1)-1).*(1.0 - ebm_con7(2,2:ebm_con1(1)-1));
    % Compute outgoing longwave radiation at ocean surface
    lwocs(2:ebm_con1(1)-1) = netswpln(2:ebm_con1(1)-1)./atmeffyr(2:ebm_con1(1)-1);
    % Compute outgoing longwave radiation of atmosphere
    lwatm(2:ebm_con1(1)-1) = (1-atmeffyr(2:ebm_con1(1)-1)).*lwocs(2:ebm_con1(1)-1);
    % Compute atmospheric temperatures
    tatmnew(2:ebm_con1(1)-1) = (lwatm(2:ebm_con1(1)-1)./max(0,ebm_con1(2))).^(1/4);
    % Calculate atmospheric surface temperatures
    tbndnew(2:ebm_con1(1)-1) = (lwocs(2:ebm_con1(1)-1)./max(0,ebm_con1(2))).^(1/4);
end
end
