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

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

%

% WE15_EBM_simple.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 

% time stepping with MATLAB's ode45.

%

% The code WE15_EBM_fast.m, on the other hand, uses a faster, but more 

% complicated formulation of the diffusion operator and Implicit Euler time

% stepping.

%

% 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

%

function [x_out,T_out] = WE15_EBM_simple

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

global  D B n cw aw ai S A F x dx 

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;      % grid resolution (number of points between equator and pole)

x  = linspace(0,1,n);

dx = 1/(n-1);

S  = S0-S2*x.^2;      % insolation [eq. (3) in WE15, with S_1 = 0]

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

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

tspan = [0 30];       % timespan in years [t_start t_end]

% time integration---------------------------------------------------------

[t,T] = ode45(@rhsWE15,tspan,T0);

if nargout==0 % plot results-----------------------------------------------

    figure(1), clf

    subplot(1,2,1)

    plot(t,T)

    xlabel('$t$ (years)'); 

    ylabel('$T$ ($^\circ$C, for different latitudes)')

    xlim([0 t(end)])

    box on

    subplot(1,2,2)

    hold all

    plot(x,T(end,:), 's-')

    xlabel('$x = sin \theta$'); 

    ylabel('$T$ ($^\circ$C)')

    box on

else % save results--------------------------------------------------------

    x_out=x;

    T_out=T;

end

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

% ODE with spatial finite differencing for ode45 --------------------------

function dTdt = rhsWE15(~,T)

global  D B n cw aw ai S A F x dx

Tdot = zeros(n,1);

% forcing

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

C = alpha.*S-A+F;                

% solve c_wdT/dt = D(1-x^2)d^2T/dx^2 - 2xDdT/dx + C - BT [cf. WE15, eq.(2)]

% use central difference

Tdot(2:n-1) = (D/dx^2)*(1-x(2:n-1).^2)'.*(T(3:n)-2*T(2:n-1)+T(1:n-2)) ...

              -(D*x(2:n-1)/dx)'.*(T(3:n)-T(1:n-2));

% boundary condition: let dT/dt(x=0) = 0

Tdot(1) = D*(2*(T(2)-T(1))/dx^2);   % WE15, eq.(7)

% use asymmetric (backward) spatial difference for the pole

Tdot(n) = - D*(2*x(n)*(T(n)-T(n-1))/dx);    %O(dt) accuracy at x=1

% temperature tendency

dTdt = (Tdot+C' - B*T)/cw;