function doall=nonormal_THC
% matlab program to look at non normal growth in THC box model.
clear all; close all
global A R X Z C

% 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,R,R_T,R_S,X,Y,F ...
      ,random_matrix,noise_correlation_time,epsilon_ic] ...
    = set_parameters;

% solve the steady state problem:
disp(' running the box model with no noise to steady state:')
sqrt_C=zeros(6);
which_noise='none';
optimal_noise_vector=[0 0 0 0 0 0]';
tic;
[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,optimal_noise_vector,which_noise);
toc;

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

disp(' variance of solution during second half of model run is:')
var_T
var_S
var_THC

disp(' Projection operator R used to define norm Kernel X: R*[T,S]=')
THC=R*[Tbar(1:3) Sbar(1:3)]'
disp(' Projection operators R_T,R_S: R_(T,S)*[T,S]=')
THC_T=R_T*[Tbar(1:3) Sbar(1:3)]', THC_S=R_S*[Tbar(1:3) Sbar(1:3)]'
disp(' Checking Kernel operator X: sqrt([T,S]*X*[T,S]^T)=')
sqrt([Tbar(1:3) Sbar(1:3)]*X*[Tbar(1:3) Sbar(1:3)]')

% error('enough for now')

% evaluate linearized model matrix analytically:
A = linearized_model_matrix ...
    (Tbar,Sbar,ubar,u0,alpha,beta,gamma_T,delta,V);

% evaluate linearized model matrix by finite difference:
which_noise='none';
A_fd = linearized_model_matrix_finite_difference ...
    (Tbar,Sbar,ubar,u0,alpha,beta,gamma_T,gamma_S ...
    ,delta,Tstar,Sstar,FW,Area,year,sqrt_C,V,H ...
    ,random_matrix,optimal_noise_vector,which_noise);

% compare the two linear matrices:
disp(' Comparing the analytical and finite difference values of A:')
A
A_fd;
A-A_fd;

% eigenvalues and eigenvectors of A:
disp(' eigen vals and vects of A:')
[Evec,Eval]=eig(A);
evec_A=Evec
diag_eval=diag(Eval)*year
disp(' 1/eigen vals (year):')
for i=1:6
  if imag(Eval(i,i))==0
    inverse_eigenvalues(i)=complex(1.0/(year*real(Eval(i,i))),0.0);
  else
    inverse_eigenvalues(i)=complex(1.0/(year*real(Eval(i,i))),1.0/(year*imag(Eval(i,i))));
  end    
  disp(sprintf('%10.4f + i %10.4f',real(inverse_eigenvalues(i)),imag(inverse_eigenvalues(i))))
end

% optimal initial conditions
optimal_initial_conditions ... 
    (A,evec_A,R,R_T,R_S,X,Y,F,ubar,Tbar,Sbar,epsilon_ic,u0,alpha,beta ...
    ,gamma_T,gamma_S,delta,year,FW,delta_t,Tstar,Sstar ...
    ,V,H,Area,sqrt_C,random_matrix,optimal_noise_vector,which_noise);

% Stochastic optimals:
disp('Calculating stochastic optimals matrix Z:')
calculate_Z=1
tic;
 [Z,Stochastic_optimals_values,Stochastic_optimals_vectors,noptimal_vector,C] = ...
    Stochastic_optimals(A,X,Y,F,year,delta_t,noise_correlation_time,Area,calculate_Z);
toc;

%dbgz:
%disp('dbgz: setting C temporarily:')
%FW_noise_var=( (20.0/year)*35.0*Area(1) )^2;
%% heat flux needs to be in units of deg/sec, or in the same units
%% of (\gamma*\delta T).  use the equation 
%% \delta H=C_p \rho_0\Delta Z \gamma\delta T to find:
%% \gamma\delta T=\delta H/(C_p \rho_0\Delta Z)
%%               = (10watts/m^2)/([4000J/(kg*deg)]*[1000kg/m^3]*[500m]):
%H_noise_var=( 10.0/(4000.0*1000.0*500.0) )^2;
%C=diag([H_noise_var H_noise_var 0.0 FW_noise_var FW_noise_var 0.0])

% run the model with stochastic forcing:
disp(' running the box model with stochastic noise:')
tic;
disp('sqrt_C:')
sqrt_C=C
if sqrt_C*sqrt_C ~= C 
  error('sqrt_C*sqrt_C ~=C')
end
which_noise='use_optimal_vector';
noise_amplitude=1.e-11;
optimal_noise_vector=Stochastic_optimals_vectors(:,6)*noise_amplitude
[var_T,var_S,var_THC,Tbar_stoch,Sbar_stoch,ubar_stoch,Tplot1,Splot1,uplot1]= ...
    box_model(Tbar,Sbar,u0,alpha,beta ...
    ,gamma_T,gamma_S,delta,year,FW,delta_t,N_max,Tstar,Sstar ...
    ,V,H,Area,sqrt_C,random_matrix,optimal_noise_vector,which_noise);
clear Tplot1 Splot1 uplot1
toc;
disp(' variance of solution during second half of model run is:')
var_T
var_S
var_THC


% run a linearized model with stochastic forcing:
disp('running linearized box model with stochastic forcing:')
tic;
T=[0 0 0]; S=[0 0 0];
which_noise='use_optimal_vector';
fid_timeseries = fopen('timeseries.data','wb');
fwrite(fid_timeseries,N_max,'integer*8');
fwrite(fid_timeseries,delta_t,'real*8');
for i=6:-1:1
  C=Stochastic_optimals_vectors(:,i)*Stochastic_optimals_vectors(:,i)';
  linear_model_title=['stochastically forced linear model, optimal vector # ',sprintf('%d',i)]
  optimal_noise_vector=Stochastic_optimals_vectors(:,i)*noise_amplitude
  N_plot=2000.0*year/delta_t;
  [var_T,var_S,var_THC,Tplot,Splot,uplot]=linearized_box_model ...
      (T,S,A,evec_A,F,R,R_T,R_S,sqrt_C,Area,V,year,delta_t,alpha,beta...
      ,N_max,random_matrix,optimal_noise_vector,which_noise,N_plot,linear_model_title);
  var_T,var_S,var_THC
  disp(['optimal noise vector #=',sprintf('%d',i) ...
	, '; expected THC variance: Tr(Z*C)=',sprintf('%d',trace(Z*C)*noise_amplitude^2) ...
	,'; linearized model variance=' ,sprintf('%d',var_THC) ...
	,'; ratio=' ,sprintf('%d',var_THC/(trace(Z*C)*noise_amplitude^2))])
  fwrite(fid_timeseries,uplot,'real*8');
end
toc;

% 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,optimal_noise_vector,which_noise)

