function doall=THC_3box
% a 3 box THC model

% comments:
% April 10, 2005: fixed a sign mistake in the ds/dt equation for
% u>0 (reported by Laure).
clear all; close all

% solution of the full model:
global Tplot1 Splot1 uplot1

% profile on

%   --------------------------------------
%   |  			   |	      	 |
%   |  1    ^  	       	======>	      2	 |
%   |	    |		   |   	       	 |
%   |-------|--------------|   	       	 |
%   |  	    |		   |   	       	 |
%   |  3     		   |   	       	 |
%   |			   |   	       	 |
%   |			   |   	       	 |
%   |		       <=======	       	 |
%   |			   |	       	 |
%   |			   |	       	 |
%   --------------------------------------


% set model parameters and initial conditions:
[L,H,delta,W,V,Area,Sv,year,u0,alpha,beta,gamma_T,gamma_S ...
      ,FW,Tstar,Sstar,delta_t,N_max,ubar,T,S,F ...
      ,random_matrix,noise_correlation_time] ...
    = set_parameters;

% run the model to a steady state:
disp(' running the box model for 10000 years with no noise to steady state:')
% can add stochastic forcing to the box model surface forcing:
sqrt_C=zeros(6); which_noise='none';
% set integration time:
N_max=floor((10000*year/delta_t)/100)*100;
% the temperature to which the model surface temperature is restored:
Tstar=[22.0,0.0]; 
% run the box model:
plot_anomalies=0; Tbar=0; Sbar=0; ubar=0;
[var_T,var_S,var_THC,Tbar,Sbar,ubar,Tplot1,Splot1,uplot1]= ...
    box_model(T,S,u0,alpha,beta ...
    ,gamma_T,gamma_S,delta,year,FW,delta_t,N_max,Tstar,Sstar ...
    ,V,H,Area,sqrt_C,random_matrix,which_noise,Tbar,Sbar,ubar);
%% save fig as pdf:
set(gcf, 'PaperUnits', 'inches'); set(gcf, 'PaperSize', [5 4]);
set(gcf, 'PaperPosition', [0 0 5 4]); % [left, bottom, width, height];
saveas(gcf,sprintf('THC_3box_figure1.pdf'));

disp(' Steady state solution is:')
Tbar
Sbar
ubar

% now run the model again for a short period, starting from the steady
% state plus a small perturbation:
disp([' running the box model for a short time from the steady state' ...
      ' plus a small perturbation:'])
% set initial conditions:
temperature_perturbation=[0.2,0.2,0.2];
salinity_perturbation=[0.2,0.2,0.2];
T=Tbar+temperature_perturbation; S=Sbar+salinity_perturbation;
% can add stochastic forcing to the box model surface forcing:
sqrt_C=zeros(6); which_noise='none';
% set integration time:
N_max=floor((600*year/delta_t)/100)*100;
% run the box model:
plot_anomalies=1;
[var_T,var_S,var_THC,Tbar,Sbar,ubar,Tplot1,Splot1,uplot1]= ...
    box_model(T,S,u0,alpha,beta ...
    ,gamma_T,gamma_S,delta,year,FW,delta_t,N_max,Tstar,Sstar ...
    ,V,H,Area,sqrt_C,random_matrix,which_noise,plot_anomalies,Tbar,Sbar,ubar);
%% save fig as pdf:
set(gcf, 'PaperUnits', 'inches'); set(gcf, 'PaperSize', [5 4]);
set(gcf, 'PaperPosition', [0 0 5 4]); % [left, bottom, width, height];
saveas(gcf,sprintf('THC_3box_figure2.pdf'));


% profile report

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [var_T,var_S,var_THC,Tbar,Sbar,ubar,Tplot,Splot,uplot]=box_model ...
    (T,S,u0,alpha,beta ...
    ,gamma_T,gamma_S,delta,year,FW,delta_t,N_max,Tstar,Sstar ...
    ,V,H,Area,sqrt_C,random_matrix,which_noise,plot_anomalies,Tbar,Sbar,ubar)

% initialize:
Tnew=T;
Snew=S;
Told=T;
Sold=S;
Tnow=T;
Snow=S;

% initialize arrays containing plot data:
Tplot=zeros(3,N_max);
Splot=zeros(3,N_max);
uplot=zeros(1,N_max);

