function doall=box_geochem

% A 3 box THC model with geochemistry.  Written as a practice problem
% to study basic geochemistry.  See documantation in enclosed latex
% notes, file paleo_co2_tutorial.tex.  
% Eli Tziperman

clear all; close all

% Solution of the full model:
global Tplot1 Splot1 uplot1

% profile on

%   Plotted circulation corresponds to a positive THC:
%
%   --------------------------------------
%   |			       	       	 |
%   |		atmosphere     	       	 |
%   |		               	       	 |
%   --------------------------------------
%   --------------------------------------
%   |  		|	    	      	 |
%   |  1      ======> 		      2	 |
%   |	  ^  	|	     | 	       	 |
%   |     |     |------------|-----------|
%   |-----|------	     | 	       	 |
%   |     |  		     v           |
%   |			       	       	 |
%   |			       	       	 |
%   |		               	       	 |
%   |		 	    	     3 	 |
%   |			    	       	 |
%   --------------------------------------
%  S. Pole                             Equator
%


% set model parameters and initial conditions:
[L,H,W,V,Area,Sv,year,u0,alpha,beta,gamma_T,gamma_S ...
 ,k_h,k_v_highlat,k_v_lowlat,PV ...
 ,FW,Tstar,Sstar,delta_t,N_max,ubar,T,S,ALK,TCO2,PO4,Hp ...
 ,Redfield_PC,Redfield_PN,hPO4,light_factor,rPO4 ...
 ,solve_physics,solve_geochemistry] ...
    =set_parameters;

% run the model to a steady state:
% run the box model:
[var_T,var_S,var_THC,Tbar,Sbar,ubar...
 ,ALKbar,PO4bar,TCO2bar,CO2_atm_bar ...
 ,Tplot1,Splot1,uplot1]= ...
    box_model(T,S,ALK,TCO2,PO4,Hp,u0,alpha,beta ...
	      ,k_h,k_v_highlat,k_v_lowlat,PV ...
	      ,gamma_T,gamma_S,year,FW,delta_t,N_max,Tstar,Sstar ...
	      ,V,H,L,Area,Redfield_PC,Redfield_PN,hPO4,light_factor,rPO4 ...
	      ,solve_physics,solve_geochemistry);
%print -depsc figure1.eps

disp(' Steady state solution is:')
Tbar
Sbar
ubar
ALKbar
PO4bar
TCO2bar
CO2_atm_bar

% profile report


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [var_T,var_S,var_THC,Tbar,Sbar,ubar...
	  ,ALKbar,PO4bar,TCO2bar,CO2_atm_bar ...
	  ,Tplot,Splot,uplot] ...
    =box_model (T,S,ALK,TCO2,PO4,Hp,u0,alpha,beta ...
		,k_h,k_v_highlat,k_v_lowlat,PV ...
		,gamma_T,gamma_S,year,FW,delta_t,N_max,Tstar,Sstar ...
		,V,H,L,Area,Redfield_PC,Redfield_PN,hPO4,light_factor,rPO4 ...
		,solve_physics,solve_geochemistry)

% Initialize:
ALKnew=ALK;   ALKnow=ALK;   ALKold=ALK;
PO4new=PO4;   PO4now=PO4;   PO4old=PO4;
Snew=S;       Snow=S;       Sold=S;
TCO2new=TCO2; TCO2now=TCO2; TCO2old=TCO2; 
Tnew=T;       Tnow=T;       Told=T;

% Initialize arrays containing plot data:
Tplot   =zeros(3,N_max);
Splot   =zeros(3,N_max);
ALKplot =zeros(3,N_max);
PO4plot =zeros(3,N_max);
TCO2plot=zeros(3,N_max);
uplot   =zeros(1,N_max);

% some silly consmetics stuff for indicating progress of run:
status_bar='[                     ]';
i_status_bar=1;
fprintf(1,'%s',['running the box model to steady state: ',status_bar])

