function doall=THC_winton_1993
% The 3 box THC model of Winton 1993, used to demonstrate relaxation
% oscillations of the THC for the climate dynamics course EPS231.  Use
% matlab ODE solver to integrate equations.  There is an option of
% adding noise, search for noise_amplitude below.  Eli Tziperman
%
%  Reference:
%  author =	 {M. Winton},
%  title =	 {Deep decoupling oscillations of the oceanic
%                  thermohaline circulation},
%  booktitle =	 {Ice in the climate system},
%  publisher =	 "Springer Verlag",
%  volume =	 12,
%  editor =	 {W. R. Peltier},
%  series =	 {NATO ASI Series I: Global Environmental Change},
%  pages =	 {417-432},
%  year =	 1993

clear all

% profile on

%   --------------------------------------
%   |  	1	      M_lh |	    2  	 |
%   |  L    ^  	       <------>	  ^   H	 |
%   |	    |    	   |   	  |    	 |
%   |-------|--------------|------|------|
%   |  	    | M_ld	       	  | M_hd |
%   |  3    v		       	  v    	 |
%   |			       	       	 |
%   |		D	       	       	 |
%   |		               	       	 |
%   |			    	       	 |
%   |			    	       	 |
%   --------------------------------------

% set model parameters and initial conditions:
p=set_parameters;
X_init=[p.T_init,p.S_init];

%% integrate the box model equations:
disp('running the box model...')
% stiff solvers: ode15s, ode23s, ode23t, and ode23tb,
[time,X]=ode15s(@(t,X) box_model_rhs(t,X,p),[0,p.integration_time],X_init);

time=time/p.year;
T=X(:,1:3);
S=X(:,4:6);
rho=S*NaN;
N=size(T(:,1));
pressure=1.0;
for n=1:N;for i=1:3;rho(n,i)=density(T(n,i),S(n,i),pressure)-1000;end;end

%% save results for additional analysis:
save('Output/timeseries.mat','time','T','S','rho');

%% plot results:
figure(1)
clf

axes('position', [0.15,0.72,0.8,0.25]); % [left, bottom, width, height])
plot(time,T(:,1),time,T(:,2),time,T(:,3))
hold on
set_axis_limits(gca,time,[T(:,1),T(:,2),T(:,3)])
h(1)=legend('T_1','T_2','T_3');
h(2)=gca;
h(3)=ylabel('T');
set(gca,'xticklabel',' ')

axes('position', [0.15,0.42,0.8,0.25]); % [left, bottom, width, height])
plot(time/1e3,S(:,1),time/1e3,S(:,2),time/1e3,S(:,3))
hold on
set_axis_limits(gca,time/1e3,[S(:,1),S(:,2),S(:,3)])
h(4)=legend('S_1','S_2','S_3');
h(5)=ylabel('S'); 
h(6)=gca;
set(gca,'xticklabel',' ')

axes('position', [0.15,0.12,0.8,0.25]); % [left, bottom, width, height])
plot(time/1e3,rho(:,1),time/1e3,rho(:,2),time/1e3,rho(:,3));%,time,rho(2,:)-rho(3,:))
h(7)=legend('\sigma_1','\sigma_2','\sigma_3');%,'\rho_2-\rho_3');
hold on
set_axis_limits(gca,time/1e3,[rho(:,1),rho(:,2),rho(:,3)])
h(8)=ylabel('density'); 
h(9)=xlabel('time (kyr)');
h(10)=gca;
set(h,'FontName','Helvetica','FontWeight','normal','FontSize',12);

%% save as pdf:
set(gcf, 'PaperUnits', 'inches'); set(gcf, 'PaperSize', [6 4]);
set(gcf, 'PaperPosition', [0 0 6 4]); % [left, bottom, width, height];
saveas(gcf,sprintf('Output/winton.pdf'));
saveas(gcf,sprintf('Output/winton.fig'));

return
figure
subplot(2,2,1)
N=length(T(:,3));
for n=floor(N/3):N
    if rho(n,2)<rho(n,3);
        plot(T(n,3),S(n,3),'r.');
        hold on
    else
        plot(T(n,3),S(n,3),'g.');
    end
end
h1=xlabel('T');
h2=ylabel('S');
h3=title('red: non convecting; green: convecting');
set([h1,h2,h3,gca],'fontsize',16)

print -dpdf Output/Winton_1993_figure2.pdf