for n=1:N_max

  % Calculate densities:
  rho=linear_density(alpha,beta,Tnow,Snow);
  % uncomment the following to use nonlinear eqn of state:
  % pressure=1.0;
  % for i=1:3
  % rho(i)=density(Tnow(i),Snow(i),pressure)/1000.0-1.0;
  % end
  
  % Calculate transport:
  u=transport(u0,rho,delta);
  
  % leap frog time stepping:

  Tnew=Told+2*delta_t*time_rate_of_change_T ...
      (n,Tnow,u,Tstar,gamma_T,V,H,sqrt_C,random_matrix,which_noise);
  Snew=Sold+2*delta_t*time_rate_of_change_S ...
      (n,Snow,u,Sstar,gamma_S,FW,sqrt_C,V,H,random_matrix,which_noise);
  
  % Robet filter:
  robert_filter_coeff=0.25;
  Tnow = Tnow + robert_filter_coeff * ( Tnew - 2*Tnow + Told );
  Snow = Snow + robert_filter_coeff * ( Snew - 2*Snow + Sold );

  Told=Tnow;
  Tnow=Tnew;
  Sold=Snow;
  Snow=Snew;
  
  % plotting data:
  Tplot(:,n)=Tnow';
  Splot(:,n)=Snow';
  uplot(n)=u/1.e12;
end

if plot_anomalies
  for ii=1:3
    Tplot(ii,:)=Tplot(ii,:)-Tbar(ii); 
    Splot(ii,:)=Splot(ii,:)-Sbar(ii); 
  end
  uplot=uplot-ubar/1.e12;
end

time=(1:n)*delta_t/year;

figure
axes('position', [0.1,0.7,0.8,0.2]); % [left, bottom, width, height])
title(' results of full box model')
hold on
plot(time,Tplot(1,:),time,Tplot(2,:),time,Tplot(3,:))
hold on
legend('T_1','T_2','T_3')
xlabel('years'); 
ylabel('T'); 
axes('position', [0.1,0.4,0.8,0.2]); % [left, bottom, width, height])
plot(time,Splot(1,:),time,Splot(2,:),time,Splot(3,:))
hold on
legend('S_1','S_2','S_3')
xlabel('years'); 
ylabel('S'); 
axes('position', [0.1,0.1,0.8,0.2]); % [left, bottom, width, height])
plot(time,uplot(:))
hold on
xlabel('years'); 
ylabel('THC'); 

Tbar=Tnow;
Sbar=Snow;
ubar=u;