% main time loop:
for n=1:N_max

  % some silly consmetics stuff for indicating progress of run:
  N_print=floor(N_max/20);
  if n == floor(n/N_print)*N_print 
    i_status_bar=i_status_bar+1;
    for ii=1:i_status_bar status_bar(ii+1)='.'; end
    % Delete old status bar and print updated one:
    fprintf(1,'\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b%s',status_bar)
  end
  if n == N_max
    disp(' ')
  end

  % Calculate densities:
  % rho=linear_density(alpha,beta,Tnow,Snow);
  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);
  
  if (solve_geochemistry)
  % Solve for carbonate system:
    Hp_init_guess=Hp;
    [CO2_g,H2CO3,OHm,Hp,pH,HCO3m,CO3m2,B_OH_3,B_OH_4m,B_T]= ...
	solve_carbonate_system(ALKnow,TCO2now,Tnow,Snow,Hp_init_guess);
    if n == N_max
      % print details of carbonate system on final time step:
      disp('The solution of carbonate system on final time step is')
      ALKnow,TCO2now,Tnow,Snow,Hp_init_guess
     [CO2_g,H2CO3,OHm,Hp,pH,HCO3m,CO3m2,B_OH_3,B_OH_4m,B_T]= ...
	solve_carbonate_system(ALKnow,TCO2now,Tnow,Snow,Hp_init_guess)
      disp('End of carbonate system solution.')
      disp('export production at end of run:')
      EP1=export_production(PO4now(1),hPO4,light_factor(1),rPO4(1));
      EP2=export_production(PO4now(2),hPO4,light_factor(2),rPO4(2));
      % our units are PO4/cm^2/sec.  Need to convert to gr C/year and to m^2:
      export_production_gr_C_per_m2=...
	  [EP1*(122+Rain(Tnow(1))),EP2*(122+Rain(Tnow(2)))] ...
	  *(12.011*10^-6)*(365*86400)*100
      global_export_production_Gt_C_per_yr=...
	  [EP1*(122+Rain(Tnow(1)))*Area(1),EP2*(122+Rain(Tnow(2)))*Area(2)] ...
	  *(12.011*10^-6)*(365*86400)*1.e-15
   end
  % solve for atmospheric co2 concentration:
  CO2_atm=(Area(1)*CO2_g(1)+Area(2)*CO2_g(2))/(Area(1)+Area(2));
  end

  % Leap frog time stepping:

  if (solve_physics)
    Tnew=Told+2*delta_t*time_rate_of_change_T ...
	 (n,Tnow,u,Tstar,gamma_T,V,Area,H,L,k_h,k_v_highlat,k_v_lowlat);
    Snew=Sold+2*delta_t*time_rate_of_change_S ...
	 (n,Snow,u,Sstar,gamma_S,FW,V,Area,H,L,k_h,k_v_highlat,k_v_lowlat);
  end
  if (solve_geochemistry)
    TCO2new=TCO2old+2*delta_t*time_rate_of_change_TCO2 ...
	    (TCO2now,PO4now,Tnow,u,V,Area,H,L,k_h,k_v_highlat,k_v_lowlat,PV,CO2_atm,CO2_g ...
	     ,Redfield_PC,Redfield_PN,hPO4,light_factor,rPO4);
    ALKnew=ALKold+2*delta_t*time_rate_of_change_ALK ...
	   (ALKnow,PO4now,Tnow,u,V,Area,H,L,k_h,k_v_highlat,k_v_lowlat ...
	     ,Redfield_PC,Redfield_PN,hPO4,light_factor,rPO4);
    PO4new=PO4old+2*delta_t*time_rate_of_change_PO4 ...
	   (PO4now,u,V,Area,H,L,k_h,k_v_highlat,k_v_lowlat ...
	     ,Redfield_PC,Redfield_PN,hPO4,light_factor,rPO4);
  end

  % Robert filter:
  robert_filter_coeff=0.25;
  if (solve_physics)
    Tnow    = Tnow    + robert_filter_coeff * ( Tnew - 2*Tnow + Told );
    Snow    = Snow    + robert_filter_coeff * ( Snew - 2*Snow + Sold );
  end
  if (solve_geochemistry)
    ALKnow  = ALKnow  + robert_filter_coeff * ( ALKnew - 2*ALKnow + ALKold );
    PO4now  = PO4now  + robert_filter_coeff * ( PO4new - 2*PO4now + PO4old );
    TCO2now = TCO2now + robert_filter_coeff * ( TCO2new - 2*TCO2now + TCO2old );
  end

  if (solve_physics)
    Told=Tnow;
    Tnow=Tnew;
    Sold=Snow;
    Snow=Snew;
  end
  if (solve_geochemistry)
    ALKold=ALKnow;
    ALKnow=ALKnew;
    PO4old=PO4now;
    PO4now=PO4new;
    TCO2old=TCO2now;
    TCO2now=TCO2new;
  end
 

  % plotting data:
  Tplot(:,n)=Tnow';
  Splot(:,n)=Snow';
  ALKplot(:,n)=ALKnow';
  PO4plot(:,n)=PO4now';
  TCO2plot(:,n)=TCO2now';
  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'); 

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,ALKplot(1,:),time,ALKplot(2,:),time,ALKplot(3,:))
hold on
legend('ALK_1','ALK_2','ALK_3')
ylabel('ALK'); 
axes('position', [0.1,0.4,0.8,0.2]); % [left, bottom, width, height])
plot(time,TCO2plot(1,:),time,TCO2plot(2,:),time,TCO2plot(3,:))
hold on
legend('TCO2_1','TCO2_2','TCO2_3')
ylabel('TCO2'); 
axes('position', [0.1,0.1,0.8,0.2]); % [left, bottom, width, height])
plot(time,PO4plot(1,:),time,PO4plot(2,:),time,PO4plot(3,:))
legend('PO4_1','PO4_2','PO4_3')
hold on
ylabel('PO4'); 