% 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=density(alpha,beta,Tnow,Snow);
  
  % 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,optimal_noise_vector,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,optimal_noise_vector,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

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')
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')
ylabel('S'); 
axes('position', [0.1,0.1,0.8,0.2]); % [left, bottom, width, height])
plot(time,uplot(:))
hold on
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=density(alpha,beta,T,S)

rho=-alpha*T+beta*S;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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,optimal_noise_vector,which_noise)
% rate of change of temperature:

my_noise=surface_noise(n,sqrt_C,random_matrix,optimal_noise_vector,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,optimal_noise_vector,which_noise)
% rate of change of temperature:

my_noise=surface_noise(n,sqrt_C,random_matrix,optimal_noise_vector,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 A = linearized_model_matrix ...
    (Tbar,Sbar,ubar,u0,alpha,beta,gamma_T,delta,V)
% The linearized model matrix using the analytic expression for the
% linearized model equations: (assuming thermally dominant flow!!!)

if (ubar<0)
  sprintf('*** Error: ubar is negative...: ubar= %d',ubar)
  error('stopping...')
end

sprintf('need to fix bug with third element of dS_Dt that had -u instead of +u')
sprintf('see file BUG in thid directory')
error('stopping...')
    
A= [ (-ubar+u0*delta*alpha*(Tbar(3)-Tbar(1)))/V(1)-gamma_T ...
    (-u0*alpha*(Tbar(3)-Tbar(1)))/V(1) ...
    (ubar+u0*(1-delta)*alpha*(Tbar(3)-Tbar(1)))/V(1) ...
    (-u0*delta*beta*(Tbar(3)-Tbar(1)))/V(1) ...
    (u0*beta*(Tbar(3)-Tbar(1)))/V(1) ...
    (-u0*(1-delta)*beta*(Tbar(3)-Tbar(1)))/V(1) ,                     
    ...
    (ubar+u0*delta*alpha*(Tbar(1)-Tbar(2)))/V(2) ...
    (-ubar-u0*alpha*(Tbar(1)-Tbar(2)))/V(2)-gamma_T ...
    (u0*(1-delta)*alpha*(Tbar(1)-Tbar(2)))/V(2) ...
    (-u0*delta*beta*(Tbar(1)-Tbar(2)))/V(2) ...
    (u0*beta*(Tbar(1)-Tbar(2)))/V(2) ...
    (-u0*(1-delta)*beta*(Tbar(1)-Tbar(2)))/V(2) ,
    ...
    (u0*delta*alpha*(Tbar(2)-Tbar(3)))/V(3) ...
    (ubar-u0*alpha*(Tbar(2)-Tbar(3)))/V(3) ...
    (-ubar+u0*(1-delta)*alpha*(Tbar(2)-Tbar(3)))/V(3) ...
    (-u0*delta*beta*(Tbar(2)-Tbar(3)))/V(3) ...
    (u0*beta*(Tbar(2)-Tbar(3)))/V(3) ...
    (-u0*(1-delta)*beta*(Tbar(2)-Tbar(3)))/V(3) ,
    ...
    (u0*delta*alpha*(Sbar(3)-Sbar(1)))/V(1) ...
    (-u0*alpha*(Sbar(3)-Sbar(1)))/V(1) ...
    (u0*(1-delta)*alpha*(Sbar(3)-Sbar(1)))/V(1) ...
    (-ubar-u0*delta*beta*(Sbar(3)-Sbar(1)))/V(1) ...
    (u0*beta*(Sbar(3)-Sbar(1)))/V(1) ...
    (ubar-u0*(1-delta)*beta*(Sbar(3)-Sbar(1)))/V(1) ,
    ...
    (u0*delta*alpha*(Sbar(1)-Sbar(2)))/V(2) ...
    (-u0*alpha*(Sbar(1)-Sbar(2)))/V(2) ...
    (u0*(1-delta)*alpha*(Sbar(1)-Sbar(2)))/V(2) ...
    (ubar-u0*delta*beta*(Sbar(1)-Sbar(2)))/V(2) ...
    (-ubar+u0*beta*(Sbar(1)-Sbar(2)))/V(2) ...
    (-u0*(1-delta)*beta*(Sbar(1)-Sbar(2)))/V(2) ,
    ...
    (u0*delta*alpha*(Sbar(2)-Sbar(3)))/V(3) ...
    (-u0*alpha*(Sbar(2)-Sbar(3)))/V(3) ...
    (u0*(1-delta)*alpha*(Sbar(2)-Sbar(3)))/V(3) ...
    (-u0*delta*beta*(Sbar(2)-Sbar(3)))/V(3) ...
    (ubar+u0*beta*(Sbar(2)-Sbar(3)))/V(3) ...
    (-ubar-u0*(1-delta)*beta*(Sbar(2)-Sbar(3)))/V(3) ];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function A_fd = linearized_model_matrix_finite_difference ...
  (Tbar,Sbar,ubar,u0,alpha,beta,gamma_T,gamma_S ...
  ,delta,Tstar,Sstar,FW,Area,year,sqrt_C,V,H ...
  ,random_matrix,optimal_noise_vector,which_noise);
% The linearized model matrix using a finite difference
% approximation for the linearized model equations:

% parameters: the size of the linear perturbation used to obtain
% the finite difference approximation for the model matrix:
epsilon=0.0001;

for j=1:6
  % setup a vector epsilon that is different from zero only in the
  % j-direction: 
  epsilon_j(1:6)=0.0;
  epsilon_j(j)=epsilon;
  T=Tbar+epsilon_j(1:3);
  S=Sbar+epsilon_j(4:6);
  rho=density(alpha,beta,T,S);
  u=transport(u0,rho,delta);
  % artificially set time step to 1, required by surface_noise that is
  % called by time_rate_of_change_S:
  n=1;
  A_fd(j,1:3) = ...
      (time_rate_of_change_T(n,T,u,Sstar,gamma_T,V,H,sqrt_C,random_matrix,optimal_noise_vector,which_noise) ...
      -time_rate_of_change_T(n,Tbar,ubar,Sstar,gamma_T,V,H,sqrt_C,random_matrix,optimal_noise_vector,which_noise))/epsilon;
  A_fd(j,4:6) = ...
      (time_rate_of_change_S ...
      (n,S,   u,   Sstar,gamma_S,FW,sqrt_C,V,H,random_matrix,optimal_noise_vector,which_noise) ...
      -time_rate_of_change_S ...
      (n,Sbar,ubar,Sstar,gamma_S,FW,sqrt_C,V,H,random_matrix,optimal_noise_vector,which_noise)) ...
      /epsilon;
end

% seems we actually need the transposed of the above:
A_fd=A_fd';


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function optimal_initial_conditions ...
    (A,evec_A,R,R_T,R_S,X,Y,F,ubar,Tbar,Sbar,epsilon_ic,u0,alpha,beta ...
    ,gamma_T,gamma_S,delta,year,FW,delta_t,Tstar,Sstar ...
    ,V,H,Area,sqrt_C,random_matrix,optimal_noise_vector,which_noise)
% optimal initial conditions that maximize initial rate of growth:
disp(' optimal i.c. that maximize \dot A(T=0):')
[evec_ATA,eval_ATA]=eig(A'*A,Y)

% optimal initial conditions that maximize growth after N years:
disp(' Optimal growth for a given number of years:')
for N=[0:0.1:2.9 3:1:10 12:2:20 25:5:50 60:5:100 ]; % 200:100:1000];
  B=expm(N*86400.0*365.0*A);
  [Evec,Eval]=eig(B'*X*B,Y);
  for i=1:6
    if abs(R*Evec(:,i)) < 1.e-10
      if abs(R*B*Evec(:,i)) < 1.e-10
	THC_amplif(i)=0.0;
      else
	THC_amplif(i)=Inf;
      end
    else
      THC_amplif(i)=(R*B*Evec(:,i))/(R*Evec(:,i));
    end
  end
  % Sort eigenvalues reflecting amplification:
  [sorted_tmp,Indx]=sort(abs(diag(Eval)));
  for i=1:6; sorted_eigenvals(i)=Eval(Indx(i),Indx(i)); end;
  for i=1:6; sorted_THC_amplif(i)=THC_amplif(Indx(i)); end;
  disp(['N=',sprintf('%6.1f',N) '; eigenvals:  ' sprintf('%14.6g ',sorted_eigenvals)])
  disp(['N=',sprintf('%6.1f',N) '; THC amplif: ' sprintf('%14.6g ',sorted_THC_amplif(:))])
  initial_THC=0.0*ubar/1.0e12+epsilon_ic*(R*  Evec(:,Indx(6)));
  final_THC  =0.0*ubar/1.0e12+epsilon_ic*(R*B*Evec(:,Indx(6)));
  disp(['N=',sprintf('%6.1f',N) '; init&final THC for eigenvec #' sprintf('%d',i) ...
	'=' sprintf(' %d ',initial_THC,final_THC)])
  % plot model evolution for two of the optimal initial conditions,
  % at the following two times:
  time_mechanism1=2.0;
  time_mechanism2=45.0;
  if ( abs(N-time_mechanism1)<0.0001 | abs(N-time_mechanism2)<0.0001 )
    % number of time steps:
    if abs(N-time_mechanism1)<0.0001
      N_steps=floor((15*year/delta_t)/5)*5;
    else
      N_steps=floor((200*year/delta_t)/10)*10;
    end
    for i=1:6
      disp(['analyzing mode #: '  sprintf(' %d ',i)])
      initial_THC=0.0*ubar/1.0e12+epsilon_ic*(R*  Evec(:,Indx(i)));
      final_THC  =0.0*ubar/1.0e12+epsilon_ic*(R*B*Evec(:,Indx(i)));
      disp([' init&final THC for eigenvec #' sprintf('%d',i) ...
	'=' sprintf(' %d ',initial_THC,final_THC)])
% $$$       disp(' running the full box model with no noise and from optimal i.c.:')
% $$$       % optimal initial conditions: small perturbation along the eigenvectors:
% $$$       T=Tbar+epsilon_ic*Evec(1:3,Indx(i))';
% $$$       S=Sbar+epsilon_ic*Evec(4:6,Indx(i))';
% $$$       sqrt_C=zeros(6);
% $$$       [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_steps,Tstar,Sstar ...
% $$$ 	  ,V,H,Area,sqrt_C,random_matrix,optimal_noise_vector,which_noise);
      disp(' running the linearized box model with no noise and from i.c=eigenvec:')
      T=epsilon_ic*Evec(1:3,Indx(i))';
      S=epsilon_ic*Evec(4:6,Indx(i))';
      linear_model_title=['linear model; optimal i.c.=eigenvec #' ...
	    sprintf('%d',i) '; N=' sprintf('%d',N)];
      N_plot=N_steps;
      [var_T,var_S,var_THC,Tplot,Splot,uplot]=linearized_box_model ...
	  (T,S,A,evec_A,F,R,R_T,R_S,sqrt_C,Area,V,year,delta_t,alpha,beta...
	  ,N_steps,random_matrix,optimal_noise_vector,which_noise,N_plot,linear_model_title);
      disp(' running the linearized box model with no noise and from i.c=eigenvec:')
    end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ...
    [Z,Stochastic_optimals_values,Stochastic_optimals_vectors,noptimal_vector,C] = ...
    Stochastic_optimals(A,X,Y,F,year,delta_t,noise_correlation_time,Area,calculate_Z) 

% maximum integration time for evaluating statistically steady
% stochastic optimals:
tau=1000.0*year;
% delta t for evaluating the integral expression for Z:
delta_t_Z=delta_t;
% integration increments when using the transformed integration
% varibles:
delta_q=delta_t;
delta_r=delta_t;
% size of truncated integration area:
T_trunc=noise_correlation_time*10.0;
% initialize Z:
Z=zeros(6);

% transform to new variables:
% r=t-s; q=t+s;   such that   s=(q-r)/2; t=(q+r)/2;
% note that the transformation is
%   [r q] = M [t s]
%   where M = [1 -1, 1 1]
% so that dt ds = det(M^{-1}) dr dq = (1/2) dr dq.
det_M=2.0;

% for efficiency, truncate integration region in the (r,q) plane based
% on the short correlation time of the noise as expressed in the D
% factor:
delta_qr=(delta_q*delta_r/det_M);


if calculate_Z == 1
  for q=0:delta_q:tau
    for r=max(-q,-T_trunc):delta_r:min(q,T_trunc)
      s=(q-r)/2.0;
      t=(q+r)/2.0;
      B1=expm(s*A);
      B2=expm(t*A);
      D_noise=exp(-abs(t-s)/noise_correlation_time);
      Z=Z+delta_qr*D_noise*F'*B1'*X*B2*F;
    end
  end
  for q=tau:delta_q:tau*2
    for r=max(q-2*tau,-T_trunc):delta_r:min(2*tau-q,T_trunc)
      s=(q-r)/2.0;
      t=(q+r)/2.0;
      B1=expm(s*A);
      B2=expm(t*A);
      D_noise=exp(-abs(t-s)/noise_correlation_time);
      Z=Z+delta_qr*D_noise*F'*B1'*X*B2*F;
    end
  end
else
Z =  1.0e+18 * ...
 [ 0.00004651338908  -0.00024591138238                  0  -0.00019703735556     0.00415063364221                  0
  -0.00024591138238   0.00132148905850                  0   0.00408046479873    -0.01941085272883                  0
                  0                  0                  0                  0                    0                  0
  -0.00019703735556   0.00408046479873                  0   0.66155223233575     0.63110403818924                  0
   0.00415063364221  -0.01941085272883                  0   0.63110403818924     2.35889216472359                  0
                  0                  0                  0                  0                    0                  0];
end
format long

[Stochastic_optimals_vectors,Stochastic_optimals_values]=eig(Z);

disp('Z=')
Z
disp('stochastic optimals eigen values')
disp(diag(Stochastic_optimals_values))
disp('stochastic optimals eigen vectors:')
disp(Stochastic_optimals_vectors)

% look at the eigenvalues of Z to decide which one is the optimal
% and set it here manually:
[Y,Index]=sort(diag(Stochastic_optimals_values));
noptimal_vector=Index(6);
disp(['choosing eigenvector # ' sprintf('%d',noptimal_vector) ' of Z as the optimal.'])

% The optimal stochastic forcing covariance matrix, 
disp(' optimal noise correlation:')
C=   Stochastic_optimals_vectors(:,noptimal_vector) ...
    *Stochastic_optimals_vectors(:,noptimal_vector)'...
    *( 1.0 )^2

disp(' the expected optimal variance from trace(Z*C):')
var_THC=trace(Z*C)

format

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [var_T,var_S,var_THC,Tplot,Splot,uplot]=linearized_box_model ...
    (T,S,A,evec_A,F,R,R_T,R_S,sqrt_C,Area,V,year,delta_t,alpha,beta...
    ,N_max,random_matrix,optimal_noise_vector,which_noise,N_plot,linear_model_title);

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

% get initial conditions from input:
P=zeros(6,1);
P(1:3)=T(1:3);
P(4:6)=S(1:3);
P0=P;

my_noise=zeros(6,1);
for n=1:N_max
  my_noise=surface_noise(n,sqrt_C,random_matrix,optimal_noise_vector,which_noise);
  P=P+delta_t*(A*P+F*my_noise);
  Tplot(:,n)=P(1:3);
  Splot(:,n)=P(4:6);
  uplot(n)=R*P;
  uplot_S(n)=R_S*P;
  uplot_T(n)=R_T*P;
  %*dbgz: I am not sure this is correct:
  modes_projection(:,n)=evec_A'*P;
end

var_THC=var(uplot);
var_T=var(Tplot');
var_S=var(Splot');

h=figure;
NP=N_plot;
time=(1:NP)*delta_t/year;
set(h,'Position', [500 520 400 400]);
axes('position', [0.1,0.7,0.8,0.25]); % [left, bottom, width, height])
title(linear_model_title)
hold on
plot(time,alpha*Tplot(1,1:NP),'k-',time,alpha*Tplot(2,1:NP),'r--',time,alpha*Tplot(3,1:NP),'b-.')
set(gca,'XTickLabel','')
legend('\alpha T_1','\alpha T_2','\alpha T_3')
grid
ylabel('\alpha T'); 
axes('position', [0.1,0.4,0.8,0.25]); % [left, bottom, width, height])
plot(time,beta*Splot(1,1:NP),'k-',time,beta*Splot(2,1:NP),'r--',time,beta*Splot(3,1:NP),'b-.')
set(gca,'XTickLabel','')
hold on
grid
legend('\beta S_1','\beta S_2','\beta S_3')
ylabel('\beta S'); 
axes('position', [0.1,0.1,0.8,0.25]); % [left, bottom, width, height])
plot(time,uplot(1:NP),'k-',time,uplot_T(1:NP),'r--',time,uplot_S(1:NP),'b-.')
legend('U','U_T','U_S')
hold on
grid
xlabel('time (years)')
ylabel('THC'); 

figure
title(linear_model_title)
hold on
plot(time,modes_projection(1,1:NP),time,modes_projection(2,1:NP)...
    ,time,modes_projection(3,1:NP),time,modes_projection(4,1:NP)...
    ,time,modes_projection(5,1:NP),time,modes_projection(6,1:NP));
legend('1','2','3','4','5','6');
xlabel('time (years)');
ylabel('projection on modes'); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function my_noise=surface_noise(n,sqrt_C,random_matrix,optimal_noise_vector,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);
elseif strcmp(which_noise,'use_optimal_vector')
  my_noise=random_matrix(1,n)*optimal_noise_vector;
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,R,R_T,R_S,X,Y,F ...
      ,random_matrix,noise_correlation_time,epsilon_ic]= ...
    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;