% profile report

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rhs=box_model_rhs(t,X,p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

T=X(1:3);
S=X(4:6);

rhs=[time_rate_of_change_T(t,T,S,p),time_rate_of_change_S(t,T,S,p)]';
  


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rho=linear_density(p,T,S)

rho=-p.alpha*T+p.beta*S;


function rsw=density(t,s,pp)
%-------------------------------------------------------------------
% PH_DENSITY $Revision: 1.0 $ $Date: 1998/12/07 $
%          Copyright (C) GEOMAR, Roger Luff 1998.
%
% USAGE: rsw = ph_density(t,s,pp)
%
% DESCRIPTION:
%   density calculates the density of seawater using the unesco
%   (1983) formulae, which are valid for a temperature range 
%   of 0-40 deg c and 0.5-43 ppt
%
% INPUT:
%  t     : temperature of seawater at the specific location [deg C]
%  s     : salinity of seawater at the specific location  [PSU]
%  pp    : pressure at the specific location  [atm, (with 1 atm
%          at sea surface)]
%
% OUTPUT:
%  rsw   : density of a "seawater" [kg/m^3]
%
% ORIGINAL AUTHOR: 
%   Bernard Boudreau 
%   Dept. Oceanography 
%   Dalhousie University
%   Halifax, NS B3H 4J1,Canada 
%   bernie@boudreau.ocean.dal.ca
% MATLAB VERSION
%   Roger Luff 98-12-07
%   GEOMAR Research Centre for Marine Geosciences
%   Wischhofstr. 1-3, D-24148 Kiel, Germany
%   visumod@technologist.com
%
% DISCLAIMER:
%   This software is provided "as it is" without warranty of any kind.
%
% REFERENCES:
%  Boudreau B. P. (1996) "A method-of-lines code for carbon and 
%   nutrient diagenesis in aquatic sediments". Computers & Geosciences
%   22(5), 479-496.
%  Luff, Haeckel and Wallmann (submitted) "Robust and fast FORTRAN
%   and MATLAB libraries to calculate pH distributions in 
%   marine systems", Computers & Geosciences.
%  Fofonoff P. and Millard R. C. (1983) "Algorithms for computation 
%   of fundamental properties of seawater." UNESCO Technical Paper 
%   in Marine Science 44, 53.
%-------------------------------------------------------------------

%-------------------------------------------------------------------
% CHECK INPUT ARGUMENTS
%-------------------------------------------------------------------
if nargin ~=3
   error('ph_density.m: Must pass 3 parameters')
end
%-------------------------------------------------------------------
%  convert pressure to bars of sea pressure only
%-------------------------------------------------------------------
p = (pp-1.0)*1.013;
%-------------------------------------------------------------------
%  density of pure reference water
%-------------------------------------------------------------------
a0 = 999.842594;
a1 = 6.793952d-02;
a2 = - 9.095290d-03;
a3 = 1.001685d-04;
a4 = - 1.120083d-06;
a5 = 6.536332d-09;
rw = ((((a5*t + a4)*t + a3)*t + a2)*t + a1)*t + a0;
%-------------------------------------------------------------------
%  density at t, s and 1 atm
%-------------------------------------------------------------------
b0 = 8.24493d-01;
b1 = - 4.0899d-03;
b2 = 7.6438d-05;
b3 = - 8.2467d-07;
b4 = 5.3875d-09;
c0 = - 5.72466d-03;
c1 = 1.0227d-04;
c2 = - 1.6546d-06;
d0 = 4.8314d-04;
rts = rw + ((((b4*t + b3)*t + b2)*t +b1)*t + b0)*s + ((c2*t + c1)*t + c0)*sqrt(s*s*s) + d0*s*s;
if p==0.0
   rsw = rts;
else
%-------------------------------------------------------------------
% density at t, s and p
%-------------------------------------------------------------------
  e0 = 19652.21;
  e1 = 148.4206;
  e2 = - 2.327105;
  e3 = 1.360477d-02;
  e4 = - 5.155288d-05;
  kw = (((e4*t + e3)*t + e2)*t + e1)*t + e0;
  h0 = 3.239908;
  h1 = 1.43713d-03;
  h2 = 1.16092d-04;
  h3 = - 5.77905d-07;
  aw = ((h3*t + h2)*t + h1)*t + h0;
  k0 = 8.50935d-05;
  k1 = - 6.12293d-06;
  k2 = 5.2728d-08;
  bw = (k2*t + k1)*t + k0;
  m0 = - 9.9348d-07;
  m1 = 2.0816d-08;
  m2 = 9.1697d-010;
  b = bw + ((m2*t + m1)*t + m0)*s;
  i0 = 2.2838d-03;
  i1 = - 1.0981d-05;
  i2 = - 1.6078d-06;
  j0 = 1.91075d-04;
  a = aw + ((i2*t + i1)*t + i0)*s + j0*sqrt(s*s*s);
  f0 = 54.6746;
  f1 = - 0.603459;
  f2 = 1.09987d-02;
  f3 = - 6.1670d-05;
  g0 = 7.944d-02;
  g1 = 1.6483d-02;
  g2 = - 5.3009d-04;
  kst0 = kw + (((f3*t + f2)*t + f1)*t + f0)*s + ((g2*t + g1)*t + g0)*sqrt(s*s*s);
  kstp = kst0 + a*p + b*p*p;
  rsw =  rts/(1.0 - p/kstp);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dT_dt=time_rate_of_change_T(t,T,S,p)
% rate of change of temperature:

% Calculate densities:
% rho=linear_density(p,Tnow,Snow);
pressure=1.0;
for i=1:3; rho(i)=density(T(i),S(i),pressure)/1000.0-1.0; end

% check first if convection occurs:
if rho(3)<rho(2) C_h=100; else C_h=1; end
if rho(3)<rho(1) C_l=100; else C_l=1; end
T_noise_forcing=p.noise_forcing(1,max(1,ceil(t/p.dt_noise)));
dT_dt=[  -p.T_noise_factor*T_noise_forcing/p.V(1)+p.gamma_T*(p.Tstar(1)-T(1))+p.M_lh*(T(2)-T(1))/p.V(1) + C_l*p.M_ld*(T(3)-T(1))/p.V(1) ...
       , -p.T_noise_factor*T_noise_forcing/p.V(2)+p.gamma_T*(p.Tstar(2)-T(2))+p.M_lh*(T(1)-T(2))/p.V(2) + C_h*p.M_hd*(T(3)-T(2))/p.V(2) ...
       , C_h*p.M_hd*(T(2)-T(3))/p.V(3) + C_l*p.M_ld*(T(1)-T(3))/p.V(3)];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dS_dt=time_rate_of_change_S(t,T,S,p)
% rate of change of salinity:

% Calculate densities:
% rho=linear_density(p,Tnow,Snow);
pressure=1.0;
for i=1:3; rho(i)=density(T(i),S(i),pressure)/1000.0-1.0; end

% check first if convection occurs:
if rho(3)<rho(2) C_h=100; else C_h=1; end
if rho(3)<rho(1) C_l=100; else C_l=1; end
% note use of zero net salt forcing:
S_noise_forcing=p.noise_forcing(2,max(1,ceil(t/p.dt_noise)));
dS_dt=[ p.S_noise_factor*S_noise_forcing/p.V(1)+p.FW(1)/p.V(1)+p.gamma_S*(p.Sstar(1)-S(1))+p.M_lh*(S(2)-S(1))/p.V(1) + C_l*p.M_ld*(S(3)-S(1))/p.V(1) ...
	, -p.S_noise_factor*S_noise_forcing/p.V(2)+p.FW(2)/p.V(2)+p.gamma_S*(p.Sstar(2)-S(2))+p.M_lh*(S(1)-S(2))/p.V(2) + C_h*p.M_hd*(S(3)-S(2))/p.V(2) ...
	, C_h*p.M_hd*(S(2)-S(3))/p.V(3) + C_l*p.M_ld*(S(1)-S(3))/p.V(3)];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function set_axis_limits(h,x,y)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set axis limits tighter around min/max of y data than the
%% typical Matlab default.
%% If plotting more than one array using more than one plot call, use:
%% set_axis_limits(gca,x,[y1,y2])
%% If y1,y2 are column vectors, use
%% set_axis_limits(gca,x,[y1',y2'])
%% If wish to add min/max values in addition to the data (e.g. if
%% data are a constant, and a wider range is desired):
%% set_axis_limits(gca,x,[ymin,ymax,y1,y2])

%% set y limits:
y=reshape(y,1,[]);
V=axis; 
a1=min(y); 
a2=max(y);
V(3)=a1-(a2-a1)*0.15; 
V(4)=a2+(a2-a1)*0.15;
if V(3)==V(4); 
  V(3)=a1*0.9;
  V(4)=a2*1.1;
  if V(3)==0.0
    V(3)=-1;
    V(4)=1;
  end
end
%% set x limits:
V(1)=x(1);
V(2)=x(length(x));

axis(h,V);
set(h,'YTickMode','auto');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p=set_parameters

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                     cgs!!                       %%
%%                     =====                       %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% create output directory if needed:
if exist('Output')~=7
  disp('creating Output/ directory');
  mkdir('Output');
end

% box geometry:
% north-south size:
p.L(1)=5000.0e5;
p.L(2)=2000.0e5;
p.L(3)=p.L(1)+p.L(2);
% depth:
p.H(1)=250.0e2;
p.H(2)=p.H(1);
p.H(3)=1000.0e2;
% width:
p.W=4000.0e5;
% volumes:
p.V=p.L(:).*p.H(:)*p.W;
% areas:
p.Area(1)=p.L(1)*p.W;
p.Area(2)=p.L(2)*p.W;

% some parameters:
p.Sv=1.e12;
p.year=3600.0*24.0*365.0;
p.day=86400;
p.week=7*p.day;
% linearized equation of state:
p.alpha=1668.0e-7;
p.beta=0.781e-3;

% diffusion coefficients:
p.M_lh=p.H(1)*p.W*0.5*(p.L(1)+p.L(2))/(5*p.year);
p.M_ld=p.Area(1)*0.5*(p.H(1)+p.H(2))/(200*p.year);
p.M_hd=p.Area(2)*0.5*(p.H(1)+p.H(2))/(200*p.year);

% SST restoring time:
p.gamma_T=1.0/(30.0*86400.0);
p.gamma_S=0.0/((50*86400.0)*(p.H(1)/50.e2));
p.FW(1)=(180.0/p.year)*35.0*p.Area(1);
p.FW(2)=-p.FW(1);
% the temperature to which the model surface temperature is restored:
p.Tstar=[15.0,0.0]; 
% the salinity to which the model surface salinity is restored (only
% if p.gamma_S is not zero):
p.Sstar=[36.0,34.0];

% parameters for numerically solving the steady state problem:
% integration time:
p.integration_time=15000.0*p.year;

% initial conditions:
p.T_init=[20.0 1 1];
p.S_init=[35.5 34.5 34.5];

% amplitude of noise added to surface boxes:
p.add_noise=0;  % 0 or 1
if p.add_noise
  p.T_noise_factor=0;
  p.S_noise_factor=2;
else
  p.T_noise_factor=2e15; % 2e15;
  p.S_noise_factor=2e13; %2e13;
end
% uncomment the following lines to test the density calculation:
% pressure=1.0;
% for i=1:3
%   rho_linear(i)=linear_density(p,T(i),S(i));
%   rho(i)=density(T(i),S(i),pressure)/1000.0-1;
% end
% disp(' test T,S,density:')
% T,S,rho
% disp(' test T,S,linear density:')
% T,S,rho_linear

%% create random forcing:
p.dt_noise=10*p.day;
N=ceil(p.integration_time/p.dt_noise);
noise_forcing=zeros(2,N);
if p.add_noise
  %% create noise at a daily resolution and specify correlation time.
  %% one for temperature and one for salinity:
  time_noise=[1:N]*p.dt_noise/p.year;
  p.noise_correlation_time=p.year*100;
  %% if any of the above 4 parameters changed, set the following to 1
  %% to generate new noise, otherwise set to 0 to read saved noise:
  generate_noise=0;
  if generate_noise 
    disp(['generating noise forcing, consider setting generate_noise=0 ' ...
          'next run to save time...']);
    %% For efficiency, set up the random numbers to be used for
    %% surface_noise (fresh water and heat flux noise) function in advance,
    %% rather than calling random() at every time step:
    %% First, build a matrix of random values:
    %%tmp_noise_forcing=rand(2,N);
    tmp_noise_forcing=random('Normal',0.0,1.0,2,N);
    %% next, create a first order markov process with the specified
    %% correlation time from the above random numbers:
    noise_forcing(:,1)=tmp_noise_forcing(:,1);
    R_noise=exp(-p.dt_noise/p.noise_correlation_time);
    tmp_R2=sqrt(1-R_noise^2);
    for n=1:N-1
      noise_forcing(:,n+1) = ...
          R_noise*noise_forcing(:,n) + tmp_R2*tmp_noise_forcing(:,n);
    end
    clear tmp_R2 tmp_noise_forcing
  
    if 1
      windowSize=3650;
      for i=1:2
        noise_forcing(i,:)=filtfilt(ones(1,windowSize)/windowSize,1,noise_forcing(i,:));
      end
    end
    N_noise=N;
    save('Output/noise.mat','noise_forcing','N_noise')
  else
    disp('loading saved noise forcing...');
    load('Output/noise.mat')
    fprintf(1,'N=%d (calculated for present run), N_noise (loaded)=%d\n',N,N_noise);
    if N>N_noise
      disp('*** error: saved noise not at right length, abort (^C) and');
      disp('    re-generate noise by setting generate_noise to 1.');
      pause
    end
  end

  %% figure
  %% plot(time_noise,p.noise_forcing(1,:));
  %% pause
end
p.noise_forcing=noise_forcing;
clear noise_forcing