Tbar=Tnow;
Sbar=Snow;
ALKbar=ALKnow;
PO4bar=PO4now;
TCO2bar=TCO2now;
ubar=u;
CO2_atm_bar=CO2_atm;

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function u=transport(u0,rho)
  u=-20.0e12;
  %dbgz: need to make this a function of meridional temp gradient...


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function EP=export_production(PO4,hPO4,light_factor,rPO4)
% In units of micro mol PO4 per cm^2 per sec:
EP=rPO4*light_factor*PO4^2/(hPO4+PO4);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Rain_ratio=Rain(Temperature)
  Rain_ratio=61.0/(1+exp(0.1*(10-Temperature)));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dT_dt=time_rate_of_change_T ...
    (n,T,u,Tstar,gamma_T,V,Area,H,L,k_h,k_v_highlat,k_v_lowlat)
% rate of change of temperature:

% beware of Matlab mess: make sure each of the elements of the vector
% dX_dt are sorrounded by ().  Otherwise continuation character at end
% of line makes matlab think that each line is a different vector
% element:
if (u>=0)
  dT_dt=[(gamma_T*(Tstar(1)-T(1)) + u*(T(3)-T(1))/V(1) ...
	+k_v_highlat*(T(3)-T(1))/H(1)+k_h*(T(2)-T(1))/L(1) ) ...
	,(gamma_T*(Tstar(2)-T(2)) + u*(T(1)-T(2))/V(2) ...
	 +k_v_lowlat*(T(3)-T(2))/H(2)+k_h*(T(1)-T(2))/L(2) )  ...
	,(u*(T(2)-T(3))/V(3) ...
	 -(k_v_lowlat*(T(3)-T(2))+k_v_highlat*(T(3)-T(1)))/H(3) ) ...
       ];
else
  dT_dt=[(gamma_T*(Tstar(1)-T(1)) - u*(T(2)-T(1))/V(1) ...
	 +k_v_highlat*(T(3)-T(1))/H(1)+k_h*(T(2)-T(1))/L(1) ) ...
	 ,(gamma_T*(Tstar(2)-T(2)) - u*(T(3)-T(2))/V(2) ...
	 +k_v_lowlat*(T(3)-T(2))/H(2)+k_h*(T(1)-T(2))/L(2) )  ...
	 ,(-u*(T(1)-T(3))/V(3) ...
	 -(k_v_lowlat*(T(3)-T(2))+k_v_highlat*(T(3)-T(1)))/H(3) ) ...
	];
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dS_dt=time_rate_of_change_S ...
    (n,S,u,Sstar,gamma_S,FW,V,Area,H,L,k_h,k_v_highlat,k_v_lowlat)
