% Reference: "How Model Complexity Influences Sea Ice Stability", 

% T.J.W. Wagner & I. Eisenman, J Clim (2015)

%

% WE15_EBM_fast.m:

% This code describes the EBM as discussed in Sec. 2b of the article above, 

% hereafter WE15. Here we use central difference spatial integration and 

% Implicit Euler time stepping.

%

% The code WE15_EBM_simple.m, on the other hand, uses a simpler formulation

% of the diffusion operator and time stepping with Matlab's ode45.

%

% Parameters are as described in WE15, table 1. Note that we do not include 

% ocean heat flux convergence or a seasonal cylce in the forcing 

% (equivalent to S_1 = F_b = 0 in WE15). This code uses an ice albedo when 

% T<0 (WE15 instead uses the condition E<0, which is appropriate for the 

% inclusion of a seasonal cycle in ice thickness). In this code, we define 

% T = Ts - Tm, where Ts is the surface temperature and Tm the melting point 

% (WE15, by contrast, defines T = Ts).

%

% Till Wagner & Ian Eisenman, Mar 15

% tjwagner@ucsd.edu or eisenman@ucsd.edu

%

%%Model parameters (WE15, Table 1 and Section 2d) -------------------------

D  = 0.6;     % diffusivity for heat transport (W m^-2 K^-1)

A  = 193;     % OLR when T = T_m (W m^-2)

B  = 2.1;     % OLR temperature dependence (W m^-2 K^-1)

cw = 9.8;     % ocean mixed layer heat capacity (W yr m^-2 K^-1)

S0 = 420;     % insolation at equator  (W m^-2)

S2 = 240;     % insolation spatial dependence (W m^-2)

a0 = 0.7;     % ice-free co-albedo at equator

a2 = 0.1;     % ice=free co-albedo spatial dependence

ai = 0.4;     % co-albedo where there is sea ice

F  = 0;       % radiative forcing (W m^-2)

% -------------------------------------------------------------------------

n  = 50;

nt = 0.5;     % can be small since time integration is implicit

dur= 30;

dt = 1/nt;

%%Spatial Grid ------------------------------------------------------------

dx = 1/n;     %grid box width

x = (dx/2:dx:1-dx/2)';  %native grid

%%Diffusion Operator (WE15, Appendix A) -----------------------------------

xb = (dx:dx:1.0-dx)';  

lambda=D/dx^2*(1-xb.^2); L1=[0; -lambda]; L2=[-lambda; 0]; L3=-L1-L2;

diffop = - diag(L3) - diag(L2(1:n-1),1) - diag(L1(2:n),-1);

 

S = S0-S2*x.^2;        % insolation: WE15, eq.3 with S_1 = 0

T = 10*ones(n,1);      % initial condition (constant temp. 10C everywhere)

aw = a0-a2*x.^2;       % albedo for open water

 

allT = zeros(dur*nt,n);

t = linspace(0,dur,dur*nt);

% integration over time using implicit difference and

% over x using central difference (through diffop)

for i = 1:dur*nt

    a = aw.*(T>0)+ai*(T<0); % WE15, eq.4

    C = a.*S-A+F;                    

    % Governing equation [cf. WE15, eq. (2)]:

    %    T(n+1) = T(n) + dt*(dT(n+1)/dt), with c_w*dT/dt=(C-B*T+diffop*T)

    % -> T(n+1) = T(n) + dt/cw*[C-B*T(n+1)+diff_op*T(n+1)]

    % -> T(n+1) = inv[1+dt/cw*(1+B-diff_op)]*(T(n)+dt/cw*C)

    T = (eye(n) + dt/cw*(B*eye(n)-diffop))\(T+dt/cw*C);

    allT(i,:)=T;

end

% plot results ------------------------------------------------------------

figure(1), clf

subplot(1,2,1)

plot(t,allT)

xlabel('$t$ (years)'); ylabel('$T$ (for different latitudes)')

xlim([0 t(end)])

subplot(1,2,2)

plot(x,T, 'x-')

xlabel('$x = sin \theta$'); ylabel('$T$ ($^\circ$C)')