% parameters for linearized model matrix:

Sv=1.e12;
year=3600.0*24.0*365.0;
u0=16*Sv/0.0004;  % 16*Sv/0.0005;
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.0*(100.0/year)*35.0*Area(1);
FW(2)=-FW(1);
Tstar=[22.0,0.0];
Sstar=[36.0,34.0];
% epsilon_ic multiplies the eigenvector of fastest growing mode,
% before using it as initial conditions in the short mode runs:
epsilon_ic=-0.0001;


% 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=5000.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;

% a norm kernel that is just an energy norm (sum of squares)
% weighted by the respective variance during a model run under
% stochastic forcing.  each element here is 1/variance:
Y=diag(1.0e3*[1/0.03268    1/0.00794    1/0.00140  1/0.1417    1/0.1286    1/0.0878]);

% A vector used to form a norm Kernel that measures amplitude of THC
% (the transport in the model, in Sv):

R=[u0*delta*alpha ...
  -u0*alpha ...
  u0*(1-delta)*alpha ...
  -u0*delta*beta ...
  u0*beta ...
  -u0*(1-delta)*beta]/1.0e12;

% The thermohaline circulation due to the Salinity/ temperature
% effects alone is given by R_S*P, R_T*P:
R_S=R*diag([0 0 0 1 1 1]);
R_T=R*diag([1 1 1 0 0 0]);

% A vector that measures the temperature gradient (T_1-T_2).
% in conjunction with R above, this is used to form a Kernel that
% measures the meridional heat transport = U*(T_1-T_2):
R1=[1 -1 0 0 0 0];

for i=1:6
  for j=1:6
    X(i,j) =R(i)*R(j);
%    X1(i,j)=R(i)*R1(j);
  end
end

% regularize the norm kernel using Petros' suggestion:
epsilon_regularize=1.0e-8;
X=X+epsilon_regularize*eye(6);

eigs_x=eig(X)
% eigs_x1=eig(X1)

%solving using a singular Y that is the same norm kernel as X:
 Y=X;

% 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=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

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