% rate of change of salinity:

% beware of Matlab mess: make sure each of the elements of the vector
% dX_dt are sorrounded by ().  Otherwise continuation character at end
% of line makes matlab think that each line is a different vector
% element:
if (u>=0)
  dS_dt=[ ((FW(1) ...
	   +gamma_S*(Sstar(1)-S(1))*V(1) +u*(S(3)-S(1)) )/V(1) ...
	   +k_v_highlat*(S(3)-S(1))/H(1)+k_h*(S(2)-S(1))/L(1) ) ...
	  ,((FW(2) ...
	    +gamma_S*(Sstar(2)-S(2))*V(2) +u*(S(1)-S(2)) )/V(2) ...
	    +k_v_lowlat*(S(3)-S(2))/H(2)+k_h*(S(1)-S(2))/L(2) ) ...
	  ,(u*(S(2)-S(3))/V(3) ...
	 -(k_v_lowlat*(S(3)-S(2))+k_v_highlat*(S(3)-S(1)))/H(3) ) ...
	];
else
  dS_dt=[ ((FW(1)...
	   +gamma_S*(Sstar(1)-S(1))*V(1) -u*(S(2)-S(1)))/V(1) ...
	   +k_v_highlat*(S(3)-S(1))/H(1)+k_h*(S(2)-S(1))/L(1) ) ...
	  ,((FW(2)...
	     +gamma_S*(Sstar(2)-S(2))*V(2) -u*(S(3)-S(2)))/V(2) ...
	    +k_v_lowlat*(S(3)-S(2))/H(2)+k_h*(S(1)-S(2))/L(2) ) ...
	  ,(-u*(S(1)-S(3))/V(3) ...
	    -(k_v_lowlat*(S(3)-S(2))+k_v_highlat*(S(3)-S(1)))/H(3) ) ...
	];
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dALK_dt=time_rate_of_change_ALK ...
    (ALK,PO4,T,u,V,Area,H,L,k_h,k_v_highlat,k_v_lowlat...
     ,Redfield_PC,Redfield_PN,hPO4,light_factor,rPO4)
% rate of change of alkalinity:

% need to take care of units here, because biochemistry is in units of
% micro mol per liter, and export production in micro mol/ (cm^2 sec).
cm3_per_liter=1000.0;
S_ALK_surface_high=-Area(1)*(2*Rain(T(1))-1.0/Redfield_PN) ...
    *cm3_per_liter*export_production(PO4(1),hPO4,light_factor(1),rPO4(1));
S_ALK_surface_low=-Area(2)*(2*Rain(T(2))-1.0/Redfield_PN) ...
    *cm3_per_liter*export_production(PO4(2),hPO4,light_factor(2),rPO4(2));
S_ALK_deep=-(S_ALK_surface_high+S_ALK_surface_low);

% beware of Matlab mess: make sure each of the elements of the vector
% dX_dt are sorrounded by ().  Otherwise continuation character at end
% of line makes matlab think that each line is a different vector
% element:
if (u>=0)
  dALK_dt=[ ((S_ALK_surface_high ...
	      +u*(ALK(3)-ALK(1)) )/V(1) ...
	     +k_v_highlat*(ALK(3)-ALK(1))/H(1)+k_h*(ALK(2)-ALK(1))/L(1) ) ...
	    ,((S_ALK_surface_low ...
	       +u*(ALK(1)-ALK(2)) )/V(2) ...
	      +k_v_lowlat*(ALK(3)-ALK(2))/H(2)+k_h*(ALK(1)-ALK(2))/L(2) ) ...
	    ,((S_ALK_deep ...
	       +u*(ALK(2)-ALK(3)))/V(3) ...
	      -(k_v_lowlat*(ALK(3)-ALK(2))+k_v_highlat*(ALK(3)-ALK(1)))/H(3) ) ...
	  ];
