function [dy,lambda,data]=model_fun(t,y)
dy = zeros(2,1);    % a column vector
% time step = one year
% constant definition
ro=1000;         % water density
dd=900;          % depth of deep water
dm=100;          % depth of mixing layer
lambda=0.67;     % climate sensitivity 
k=1e-4;          % water Diffusivity 
Cp=4218;         % specific heat copacity for water
Cd=ro*Cp*dd;     % heat capacity for deep leyer
Cm=ro*Cp*dm;     % heat capacity for mixing layer
data=forcingdata;
year=data.year;
%data.TOT=0*data.TOT;
%data.TOT(5:end)=-2;
%data.TOT(5:end)=(year(5:end)-year(5))/100;
FOR=data.TOT;
dF=interp1(year,FOR,t/(365*24*3600));
dy(2)=k*ro*Cp/Cd*(y(1)-y(2))/(2*(-dm+dd));
dy(1)=dF/Cm-y(1)/(Cm*lambda)-Cd/Cm*dy(2);
