%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Alkalinity_tutorial
%% demonstrate what happens to difference carbonate species when
%% alkalinity is changed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%close all;

%profile on
p=set_parameters;

response_to_Alkalinity_change(p); return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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(x,[y1,y2])
%% If y1,y2 are column vectors, use
%% set_axis_limits(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(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 Pos=my_subplot(nrows,ncols,ifig)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Replace Matlab's subplot, to allow a better control of spacing
% between panels.  Eli, June 12, 2007.

if ifig>nrows*ncols;
  fprintf('** problem in my_subplot: ifig=%d > nrows*ncols=%d*%d=%d.\n'...
          ,ifig,nrows,ncols,nrows*ncols);
  disp('pausing.');
  pause
end

%fprintf(1,'ir,jc,nx_panels,ny_panels=%d, %d, %d, %d\n',ir,jc,nx_panels,ny_panels)
% left and right margins of page:
left_margin=0.12; 
right_margin=0;
% top margin of page:
y_top_margin=0.05;
% bottom margin of page:
y_bot_margin=0.1;
% left and right margins of each panel (relative to panel width):
dx=0.2;
% top and bottom margins of each panel (relative to panel height):
dy=0.22;
% width and height of entire plot on page:
W=1-(left_margin+right_margin); 
H=1-(y_top_margin+y_bot_margin);
% separation in x and y between panels:
xsep=W/ncols; ysep=H/nrows;
% width and height of each pannel:
width=xsep*(1-dx); height=ysep*(1-dy);
% indices of row and column of the panel being drawn:
ir=floor((ifig-1)/ncols)+1; jc=mod(ifig,ncols); if jc==0; jc=ncols; end
% position:
left=left_margin+xsep*(jc-1); bottom=(1-y_top_margin)-ysep*ir;
Pos=[left bottom width height];
axes('Position',Pos);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function a=temperature(t,p);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Celsius!
a=p.Temperature_C_0;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The following functions calculate the dissociation coefficients
%% needed for calculation of pCO2, and the solubility of oxygen,
%% taken from paleo_co2_tutorial.m, with slight naming changes.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function K0=K_0(p,P,T,S)
%     solubility of CO2 in seawater.  
%     by definition, CO2*=H2CO3*=CO2(aq)+H2CO3.
%         [CO2*]
%     K_0=-------- 
%         [CO2g]
%     from Zeebe's book
% 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;
K0=exp( ...
    (9345.17/T)-60.2409+23.3585*log(T/100.) ...
    +S*( 0.023517-0.00023656*T  +  ...
         0.00000047036*T*T ));


function K1=K_1(p,P,T,S)
%        [H+][HCO3-]  
%     K_1=-----------      
%          [CO2(aq)]
%     from Millero (1987)

%% Old:
%% T=T+273.15;
%% K1=10.^(-((6320.18/T)-126.3405+19.568*log(T)+ ...
%% 	      (19.894-(840.39/T)-3.0189*log(T))* ...
%% 	      sqrt(S)+0.00678*S) );

%% From Zeebe's book, Eli:
T=T+273.15;
tmp1 = 2.83655 - 2307.1266 ./ T - 1.5529413 .* log(T);
tmp2 =         - (0.20760841 + 4.0484 ./ T) .* sqrt(S);
tmp3 =         + 0.08468345 .* S - 0.00654208 .* S .* sqrt(S);   
tmp4 =         + log(1 - 0.001005 .* S);
lnK1roy = tmp1 + tmp2 + tmp3 + tmp4;
K1 = exp(lnK1roy);

%% pressure effect:
K1 = K1*exp(p.lnkpok0(1));

% convert units to micromole:
K1=K1*1.e6;


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

T=T+273.15;
%% Old:
%K2=10.^(-((5143.69/T)-90.1833+14.613*log(T)+ ...
%	      (17.176-(690.59/T)-2.6719*log(T))* ...
%	      sqrt(S)+0.02169*S) );

%% Zeebe's book, Eli:
tmp1 = -9.226508 - 3351.6106 ./ T - 0.2005743 .* log(T);
tmp2 = (-0.106901773 - 23.9722 ./ T) .* sqrt(S);
tmp3 = 0.1130822 .* S - 0.00846934 .* S.^1.5 + log(1 - 0.001005 * S);
lnK2roy = tmp1 + tmp2 + tmp3;
K2 = exp(lnK2roy);

%% pressure effect:
K2 = K2*exp(p.lnkpok0(2));

% convert units to micromole:
K2=K2*1.e6;


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

T=T+273.15;
%% Old:
%output=10.^(-((1030.5/T)+5.5076-0.015469*S+ ...
%	      1.5339d-4*S*S ) );

%% Zeebe's book, Eli:
tmp1 =  (-8966.90-2890.53*sqrt(S)-77.942*S+1.728*S.^(3./2.)-0.0996*S.*S);
tmp2 =   +148.0248+137.1942*sqrt(S)+1.62142*S;
tmp3 = +(-24.4344-25.085*sqrt(S)-0.2474*S).*log(T);
lnKb = tmp1 ./ T + tmp2 + tmp3 + 0.053105*sqrt(S).*T;
Kb = exp(lnKb);

%% pressure effect:
Kb = Kb*exp(p.lnkpok0(3));

% convert units to micromole:
Kb=Kb*1.e6;


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

%% Old:
%%T=T+273.15;
%%output=10.^(-((3441.0/T)+2.241-0.09415*sqrt(S) ));

%% Zeebe's book, Eli:
T=T+273.15;
Kw=exp(148.96502-13847.26/T-23.6521*log(T) ...
           +(118.67/T-5.977+1.0495*log(T))*sqrt(S)-0.01615*S);

%% pressure effect:
Kw = Kw*exp(p.lnkpok0(4));

% convert units to micromole:
Kw=Kw*1.e12;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function K=K_sp(p,P,T,S)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if p.do_use_Ksp_aragonite==1;
    K=K_sp_aragonite(p,P,T,S);
  else
    K=K_sp_calcite(p,P,T,S);
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function output=K_sp_calcite(p,P,T,S)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The Calcite solubility product according to Mucci (1983) on the molal scale.
%% from Jack Silverman.
%P=p.K_sp_pressure; % atm
%dVcal=-(47.5-0.23.*T);
%Kcal=10.^(-171.9065-0.077993.*(T+273.15)+2839.319./(T+273.15)+...
%   71.595.*log10(T+273.15)+(-0.77712+0.0028426.*(T+273.15)+178.34./(T+...
%   273.15)).*S.^0.5-0.07711.*S+0.0041249.*S.^1.5);
%Kcalp=Kcal.*10.^-(dVcal.*(P-1)./(2.303*82.057834.*(273.15+T)));

%% Zeebe's book, Eli:
T=T+273.15;
tmp1 = -171.9065-0.077993.*T+2839.319./T+71.595.*log10(T);
tmp2 = +(-0.77712+0.0028426.*T+178.34./T).*sqrt(S);
tmp3 = -0.07711.*S+0.0041249.*S.^1.5;
log10Kspc = tmp1 + tmp2 + tmp3;
Kspc = 10.^(log10Kspc);

%% pressure effect:
Kspc = Kspc*exp(p.lnkpok0(7));

%% Eli: convert to micromole^2:
output=Kspc*1.e12;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Kspa=K_sp_aragonite(p,P,T,S)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The agonite solubility product according to Mucci (1983) on the molal scale.
%% from Jack Silverman.
P=p.K_sp_pressure; % atm

%% Old:
% dVarg=-(45-0.23.*T);
% Kspa=10.^(-171.945-0.077993.*(T+273.15)+2903.293./(T+273.15)+...
%    71.595.*log10(T+273.15)+(-0.068393+0.0017276.*(T+273.15)+...
%    88.135./(T+273.15)).*S.^0.5-0.10018.*S+0.0059415.*S.^1.5);
% Kspa=Kspa.*10.^-(dVarg.*(P-1)./(2.303*82.057834.*(273.15+T)));

%% Zeebe's book, Eli:
T=T+273.15;
tmp1 = -171.945-0.077993.*T+2903.293./T+71.595.*log10(T);
tmp2 = +(-0.068393+0.0017276.*T+88.135./T).*sqrt(S);
tmp3 = -0.10018.*S+0.0059415.*S.^1.5;
log10Kspa = tmp1 + tmp2 + tmp3;
Kspa = 10.^(log10Kspa);

%% pressure effect:
Kspa = Kspa*exp(p.lnkpok0(8));

% Eli: convert to micromole^2:
Kspa=Kspa*1.e12;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function output=K_HSO4(p,P,T,S)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% P in bar, T celsius, S ppt
%            [H+][SO4--]  
%     K_HSO4=-----------      
%              [HSO4-]
P=p.K_sp_pressure; % atm

if 1
  %% from aquaEnv: http://r-forge.r-project.org/projects/aquaenv/
  %% AquaEnv - Aquatic modelling Environment in R, Andreas by
  %% F. Hofmann, Karline Soetaert, Filip J. R. Meysman, 2009
  p=P;
  T         = T + 273.15;
  %% ionic strength in mol/kg-solution (molinity)
  %% salinity S in practical salinity units:
  I=(19.924*S/(1000-1.005*S));
  Iterms=[I,I^2,sqrt(I),I^1.5];
  A =  324.57*Iterms(3) - 771.54*Iterms(1) + 141.328;
  B =   35474*Iterms(1) +   1776*Iterms(2) - 13856*Iterms(3) - 2698*Iterms(4) - 4276.1;
  C = 114.723*Iterms(1) - 47.986*Iterms(3) - 23.093;
  D = 0;
  E = 0;
  lnK = A + (B/T) + C*log(T) + D*T + E*T^2;
  t         = T;
  t2        = t^2;
  coeff = [-18.03, 0.0466, 0.3160e-3,- 4.53, 0.0900,0];
  deltaV    = coeff(1) + coeff(2)*t + coeff(3)*t2;
  deltaK    = (coeff(4) + coeff(5)*t + coeff(6)*t2)/1000; % the 1000 is correction from Lewis1998
  R=83.14472;        % (bar*cm3)/(mol*K) the gas constant 
  deltaPlnK = -(deltaV/(R*T))*p + (0.5*(deltaK/(R*T)))*(p^2); 
  %% salinity S in practical salinity units the conversion factor from
  %% molality (mol/kg-H2O) to molinity (mol/kg-solution)
  molal2molin=(1-0.001005*S);
  K_HSO4 = exp(lnK + deltaPlnK + log(molal2molin));
  % conversion factor to pH_{NBS} scale (Pilson table F.2, p 405):
  f_H=0.739+0.0307*S+0.0000794*S^2+0.00006443*T-0.000117*S*T;
  output=K_HSO4*f_H;
  %% convert units to micromole:
  output=output*1.e6;
end

if 0
  %% Khoo 1977 from UNESCO report:
  T=T+273.15;
  lnT=log(T);
  lnBeta=1121.1/T-45.2167+7.484*lnT+0.0011984*S-1.2613e-4*S^2;
  beta=exp(lnBeta);
  output=1/beta;
end

if 0
  %% From "An introduction to the chemistry of the sea" by Michael
  %% E. Q. Pilson; based on Dickson 1990: j.chem.thermodynamics, vol 22,
  %% pp 113-127. "Standard potential of the reaction : AgCl(s) + H2(g)
  %% = Ag(s) + HCl(aq), and the standard acidity constant of the ion
  %% HSO4- in synthetic sea water from 273.15 to 318.15 K".
  T=T+273.15;
  lnT=log(T);
  I=19.92*S/(1000-1.00488*S);
  
  lnK_T=-4276.1/T+141.328-23.093*lnT+(-13856/T+324.57+47.986*lnT)*sqrt(I) ...
        +(35474/T-771.54+114.723*lnT)*I+(-2698/T)*I^(1.5)+(1776/T)*I^2;
  
  lnK_T
  output=exp(lnK_T);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%function output=K0_H2S(p,P,T,S)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         [H2S(aq)]
%     K0=-------- 
%         [H2S(g)]
% Henry constant for H2S
% documentation: Doc/K0_H2S.tex
% output = 0.100908; % micromole/(kg*ppm atm);
  
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function output=K1_H2S(p,P,T,S)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% first dissociation constant for H2S
%% From AquaEnv, compared for P=0 to Table F.6 p 408, Pilson's book.
P=p.K_sp_pressure; % atm
T = T + 273.15;
sqrtS  = sqrt(S);
Sterms = [S^2, sqrtS, S*sqrtS];
A = 225.838  + 0.3449*Sterms(2) - 0.0274*S;
B = -13275.3;
C = -34.6435;
D = 0;
E = 0;
lnK = A + (B/T) + C*log(T) + D*T + E*T^2;
R=83.14472;        % (bar*cm3)/(mol*K) the gas constant 
t = T;
t2        = t^2;
coeff=[-14.80, 0.0020, -0.4000e-3,  2.89, 0.0540, 0];
deltaV    = coeff(1) + coeff(2)*t + coeff(3)*t2;
deltaK    = (coeff(4) + coeff(5)*t + coeff(6)*t2)/1000; %the 1000 is correction from Lewis1998
deltaPlnK = -(deltaV/(R*T))*P + (0.5*(deltaK/(R*T)))*(P^2); 

K_H2S = exp(lnK + deltaPlnK);
%% + log(scaleconvert(S, t, P=0, SumH2SO4, SumHF, khf=khf, khso4=khso4)$tot2sws) +
%% + log(scaleconvert(S, t, P, SumH2SO4, SumHF, khf=khf, khso4=khso4)$sws2free))

output=K_H2S;
%% conversion factor to pH_{NBS} scale (Pilson table F.2, p 405):
%f_H=0.739+0.0307*S+0.0000794*S^2+0.00006443*T-0.000117*S*T;
%output=output*f_H;

%% convert units to micromole:
output=output*1.e6;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [pH,Hp,OHm,CO2g,H2CO3,HCO3m,CO3m2,BOH3,BOH4m,HSO4m,SO4m2,H2Saq,HSm] ...
  =solve_carbonate_system_with_sulfur(t,p,Alk,C_T,B_T,SO4_T,H2S_T_ell,T,S,init_guess)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

p_fsolve.C_T=C_T;
p_fsolve.B_T=B_T;
p_fsolve.SO4_T=SO4_T;
p_fsolve.H2S_T_ell=H2S_T_ell;
p_fsolve.Alk=Alk;
p_fsolve.T=T;
p_fsolve.S=S;

% misc reaction rate constants that need to be calculated from T,S,P
% only once per time step, so are better done here than in
% carbonate_sulfate_eqns:
P=p.K_sp_pressure;
T=temperature(t,p);
p=pressure_effect_on_K(t,p);
S=p.Salinity0;
K.K_w=K_w(p,P,T,S);
K.K_0=K_0(p,P,T,S);
K.K_1=K_1(p,P,T,S);
K.K_2=K_2(p,P,T,S);
K.K_HSO4=K_HSO4(p,P,T,S);
K.K_B=K_B(p,P,T,S);
K.K1_H2S=K1_H2S(p,P,T,S);

%% nondimensionalize initial guess:
init_guess=(init_guess-p.carb_sulf_offset)./p.carb_sulf_scale;
[x,fval,exitflag,output,J]=fsolve(@(X)carbonate_sulfate_eqns(X,t,p,p_fsolve,K),init_guess,p.optim_options);
%% dimensionalize the fsolve solution:
X=x.*p.carb_sulf_scale+p.carb_sulf_offset;

%fprintf(1,'iterations=%d\n',output.iterations);

if 0 %% for debugging:
  format long
  disp('in solve_carbonate_system_with_sulfur:, nondim x=');
  disp(x);
  disp('dimensional X=');
  disp(X);
end

if exitflag<1
  output
  output.message
  exitflag
  fval
  fprintf(1,'t=%gyrs\n',t);
  fprintf(1,'F='); fprintf(1,'%g; ',fval); fprintf(1,'\n');
  fprintf(1,'X='); fprintf(1,'%g; ',X); fprintf(1,'\n');
  t,p,Alk,C_T,B_T,SO4_T,H2S_T_ell,T,S,init_guess
  disp('*** fsolve problem, pausing in solve_carbonate_system_with_sulfur.');
  pause
end

% assign output arguments:
Hp      = X(1);
OHm     = X(2);
CO2g    = X(3);
H2CO3   = X(4);
HCO3m   = X(5);
CO3m2   = X(6);
BOH3    = X(7);
BOH4m   = X(8);
HSO4m   = X(9);
SO4m2   = X(10);
H2Saq   = X(11);
HSm     = X(12);
pH=-log10(Hp)+6;

if 0
  fprintf(1,'solution:\n')
  for i=1:length(X)
    lab=p.carb_sulf_label(i,:);
    fprintf(1,'%s=%g,\n',lab,X(i));
  end
  fprintf(1,'pH   =%g\n',-log10(X(1))+6);
  fprintf(1,'\n')
end
%fprintf(1,'J at end of optimization:\n'); printJ(J);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [F,J]=carbonate_sulfate_eqns(X,t,p,p_fsolve,K)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Called by fsolve which is invoked in solve_carbonate_sulfate_system;
% returns the equations of the carbonate+sulfate system and the
% Jacobian.

  %% get variables from unknown vector X:  X is nondimensional, so
  %% for example X(1)=(Hp-p.carb_sulf_offset(1))/p.carb_sulf_scale(1):
  Hp      = X(1) *p.carb_sulf_scale(1) +p.carb_sulf_offset(1);
  OHm     = X(2) *p.carb_sulf_scale(2) +p.carb_sulf_offset(2);
  CO2g    = X(3) *p.carb_sulf_scale(3) +p.carb_sulf_offset(3);
  H2CO3   = X(4) *p.carb_sulf_scale(4) +p.carb_sulf_offset(4);
  HCO3m   = X(5) *p.carb_sulf_scale(5) +p.carb_sulf_offset(5);
  CO3m2   = X(6) *p.carb_sulf_scale(6) +p.carb_sulf_offset(6);
  BOH3    = X(7) *p.carb_sulf_scale(7) +p.carb_sulf_offset(7);
  BOH4m   = X(8) *p.carb_sulf_scale(8) +p.carb_sulf_offset(8);
  HSO4m   = X(9) *p.carb_sulf_scale(9) +p.carb_sulf_offset(9);
  SO4m2   = X(10)*p.carb_sulf_scale(10)+p.carb_sulf_offset(10);
  H2Saq   = X(11)*p.carb_sulf_scale(11)+p.carb_sulf_offset(11);
  HSm     = X(12)*p.carb_sulf_scale(12)+p.carb_sulf_offset(12);

  %% get some additional needed quantities:
  C_T   = p_fsolve.C_T;
  B_T   = p_fsolve.B_T;
  SO4_T   = p_fsolve.SO4_T;
  H2S_T_ell = p_fsolve.H2S_T_ell;
  Alk   = p_fsolve.Alk;

  %% Calculate equations:
  F=[ ...
    K.K_w-Hp*OHm,
    K.K_0*CO2g-H2CO3,
    K.K_1*H2CO3-Hp*HCO3m,
    K.K_2*HCO3m-Hp*CO3m2,
    K.K_HSO4*HSO4m-Hp*SO4m2,
    K.K_B*BOH3-Hp*BOH4m,
    K.K1_H2S*H2Saq-Hp*HSm,
    H2CO3+HCO3m+CO3m2-C_T,
    BOH4m+BOH3-B_T,
    HSO4m+SO4m2-SO4_T,
    HSm+H2Saq-H2S_T_ell,
    HCO3m+2*CO3m2+BOH4m+HSO4m+2*SO4m2+HSm+OHm-Hp-Alk ...
    ]';
  
  %% Calculate Jacobian:
  
  J=[ -OHm,	-Hp,	0,	0,	0,	0,	0,	0,	0,		0,	0,		0 ;
      0,	0,	K.K_0,	-1,	0,	0,	0,	0,	0,		0,	0,		0 ;
      -HCO3m,	0,	0,	K.K_1,	-Hp,	0,	0,	0,	0,		0,	0,		0 ;
      -CO3m2,	0,	0,	0,	K.K_2,	-Hp,	0,	0,	0,		0,	0,		0 ;
      -SO4m2,	0,	0,	0,	0,	0,	0,	0,	K.K_HSO4,	-Hp,	0,		0 ;
      -BOH4m,	0,	0,	0,	0,	0,	K.K_B,	-Hp,	0,		0,	0,		0 ;
      -HSm,	0,	0,	0,	0,	0,	0,	0,	0,		0,	K.K1_H2S,	-Hp ;
      0,	0,	0,	1,	1,	1,	0,	0,	0,		0,	0,		0 ;
      0,	0,	0,	0,	0,	0,	1,	1,	0,		0,	0,		0 ;
      0,	0,	0,	0,	0,	0,	0,	0,	1,		1,	0,		0 ;
      0,	0,	0,	0,	0,	0,	0,	0,	0,		0,	1,		1 ;
      -1,	1,	0,	0,	1,	2,	0,	1,	1,		2,	0,		1 ];

  if 0
    %% Add inequality penalties in case the solver gets to negative
    %% carbonate variables:
    %% dimensional variables:
    Xdim=X.*p.carb_sulf_scale+p.carb_sulf_offset;
    %% find indices of negative variables:
    I_penalty=find(Xdim<0);
    %% add penalty:
    F(I_penalty)=F(I_penalty)+Xdim(I_penalty).^2;
    %% add to Jacobian:
    for i=I_penalty
      J(i,i)=J(i,i)+2*Xdim(i);
    end
  end
  
  %% scale the Jacobian:
  N=length(F);
  for i=1:N
    J(i,:)=J(i,:).*p.carb_sulf_scale(:)';
  end

%%fprintf(1,'analytical J:\n'); printJ(J); pause

%fprintf(1,'F='); fprintf(1,'%g; ',F); fprintf(1,'\n');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function printJ(J)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N=length(J(:,1));
for i=1:N
  fprintf(1,'%g, ',J(i,:)); fprintf(1,'\n');
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function response_to_Alkalinity_change(p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% reproduce figure 1.1.3, p 7, from Zeebe and Wolf-Gladrow
%% or firgure 4-5 from Stum & Morgan, p 133
  t=0;
  T=temperature(t,p); % cannot change temp here, must change the
                      % temperature function.
  S=35;
  init_guess=p.carb_sulf_init_guess;
  P=500;
  p.K_sp_pressure=P;
  p=pressure_effect_on_K(t,p);

  %% ------------------------------------------------------------------------
  %% first, print carbonate solution for a specific choice of the
  %% input parameters, to compare with Zeebe's csys3 matlab program:
  %% ------------------------------------------------------------------------
  C_T=2000;
  Alk=2150+p.SO40*2;
  B_T=416;
  SO4_T=p.SO40;
  H2S_T_ell=0;
  disp('carbonate system solution:');
  [pH,Hp,OHm,CO2g,H2CO3,HCO3m,CO3m2,BOH3,BOH4m,HSO4m,SO4m2,H2Saq,HSm] ...
    =solve_carbonate_system_with_sulfur(t,p,Alk,C_T,B_T,SO4_T,H2S_T_ell,T,S,init_guess);
%    =solve_carbonate_system_carbonate_only_iterative(t,p,Alk,C_T,B_T,SO4_T,H2S_T_ell,T,S,init_guess);
  fprintf(1,[' t=%gMyr,\n' ...
             ' Alk=%g,\n C_T=%g,\n B_T=%g,\n SO4_T=%g,\n H2S_T_ell=%g,\n ' ...
             'CO2g=%g,\n H2CO3=%g,\n OHm=%g,\n Hp=%g,\n pH=%g,\n ' ...
             'HCO3m=%g,\n CO3m2=%g,\n BOH3=%g,\n BOH4m=%g\n'] ...
          ,t/1.e6,Alk,C_T,B_T,SO4_T,H2S_T_ell,CO2g,H2CO3,OHm,Hp,pH, ...
          HCO3m,CO3m2,BOH3,BOH4m);  

  CO3m2_sat=K_sp(p,P,T,S)/p.Ca2p0;
  fprintf(' CO3m2(sat)=%11.8g\n',CO3m2_sat);


  %% set input values for contour plot:  
  B_T=p.B_T0;
  SO4_T=p.SO40;
  H2S_T_ell=p.H2S0;
  
  %% ------------------------------------------------------------------------
  %% Contour co2 and pH as function of Alk and DIC:
  %% ------------------------------------------------------------------------
  %% Zeebe:
  DIC_axis=1950:20:2150;
  Alk_axis=2250:20:2450;
  %% Stum & Morgan:
  %DIC_axis=0:100:2500;
  %Alk_axis=-500:100:3500;

  N_DIC=length(DIC_axis);
  N_Alk=length(Alk_axis);
  CO2g_contour=zeros(N_DIC,N_Alk);
  H2CO3_contour=zeros(N_DIC,N_Alk);
  pH_contour=zeros(N_DIC,N_Alk);
  fprintf(1,'working ');
  for i=1:N_DIC
    fprintf(1,'.');
    for j=1:N_Alk
      C_T=DIC_axis(i);
      Alk=Alk_axis(j)+p.SO40*2;
      [pH,Hp,OHm,CO2g,H2CO3,HCO3m,CO3m2,BOH3,BOH4m,HSO4m,SO4m2,H2Saq,HSm] ...
        =solve_carbonate_system_with_sulfur(t,p,Alk,C_T,B_T,SO4_T,H2S_T_ell,T,S,init_guess);
      H2CO3_contour(i,j)=H2CO3;
      CO2g_contour(i,j)=CO2g;
      pH_contour(i,j)=pH;
    end
  end
  fprintf(1,'  done.\n');
  fig1=figure(1);
  clf
  V=[100,150,200:100:800]; [c,h]=contour(DIC_axis,Alk_axis,CO2g_contour',V,'r');
  clabel(c,h,'fontsize',14)
  set(gca,'xtick',min(DIC_axis):25:max(DIC_axis));
  set(gca,'ytick',min(Alk_axis):25:max(Alk_axis));
  hold on
  V=[7.8,8,8.2,8.4];
  [c,h]=contour(DIC_axis,Alk_axis,pH_contour',V,'--b');
  clabel(c,h,'fontsize',14)
  h1=xlabel('DIC');
  h2=ylabel('Alk');
  h3=title('pH (b) & CO2 (r)');
  set([gca,h1,h2,h3],'fontsize',14);
  axis square
  grid
  
  set(fig1, 'PaperUnits', 'inches'); set(fig1, 'PaperSize', [6 6]);
  set(fig1, 'PaperPosition', [0 0 6 6]); % [left, bottom, width, height];
  saveas(fig1,sprintf('Output/CO2g-as-func-of-DIC-and-Alk.pdf'));
  saveas(fig1,sprintf('Output/CO2g-as-func-of-DIC-and-Alk.fig'));

  %% ------------------------------------------------------------------------
  %% calculate and plot carbonate species, CO2 & pH as function of Alk:
  %% ------------------------------------------------------------------------
  %% Zeebe:
  DIC=2050;
  Alk_axis=2250:20:2450;
  %% Stum & Morgan:
  %DIC_axis=0:100:2500;
  %Alk_axis=-500:100:3500;

  N_Alk=length(Alk_axis);
  o.CO2g=zeros(N_Alk,1);
  o.H2CO3=zeros(N_Alk,1);
  o.pH=zeros(N_Alk,1);
  o.Hp=zeros(N_Alk,1);
  o.OHm=zeros(N_Alk,1);
  o.HCO3m=zeros(N_Alk,1);
  o.CO3m2=zeros(N_Alk,1);
  o.BOH3=zeros(N_Alk,1);
  o.BOH4m=zeros(N_Alk,1);
  o.HSO4m=zeros(N_Alk,1);
  o.SO4m2=zeros(N_Alk,1);
  o.H2Saq=zeros(N_Alk,1);
  o.HSm=zeros(N_Alk,1);
  for j=1:N_Alk
      C_T=DIC_axis(i);
      Alk=Alk_axis(j)+p.SO40*2;
      [pH,Hp,OHm,CO2g,H2CO3,HCO3m,CO3m2,BOH3,BOH4m,HSO4m,SO4m2,H2Saq,HSm] ...
        =solve_carbonate_system_with_sulfur(t,p,Alk,C_T,B_T,SO4_T,H2S_T_ell,T,S,init_guess);
      o.CO2g(j)=CO2g;
      o.H2CO3(j)=H2CO3;
      o.pH(j)=pH;
      o.Hp(j)=Hp;
      o.OHm(j)=OHm;
      o.HCO3m(j)=HCO3m;
      o.CO3m2(j)=CO3m2;
      o.BOH3(j)=BOH3;
      o.BOH4m(j)=BOH4m;
      o.HSO4m(j)=HSO4m;
      o.SO4m2(j)=SO4m2;
      o.H2Saq(j)=H2Saq;
      o.HSm(j)=HSm;
  end

  fig2=figure(2);
  clf
  hl1=semilogy(Alk_axis,o.CO2g,'r');
  hold on
  hl2=semilogy(Alk_axis,o.H2CO3,'g');
  hl3=semilogy(Alk_axis,o.HCO3m,'b');
  hl4=semilogy(Alk_axis,o.CO3m2,'c');
  hl5=semilogy(Alk_axis,o.OHm,'m-');
  hl6=semilogy(Alk_axis,o.Hp,'m--');
  set([hl1,hl2,hl3,hl4,hl5,hl6],'linewidth',2)
  h1=xlabel('Alk');
  h2=ylabel(' ');
  h3=title('Response to Alkalinity change');
  h4=legend('CO_2(g)','H_2CO_3','HCO_3^-','CO_3^{2-}','OH^-','H^+' ...
            ,'location','NorthEastOutside');
  set([gca,h1,h2,h3,h4],'fontsize',14);

  
  set(fig2, 'PaperUnits', 'inches'); set(fig2, 'PaperSize', [6 4]);
  set(fig2, 'PaperPosition', [0 0 6 4]); % [left, bottom, width, height];
  saveas(fig2,sprintf('Output/carbonate-response-to-Alk-changes.pdf'));
  saveas(fig2,sprintf('Output/carbonate-response-to-Alk-changes.fig'));

  fid=fopen('Output/carbonate-response-to-Alk-changes.txt','w');
  fprintf(fid,['     CO2g,    H2CO3,       pH,       Hp,      OHm,    ' ...
               'HCO3m,    CO3m2,     BOH3,    BOH4m,    HSO4m,    ' ...
               'SO4m2,    H2Saq,    HSm\n']);  
  for j=1:N_Alk
    fprintf(fid,'%9.4g,',o.CO2g(j),o.H2CO3(j),o.pH(j),...
            o.Hp(j),o.OHm(j),o.HCO3m(j),o.CO3m2(j),...
            o.BOH3(j),o.BOH4m(j),o.HSO4m(j),...
            o.SO4m2(j),o.H2Saq(j),o.HSm(j));
    fprintf(fid,'\n');
  end
  fclose('all');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
function p=pressure_effect_on_K(t,p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% effect of pressure on dissociation and solubility product constants.
%% From Zeebe's csys3

%---------------------- Pressure effect on K's (Millero, 95) ----------%
TC=temperature(t,p); % Celsius!
T=TC+273.15;
P=p.K_sp_pressure;

RGAS = 8.314510;        % J mol-1 deg-1 (perfect Gas)  
R = 83.131;             % mol bar deg-1 
                        % conversion cm3 -> m3          *1.e-6
                        %            bar -> Pa = N m-2  *1.e+5
                        %                => *1.e-1 or *1/10

% index: K1 1, K2 2, Kb 3, Kw 4, Ks 5, Kf 6, Kspc 7, Kspa 8,
%        K1P 9, K2P 10, K3P 11

%----- note: there is an error in Table 9 of Millero, 1995.
%----- The coefficients -b0 and b1
%----- have to be multiplied by 1.e-3!

%----- there are some more errors! 
%----- the signs (+,-) of coefficients in Millero 95 do not
%----- agree with Millero 79


a0 = -[25.5   15.82  29.48  25.60  18.03    9.78  48.76   46. ...
	14.51 23.12 26.57];
a1 =  [0.1271 -0.0219 0.1622 0.2324 0.0466 -0.0090 0.5304  0.5304 ...
	0.1211 0.1758 0.2020];
a2 =  [0.0     0.0    2.608 -3.6246 0.316  -0.942  0.0     0.0 ...
	-0.321 -2.647 -3.042]*1.e-3;
b0 = -[3.08   -1.13   2.84   5.13   4.53    3.91  11.76   11.76 ...
	2.67 5.15 4.08]*1.e-3;
b1 =  [0.0877 -0.1475 0.0    0.0794 0.09    0.054  0.3692  0.3692 ...
	0.0427 0.09 0.0714]*1.e-3;
b2 =  [0.0     0.0    0.0    0.0    0.0     0.0    0.0     0.0 ...
	0.0 0.0 0.0];

for ipc=1:length(a0);
  deltav(ipc)  =  a0(ipc) + a1(ipc).*TC + a2(ipc).*TC.*TC;
  deltak(ipc)  = (b0(ipc) + b1(ipc).*TC + b2(ipc).*TC.*TC);  
  p.lnkpok0(ipc) = -(deltav(ipc)./(R.*T)).*P + (0.5*deltak(ipc)./(R.*T)).*P.*P;
end;

% K1 = K1*exp(p.lnkpok0(1));
% K2 = K2*exp(p.lnkpok0(2));
% Kb = Kb*exp(p.lnkpok0(3));
% Kw = Kw*exp(p.lnkpok0(4));
% Ks = Ks*exp(p.lnkpok0(5));
% Kf = Kf*exp(p.lnkpok0(6));
% Kspc = Kspc*exp(p.lnkpok0(7));
% Kspa = Kspa*exp(p.lnkpok0(8));
% K1P = K1P*exp(p.lnkpok0(9));
% K2P = K2P*exp(p.lnkpok0(10));
% K3P = K3P*exp(p.lnkpok0(11));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p=set_parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  t=0;
  p.scenario=1;
  fprintf(1,'running scenario %2.2d\n',p.scenario);

  %% read scenario parameters:
  %% integration time (yrs):
  p.integration_time=2.5e6;

  % units-related constants: 
  % one mole of C weights 12g; 
  %  ==> 1g C=(1/12) mole ==> 1kg C=(1000/12) mole = (1.e9/12)micro mole
  p.micromole_per_kg_C=1.e9/12;
  p.Carbon_mass_to_molality=60.8;
  p.ocean_mass_kg=13700*1e17; % also roughly equal to volume in liters
  p.ocean_area_meter=p.ocean_mass_kg*1.e-3/4000; %% meter^2

  % parameters:
  p.delta_in=-6;
  p.alpha_remin_frac=0.99;
  p.epsilon=-28; %% fractionation during organic matter production

  %% specified adjustment times in years:
  p.tau_ep=0.3e6;
  p.tau_CaCO3_dissolution=5.0e3; 
  p.tau_mixing=3.e3;
  p.tau_pyrite=3.e3;
  p.tau_siderite=3.e3;

  %% functional form of the forcing function: 'gaussian' or 'step'
  p.H_form='gaussian';

  %% specify characteristics of change during decrease in productivity:
  %% fraction of exported organic matter undergoing sulfate reduction:
  p.t0=1.5e6;
  %% amplitude of perturbation to exported organic matter:
  p.Delta_ep=13.0;
  %% amplitude of perturbation to volcano co2 flux:
  p.Delta_F_in=0.0;
  %% fraction of exported organic matter undergoing sulfate reduction:
  p.alpha_scenario=0.0;
  %% drop in aerobic remineralization as fraction of exported organic matter:
  p.beta_scenario=0.0;
  %% drop in aerobic production as fraction of exported organic matter:
  p.gamma_scenario=1.0;
  %% increase in iron remineralization as fraction of exported organic matter:
  p.eta_scenario=0.0;
  %% an initial perturbation to DIC, mostly for looking at
  %% carbonate system adjustment time scale, a factor that's equal
  %% to a small nondim number.  Set to 0 for no perturbation: 
  p.Delta_M_inorg_initial_perturbation=0.0;

  %% initial conditions:
  p.M_org_init=250.0;
  p.M_inorg_init=35.0;
  p.F_prod0=3.6e-3; % present day value is about 3.6e-2 Gt C/yr
  p.F_remin0=p.alpha_remin_frac*p.F_prod0; % Gt C/yr
  p.F_in0=0.15e-3; % 0.1e-3 Gt C/yr
  p.F_ep0=(1-p.alpha_remin_frac)*p.F_prod0; % Gt C/yr
  p.F_b_inorg0=(1-p.alpha_remin_frac)*p.F_prod0;
  p.F_in_Ca2p_ell0=1.14e-4*p.Carbon_mass_to_molality;

  %% parameters for geochemistry part of the model:
  %% Pressure used to calculate CaCO3 solubility produce constant:
  p.K_sp_pressure=500; %% in atm:
  p.Alk0=2300; % micromole/kg = micro molality
  p.B_T0=400;
  p.SO40=28000;
  p.H2S0=100;
  p.Ca2p0=10000.0; % micromole/kg
  p.Fe3p0=20000;
  p.Fe2p0=200;
  p.Temperature_C_0=5; % Celsius
  p.Salinity0=35;

  %% constant used in parameterization of fraction of volcanic CO2
  %% which goes into the ocean via silicate weathering:
  p.alpha_wthr0=0.8;
  p.alpha_SO4_remin=3.0;
  p.Gamma_Berner=6;
  p.Omega0=1;  % CaCO3 reference degree of saturation
  p.CaCO3_exponent=4.5;

  %% scaling and offseting the carbonate system variables:
  p.carb_sulf_scale=[1.0e-2,1.0,100,10,200,100 ...
          ,100,100 ...
          ,0.01,1,100,100];
  p.carb_sulf_offset=[0,0,300,10,2000,100 ...
          ,300,100 ...
          ,0,p.SO40,0,1000];
  p.carb_sulf_label=[ ...
    'Hp   ';'OHm  ';'CO2g ';'H2CO3';'HCO3m';'CO3m2' ...
    ;'BOH3 ';'BOH4m' ...
    ;'HSO4m';'SO4m2';'H2Saq';'HSm  '];
  p.carb_sulf_init_guess=[0.006, 3.6, 300, 10, 1800, 200, ...
                300, 100, ...
                0, p.SO40, 0, 0];

  indx=0;
  indx=indx+1; p.prognostic_index.delta_org=indx;
  indx=indx+1; p.prognostic_index.delta_inorg=indx;
  indx=indx+1; p.prognostic_index.M_org=indx;
  indx=indx+1; p.prognostic_index.M_inorg=indx;
  indx=indx+1; p.prognostic_index.Alk=indx;
  indx=indx+1; p.prognostic_index.SO4_T=indx;
  indx=indx+1; p.prognostic_index.H2S_T=indx;
  indx=indx+1; p.prognostic_index.Fe3p=indx;
  indx=indx+1; p.prognostic_index.Fe2p=indx;
  indx=indx+1; p.prognostic_index.Ca2p=indx;

  p.prognostic_offset(p.prognostic_index.delta_org)=-28;
  p.prognostic_offset(p.prognostic_index.delta_inorg)=0;
  p.prognostic_offset(p.prognostic_index.M_org)=0;
  p.prognostic_offset(p.prognostic_index.M_inorg)=0;
  p.prognostic_offset(p.prognostic_index.Alk)=2.0*p.SO40;
  p.prognostic_offset(p.prognostic_index.SO4_T)=p.SO40;
  p.prognostic_offset(p.prognostic_index.H2S_T)=p.H2S0;
  p.prognostic_offset(p.prognostic_index.Fe3p)=0;
  p.prognostic_offset(p.prognostic_index.Fe2p)=0;
  p.prognostic_offset(p.prognostic_index.Ca2p)=p.Ca2p0;

  p.prognostic_scale(p.prognostic_index.delta_org)=10;
  p.prognostic_scale(p.prognostic_index.delta_inorg)=10;
  p.prognostic_scale(p.prognostic_index.M_org)=30;
  p.prognostic_scale(p.prognostic_index.M_inorg)=30;
  p.prognostic_scale(p.prognostic_index.Alk)=100;
  p.prognostic_scale(p.prognostic_index.SO4_T)=100;
  p.prognostic_scale(p.prognostic_index.H2S_T)=100;
  p.prognostic_scale(p.prognostic_index.Fe3p)=100;
  p.prognostic_scale(p.prognostic_index.Fe2p)=100;
  p.prognostic_scale(p.prognostic_index.Ca2p)=100;

  p.prognostic_label(p.prognostic_index.delta_org,:)  ='delta_org  ';
  p.prognostic_label(p.prognostic_index.delta_inorg,:)='delta_inorg';
  p.prognostic_label(p.prognostic_index.M_org,:)      ='M_org      ';
  p.prognostic_label(p.prognostic_index.M_inorg,:)    ='M_inorg    ';
  p.prognostic_label(p.prognostic_index.Alk,:)        ='Alk        ';
  p.prognostic_label(p.prognostic_index.SO4_T,:)      ='SO4_T      ';
  p.prognostic_label(p.prognostic_index.H2S_T,:)      ='H2S_T      ';
  p.prognostic_label(p.prognostic_index.Fe3p,:)       ='Fe3p       ';
  p.prognostic_label(p.prognostic_index.Fe2p,:)       ='Fe2p       ';
  p.prognostic_label(p.prognostic_index.Ca2p,:)       ='Ca2p       ';

  % options for fsolve, for carbonate system solution:
  p.optim_options=optimset('MaxIter',1000,'MaxFunEvals',1000 ...
                           ,'DerivativeCheck','off','Display','off' ...
                           ,'Jacobian','on' ...
                           );
  %                 ,'DiffMaxChange',1.e-7 ...
  %                 ,'TolFun',1.e-7,'TolX',1.e-7 ...

  % options for ode solver:
  p.ode_options=odeset('RelTol',1.0e-6,'AbsTol',1.e-6,'MaxStep',1e3);
  
  %% plotting controls:
  p.fig1_plot_panels=[1 0 1 0 0 0 0 0 1];
  p.fig1_panel_order=[1 0 2 0 0 0 0 0 3];
  p.fig1_nrows_ncols=[3 1];
  p.fig1_page_size_inches=[4.5 6]; % width, height
  
  %% logical switches:
  p.do_calc_pCO2=1;
  p.carbonate_system_solver='sulfur';  % 'standard', 'sulfur' or 'approx'
  p.do_CaCO3_dissolution=1;
  p.do_interactive_alpha_wthr=1;
  p.do_write_restart=1;
  p.do_read_restart=1;
  p.do_use_Ksp_aragonite=0;
  p.do_constant_Alk=0;
  p.do_prognostic_Ca2p=1;
  p.do_sulfate_remineralization=0;
  p.do_H2S_mixing_and_oxidation=0;
  p.do_pyrite_formation=0;
  p.do_iron_remineralization=0;
  p.do_Fe2p_mixing_and_oxidation=0;
  p.do_siderite_formation=0;
  p.do_run_model=1;
  p.do_save_figures_as_pdf=1;
  p.do_save_figures_as_Fig_files=1;

  %% add to p the effects of pressure on solubility and
  %% dissociation constants:
  p=pressure_effect_on_K(t,p);