else
  dALK_dt=[ ((S_ALK_surface_high ...
	      - u*(ALK(2)-ALK(1)))/V(1) ...
	     +k_v_highlat*(ALK(3)-ALK(1))/H(1)+k_h*(ALK(2)-ALK(1))/L(1) ) ...
	    ,((S_ALK_surface_low ...
	       - u*(ALK(3)-ALK(2)))/V(2) ...
	      +k_v_lowlat*(ALK(3)-ALK(2))/H(2)+k_h*(ALK(1)-ALK(2))/L(2) ) ...
	    ,((S_ALK_deep ...
	       -u*(ALK(1)-ALK(3)))/V(3) ...
	      -(k_v_lowlat*(ALK(3)-ALK(2))+k_v_highlat*(ALK(3)-ALK(1)))/H(3) ) ...
	  ];
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dTCO2_dt=time_rate_of_change_TCO2 ...
    (TCO2,PO4,T,u,V,Area,H,L,k_h,k_v_highlat,k_v_lowlat,PV,CO2_atm,CO2_g ...
     ,Redfield_PC,Redfield_PN,hPO4,light_factor,rPO4)
% rate of change of total co2:

% need to take care of units here, because biochemistry is in units of
% micro mol per liter, and export production in micro mol/ (cm^2 sec).
cm3_per_liter=1000.0;
S_TCO2_surface_high=-Area(1)*(1.0/Redfield_PC+Rain(T(1)))...
    *cm3_per_liter*export_production(PO4(1),hPO4,light_factor(1),rPO4(1));
S_TCO2_surface_low=-Area(2)*(1.0/Redfield_PC+Rain(T(2)))...
    *cm3_per_liter*export_production(PO4(2),hPO4,light_factor(2),rPO4(2));
S_TCO2_deep=-(S_TCO2_surface_high+S_TCO2_surface_low);

% beware of Matlab mess: make sure each of the elements of the vector
% dX_dt are sorrounded by ().  Otherwise continuation character at end
% of line makes matlab think that each line is a different vector
% element:
if (u>=0)
  dTCO2_dt=[ ((S_TCO2_surface_high + Area(1)*PV*(CO2_atm-CO2_g(1)) ...
	       +u*(TCO2(3)-TCO2(1)) )/V(1)...
	      +k_v_highlat*(TCO2(3)-TCO2(1))/H(1)+k_h*(TCO2(2)-TCO2(1))/L(1) ) ...
	     ,((S_TCO2_surface_low + Area(2)*PV*(CO2_atm-CO2_g(2)) ...
		+u*(TCO2(1)-TCO2(2)) )/V(2)...
	       +k_v_lowlat*(TCO2(3)-TCO2(2))/H(2)+k_h*(TCO2(1)-TCO2(2))/L(2) ) ...
	     ,((S_TCO2_deep ...
		+u*(TCO2(2)-TCO2(3)))/V(3) ...
	       -(k_v_lowlat*(TCO2(3)-TCO2(2))+k_v_highlat*(TCO2(3)-TCO2(1)))/H(3) ) ...
	  ];
else
  dTCO2_dt=[ ((S_TCO2_surface_high + Area(1)*PV*(CO2_atm-CO2_g(1)) ...
	       -u*(TCO2(2)-TCO2(1)))/V(1) ...
	      +k_v_highlat*(TCO2(3)-TCO2(1))/H(1)+k_h*(TCO2(2)-TCO2(1))/L(1) ) ...
	     ,((S_TCO2_surface_low + Area(2)*PV*(CO2_atm-CO2_g(2)) ...
		-u*(TCO2(3)-TCO2(2)))/V(2) ...
	       +k_v_lowlat*(TCO2(3)-TCO2(2))/H(2)+k_h*(TCO2(1)-TCO2(2))/L(2) ) ...
	     ,((S_TCO2_deep ...
		-u*(TCO2(1)-TCO2(3)))/V(3) ...
	       -(k_v_lowlat*(TCO2(3)-TCO2(2))+k_v_highlat*(TCO2(3)-TCO2(1)))/H(3) ) ...
	   ];
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dPO4_dt=time_rate_of_change_PO4 ...
    (PO4,u,V,Area,H,L,k_h,k_v_highlat,k_v_lowlat ...
     ,Redfield_PC,Redfield_PN,hPO4,light_factor,rPO4)
% rate of change of nutrient:

% need to take care of units here, because biochemistry is in units of
% micro mol per liter, and export production in micro mol/ (cm^2 sec).
cm3_per_liter=1000.0;
S_PO4_surface_high=-Area(1)*cm3_per_liter*export_production(PO4(1),hPO4,light_factor(1),rPO4(1));
S_PO4_surface_low=-Area(2)*cm3_per_liter*export_production(PO4(2),hPO4,light_factor(2),rPO4(2));
S_PO4_deep=-(S_PO4_surface_high+S_PO4_surface_low);

% beware of Matlab mess: make sure each of the elements of the vector
% dX_dt are sorrounded by ().  Otherwise continuation character at end
% of line makes matlab think that each line is a different vector
% element:
if (u>=0)
  dPO4_dt=[ ((S_PO4_surface_high ...
	     +u*(PO4(3)-PO4(1)) )/V(1)...
	     +k_v_highlat*(PO4(3)-PO4(1))/H(1)+k_h*(PO4(2)-PO4(1))/L(1) ) ...
	    ,((S_PO4_surface_low ...
	       +u*(PO4(1)-PO4(2)) )/V(2)...
	      +k_v_lowlat*(PO4(3)-PO4(2))/H(2)+k_h*(PO4(1)-PO4(2))/L(2) ) ...
	    ,((S_PO4_deep ...
	       +u*(PO4(2)-PO4(3)))/V(3) ...
	      -(k_v_lowlat*(PO4(3)-PO4(2))+k_v_highlat*(PO4(3)-PO4(1)))/H(3) ) ...
	  ];
else
  dPO4_dt=[ ((S_PO4_surface_high ...
	      -u*(PO4(2)-PO4(1)))/V(1) ...
	     +k_v_highlat*(PO4(3)-PO4(1))/H(1)+k_h*(PO4(2)-PO4(1))/L(1) ) ...
	    ,((S_PO4_surface_low ...
	       -u*(PO4(3)-PO4(2)))/V(2) ...
	      +k_v_lowlat*(PO4(3)-PO4(2))/H(2)+k_h*(PO4(1)-PO4(2))/L(2) ) ...
	    ,((S_PO4_deep ...
	       -u*(PO4(1)-PO4(3)))/V(3) ...
	      -(k_v_lowlat*(PO4(3)-PO4(2))+k_v_highlat*(PO4(3)-PO4(1)))/H(3) ) ...
	  ];
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [CO2_g,H2CO3,OHm,Hp,pH,HCO3m,CO3m2,B_OH_3,B_OH_4m,B_T]= ...
      solve_carbonate_system(ALK,TCO2,T,S,Hp_init_guess)

% Inputs:
% -------
% alkalinity, total co_2, temperature, salinity, pressure.

% Outputs:
% --------
% CO2_g    = co_2(g)
% H2CO3    = H_2co_3^*  = H_2co_3+co_2(aq)
% OHm      = OH^-
% Hp       = H^+
% pH       = -log10(Hp)
% HCO3m    = HCO_3^-
% CO3m2    = CO_3^=
% B_OH_3   = B(OH)_3
% B_OH_4m  = B(OH)_4^-

% solve for all three boxes:
for i=1:3
  % calculate dissolution coefficients to save operations during
  % the iterations below:
  K0(i)=K_0(T(i),S(i));
  K1(i)=K_1(T(i),S(i));
  KB(i)=K_B(T(i),S(i));
  K2(i)=K_2(T(i),S(i));
  Kw(i)=K_w(T(i),S(i));
  % calculate total Borate from salinity; multiply by 1.e6 to get
  % it in micromole:
  rho_0       = 1.024;
  B_T(i)      = 1.e6*1.212*10^-5*S(i)*rho_0;
  % initialize a=1/[H^+] and ALK_C from previous time step:
  a=1/Hp_init_guess(i);
  max_carbonate_iterates=20;
  min_carbonate_iterates=5;
  iterate=1;
  epsilon_a=1.0;
  while iterate<min_carbonate_iterates || ( iterate< max_carbonate_iterates && epsilon_a > 1.0e-6 )
    iterate=iterate+1;
    %% Calculate carbonate alkalinity:
    ALK_C=ALK(i)-( (a*KB(i)*B_T(i))/(1+a*KB(i)) + a*Kw(i) - 1/a );
    %% gamma:
    gamma=ALK_C/TCO2(i);
    %% solve for a=1/[H^+]:
    a_new=(-(1-gamma)*K1(i)...
       +sqrt((1-gamma)^2*K1(i)^2 ...
	     + 4*(2-gamma)*K1(i)*K2(i)*gamma)...
       )/(2*(2-gamma)*K1(i)*K2(i));
    %% evaluate relative increment to a in this iteration:
    epsilon_a=abs(a_new-a)/a;
    a=a_new;
  end
  if iterate == max_carbonate_iterates
    fprintf(1,['reached max carbonate system iterates,' ...
               ' t=%g, iterate=%d, Alk=%8.2f, DIC=%8.2f\n'],t,iterate,ALK,TCO2);
  end
  
  