var_THC=var(uplot(n/2:n));
var_T=var(Tplot(:,n/2:n)');
var_S=var(Splot(:,n/2:n)');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rho=linear_density(alpha,beta,T,S)

rho=-alpha*T+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 u=transport(u0,rho,delta)
  u=-u0*(rho(1)*delta+rho(3)*(1-delta)-rho(2));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dT_dt=time_rate_of_change_T ...
    (n,T,u,Tstar,gamma_T,V,H,sqrt_C,random_matrix,which_noise)
% rate of change of temperature:

my_noise=surface_noise(n,sqrt_C,random_matrix,which_noise);
if (u>=0)
  dT_dt=[my_noise(1)+gamma_T*(Tstar(1)-T(1)) + u*(T(3)-T(1))/V(1), ...
	my_noise(2)+gamma_T*(Tstar(2)-T(2)) + u*(T(1)-T(2))/V(2), ...
	u*(T(2)-T(3))/V(3)];
else
  dT_dt=[my_noise(1)+gamma_T*(Tstar(1)-T(1)) - u*(T(2)-T(1))/V(1), ...
	my_noise(2)+gamma_T*(Tstar(2)-T(2)) - u*(T(3)-T(2))/V(2), ...
	-u*(T(1)-T(3))/V(3)];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dS_dt=time_rate_of_change_S ...
    (n,S,u,Sstar,gamma_S,FW,sqrt_C,V,H,random_matrix,which_noise)
% rate of change of temperature:

my_noise=surface_noise(n,sqrt_C,random_matrix,which_noise);
if (u>=0)
  dS_dt=[ (FW(1)+my_noise(3) ...
	+gamma_S*(Sstar(1)-S(1))*V(1) + u*(S(3)-S(1)) )/V(1), ...
        (FW(2)+my_noise(4) ...
	+gamma_S*(Sstar(2)-S(2))*V(2) + u*(S(1)-S(2)) )/V(2), ...
	(u*(S(2)-S(3)))/V(3) ];
else
  dS_dt=[ (FW(1)+my_noise(3) ...
	+gamma_S*(Sstar(1)-S(1))*V(1) - u*(S(2)-S(1)))/V(1), ...
        (FW(2)+my_noise(4) ...
	+gamma_S*(Sstar(2)-S(2))*V(2) - u*(S(3)-S(2)))/V(2), ...
	(-u*(S(1)-S(3)))/V(3) ];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function my_noise=surface_noise(n,sqrt_C,random_matrix,which_noise)
% returns noise forcing for surface boxes: elements 1:2 are heat
% flux forcing for boxes 1:2; elements 4:5 are FW forcing for boxes
% 1:2. 

if strcmp(which_noise,'none')
  my_noise=[0 0 0 0 0 0]';
elseif strcmp(which_noise,'use_sqrt_C')
  my_noise=sqrt_C*random_matrix(:,n);
else
  which_noise
  error('unrecognized value of which_noise in function surface_noise.')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [L,H,delta,W,V,Area,Sv,year,u0,alpha,beta,gamma_T,gamma_S ...
      ,FW,Tstar,Sstar,delta_t,N_max,ubar,T,S,F ...
      ,random_matrix,noise_correlation_time]= ...
    set_parameters

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

% box geometry:
% north-south size:
L(1)=5000.0e5;
L(2)=1000.0e5;
L(3)=L(1);
% depth:
H(1)=1.0e5;
H(3)=4.0e5;
H(2)=H(1)+H(3);
delta= H(1)/(H(1)+H(3));
% width:
W=4000.0e5;
% volumes:
V(1)=L(1)*H(1)*W;
V(2)=L(2)*(H(1)+H(3))*W;
V(3)=L(3)*H(3)*W;

% areas:
Area(1)=L(1)*W;
Area(2)=L(2)*W;

% additional parameters:

Sv=1.e12;
year=3600.0*24.0*365.0;
% parameter for the relation between u (THC) and the density
% gradients: 
u0=16*Sv/0.0004;  % 16*Sv/0.0005;
% parameters for linearized equation of state:
alpha=1668.0e-7;
beta=0.781e-3;
% set SST restoring time in days per 50 meters upper layer
% thickness:
days_per_50m_mixed_layer=50.0;
gamma_T=1.0/(300.0*86400.0); %*dbgz:   1.0/((days_per_50m_mixed_layer*86400.0)*(H(1)/50.e2));
gamma_S=0.0/((days_per_50m_mixed_layer*86400.0)*(H(1)/50.e2));
FW(1)=1.2*(100.0/year)*35.0*Area(1);
FW(2)=-FW(1);

% the temperature to which the model surface temperature is restored:
Tstar=[22.0,0.0]; 
% the salinity to which the model surface salinity is restored (only
% if gamma_S is not zero):
Sstar=[36.0,34.0];

% parameters for numerically solving the steady state problem:
noise_correlation_time=86400.0*15.0;
% time step:
delta_t=5*86400.0;
% integration time:
integration_time=1500.0*year;
% number of time steps:
N_max=floor((integration_time/delta_t)/100)*100;

% steady state initial conditions:
ubar= 20*Sv;
T(1)= 20.0;
T(2)= 1.0;
T(3)= 1.0;
S(1)= 35.5;
S(2)= 34.5;
S(3)= 34.5;

% uncomment the following lines to test the density calculation:
% pressure=1.0;
% for i=1:3
%   rho_linear(i)=linear_density(alpha,beta,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


% an operator that projects the noise forcing on the surface
% temperature and salinity equations only:
F=zeros(6);
F(1,1)=1.0;
F(2,2)=1.0;
F(4,4)=1.0;
F(5,5)=1.0;

% 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_random_matrix=rand(6,N_max);% random('Normal',0.0,1.0,6,N_max);
% next, create a first order markov process with a one week
% correlation time from the above random numbers:
random_matrix=zeros(6,N_max);
random_matrix(:,1)=tmp_random_matrix(:,1);
R_noise=exp(-delta_t/noise_correlation_time);
tmp_R2=sqrt(1-R_noise^2);
for n=1:N_max-1
  random_matrix(:,n+1) = ...
      R_noise*random_matrix(:,n) + tmp_R2*tmp_random_matrix(:,n);
end
clear tmp_R2 tmp_random_matrix

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