%% Now solve for all other carbonate system variables:
Hp(i)       = 1/a;
% add 6 to pH because units of Hp are micro mole/liter:
pH(i)       = -log10(Hp(i))+6;
OHm(i)      = Kw(i)*a;
H2CO3(i)    = (1/(a*K1(i)+a^2*K1(i)*K2(i)+1))*TCO2(i);
HCO3m(i)    = ((a*K1(i))/(a*K1(i)+a^2*K1(i)*K2(i)+1))*TCO2(i);
CO3m2(i)    = ((a^2*K1(i)*K2(i))...
               /(a*K1(i)+a^2*K1(i)*K2(i)+1))*TCO2(i);
CO2_g(i)    = H2CO3(i)/K0(i);
B_OH_4m(i)  = ((a*KB(i))/(1+a*KB(i)))*B_T(i);
B_OH_3(i)   = B_OH_4m(i)/(a*KB(i));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [L,H,W,V,Area,Sv,year,u0,alpha,beta,gamma_T,gamma_S ...
	  ,k_h,k_v_highlat,k_v_lowlat,PV ...
	  ,FW,Tstar,Sstar,delta_t,N_max,ubar,T,S,ALK,TCO2,PO4,Hp ...
	  ,Redfield_PC,Redfield_PN,hPO4,light_factor,rPO4 ...
	  ,solve_physics,solve_geochemistry] ...
    =set_parameters

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

solve_physics=0;
solve_geochemistry=1;
% box geometry:
% north-south size:
L(1)=5000.0e5;
L(2)=1000.0e5;
L(3)=L(1);
% depth:
H(1)=1.0e5;
H(2)=H(1);
H(3)=4.0e5;
if (H(1) ~= H(2))
  disp(['need to check that advection/ diffusion equation are' ...
	' consistent with the case of h1!=h2!!']) 
  quit
end
% width:
W=4000.0e5;
% volumes:
V(1)=L(1)*H(1)*W;
V(2)=L(2)*H(2)*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;
% 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);
gamma_S=0.0/((days_per_50m_mixed_layer*86400.0)*(H(1)/50.e2));
FW(1)=-1.1*(100.0/year)*35.0*Area(1);
FW(2)=-FW(1);
% the temperature to which the model surface temperature is restored:
Tstar=[0.0,22.0]; 
% the salinity to which the model surface salinity is restored (only
% if gamma_S is not zero):
Sstar=[34.0,36.0];

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

% Initial conditions:
ubar= 20*Sv;
T    = [1.0,20.0 ,2.0];
S    = [34.0,35.5,34.5];
PO4  = [2.09,2.09,2.09];
ALK  = [2371,2371,2371];
TCO2 = [2258,2258,2258];
Hp=[10^-8.2,10^-8.2,10^-8.2];

% some biochemistry parameters:
Redfield_PC=1.0/122;
Redfield_PN=1.0/16;
% paraemters for export production parameterization:
hPO4=1.0;
light_factor=[350.0,200.0];
rPO4=1.0e-10*[0.25,250.0];
% piston velocity for air-sea co2 exchanges
PV=1.0*0.59e-2 ; % \approx 5 m/day


% mixing coefficients:
k_h        =  1.0e7/(0.5*(L(1)+L(2)));
k_v_highlat= 200.0e0/(0.5*(H(1)+H(3)));
k_v_lowlat =  0.0e0/(0.5*(H(1)+H(3)));

% Test the density calculation:
if 0
  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
end

% Test the geochemical coefficient calculation and the solution of the
% carbonate system:
if 0
  disp('Testing geochemical coefficient calculation:')
  T1=18,S1=35
  Kw=K_w(T1,S1),k0=K_0(T1,S1),k1=K_1(T1,S1),k2=K_2(T1,S1),kb=K_B(T1,S1)
  disp('Testing solution of carbonate system:')
  Hp_init_guess=Hp;
  T,S
  [CO2_g,H2CO3,OHm,Hp,pH,HCO3m,CO3m2,B_OH_3,B_OH_4m,B_T]= ...
      solve_carbonate_system(ALK,TCO2,T,S,Hp_init_guess)
  disp('Testing export production calculation (printed in units of gr C/yr/m^2):')
  EP1=export_production(PO4(1),hPO4,light_factor(1),rPO4(1));
  EP2=export_production(PO4(2),hPO4,light_factor(2),rPO4(2));
  % our units are PO4/cm^2/sec.  Need to convert to gr C/year and to m^2:
  export_production_gr_C_per_m2=...
      [EP1*(122+Rain(Tnow(1))),EP2*(122+Rain(Tnow(2)))] ...
      *(12.011*10^-6)*(365*86400)*100
  export_production_Gt_C_per_yr=...
      [EP1*(122+Rain(Tnow(1)))*Area(1),EP2*(122+Rain(Tnow(2)))*Area(2)] ...
      *(12.011*10^-6)*(365*86400)*1.e-15
  disp('Testing rain ratio paramterization:')
  T,rain_T1=Rain(T(1)),rain_T2=Rain(T(2)),rain_T3=Rain(T(3))
end

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


%     These functions calculate the dissociation coefficients needed for
%     calculation of pCO2, and the solubility of oxygen.

function output=K_0(T,S)
%     solubility of CO2 in seawater
%         [CO2*]
%     K_0=-------- 
%         [CO2g]
%     from Weiss (1974)
% Eli: looks like units are mole CO2/(kg Sea Water*atm)
% which is the same as the units used here: micro mole CO2/(kg SW*ppm)
T=T+273.15;
output=exp( ...
    (9345.17/T)       -60.2409+23.3585*log(T/100.) ...
    +S*( 0.023517-0.00023656*T  +  ...
         0.00000047036*T*T ));


function output=K_1(T,S)
%        [H+][HCO3-]  
%     K_1=-----------      
%          [CO2*]
%     from Millero (1987)

T=T+273.15;
output=10.^(-((6320.18/T)-126.3405+19.568*log(T)+ ...
	      (19.894-(840.39/T)-3.0189*log(T))* ...
	      sqrt(S)+0.00678*S) );
% convert units to micromole:
output=output*1.e6;


function output=K_2(T,S)
%               [H+][CO3--]                         
%        and K_2=-----------  from Millero (1987)
%                 [HCO3-]

T=T+273.15;
output=10.^(-((5143.69/T)-90.1833+14.613*log(T)+ ...
	      (17.176-(690.59/T)-2.6719*log(T))* ...
	      sqrt(S)+0.02169*S) );
% convert units to micromole:
output=output*1.e6;


function output=K_B(T,S)
%     for Boric acid,        [H+][B(OH4)-]
%                    K_B=-------------- is from Johansson & Wedborg (1982)
%                          [B(OH)3 

T=T+273.15;
output=10.^(-((1030.5/T)+5.5076-0.015469*S+ ...
	      1.5339d-4*S*S ) );
% convert units to micromole:
output=output*1.e6;


function output=K_w(T,S)
%     for water. K_w=[h+]*[OH-], from Dickson & Riley (1979)
%     pKw=3441.0*(1./T)+2.241-0.09415*sqrt(S)

T=T+273.15;
output=10.^(-((3441.0/T)+2.241-0.09415*sqrt(S) ));
% convert units to micromole:
output=output*1.e12;
