% PS06p1_main Main m-file program for integrating the five-box, global 
%   model of phosphorus distributions. The model idea was "borrowed" from
%   Broecker and Peng's (1982) super-problem #7.
%
% Started 2004:10:21 D. Glover, WHOI
% Modif'd 2008:11:06 DMG to change the names to PS06_1

% Start the stopwatch

t_init_ps06p1=clock;

% Load colormap

load purp_map4;

% Declare some GLOBAL variables

global TAU REFF V SPYR;

SPYR=86400*365;                         % secs per year
TAU=10*SPYR;                            % Residence time scale (yrs) in surface box
REFF=0.99;                              % Remineralization efficiency
V=[3e16 1.5e17 9e16 8.1e17 2.7e17];     % Volumes

Pinit=[0 0 0 0 0];      % Initial phosphorus concentrations
Trang=[0 5e5];          % Integration time range

% Do the initial run

[T P]=ode15s('phos5box',Trang,Pinit);

% Calculate total phosphorus evolving with time

Ptot=sum(P');

% Plot inital results

figure(1);

subplot(321);plot(T,P(:,1));    % SAT
title('Surface Atlantic');
subplot(322);plot(T,P(:,3));    % SIP
title('Surface Indo-Pacific');
subplot(323);plot(T,P(:,5));    % DAT
title('Deep Atlantic');
subplot(324);plot(T,P(:,4));    % DIP
title('Deep Indo-Pacific');
subplot(325);plot(T,P(:,2));    % ANT
title('Antarctic');
subplot(326);plot(T,Ptot);      % Evolution of total P
title('Total Phosphorus');

format compact;
disp(' ');
disp('The final phosphorus concentrations for each box:');

format short g
fin=length(P);
P(fin,:)

disp(' ');
disp('The global average phosphorus cconcentration:');
P_avg=sum(P(fin,:).*V)/sum(V)

disp(' ');
disp('Broecker''s horizontal segregation diagnostic:');
horzseg=(P(fin,4)-P(fin,3))/(P(fin,5)-P(fin,1))

% Start sensitivity analysis

% Initialize diagnostics

pavg=zeros(10,10);
horz=zeros(10,10);
ratatl=zeros(10,10);
ratpac=zeros(10,10);

%tauvec=[0.1 0.3 0.5 1.0 3.0 5.0 10.0 30.0 50.0 100.0];
%tauvec=[1e-2 1e-1 1e0 1e1 1e2 1e3 1e4 1e5 1e6 1e7];
tauvec=[3e0 1e1 3e1 1e2 3e2 1e3 3e3 1e4 3e4 1e5];
%Reffvec=[0.95 0.955 0.96 0.965 0.97 0.975 0.98 0.985 0.99 0.995];
Reffvec=[0.91 0.92 0.93 0.94 0.95 0.96 0.97 0.98 0.99 0.999];
%Reffvec=[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9];

for i=1:10;
    REFF=Reffvec(i);
    for j=1:10;
        TAU=tauvec(j)*SPYR;
        [T,P]=ode15s('phos5box',Trang,Pinit);
        fin=length(P);
        pavg(i,j)=sum(P(fin,:).*V)/sum(V);
        horz(i,j)=(P(fin,4)-P(fin,3))/(P(fin,5)-P(fin,1));
        ratatl(i,j)=P(fin,1)/P(fin,5);
        ratpac(i,j)=P(fin,3)/P(fin,4);
        PsurfAtl(i,j)=P(fin,1);
        PdeepAtl(i,j)=P(fin,5);
        PsurfPac(i,j)=P(fin,3);
        PdeepPac(i,j)=P(fin,4);
        Antarctic(i,j)=P(fin,2);
    end
end

% Make contour plots

figure(2);colormap(purp_map4);
contourf(log10(tauvec),Reffvec,pavg);colorbar;
title('Average P Concentration (mmol m^{-3})','fontsize',16);

figure(3);colormap(purp_map4);
contourf(log10(tauvec),Reffvec,horz);colorbar;
title('Log_{10}(Horizontal Segregation) (Broecker and Peng, 1982)','fontsize',16);

figure(4);colormap(purp_map4);
contourf(log10(tauvec),Reffvec,ratatl);colorbar;
title('Log_{10}(P_{surf}:P_{deep}) Atlantic','fontsize',16);

figure(5);colormap(purp_map4);
contourf(log10(tauvec),Reffvec,ratpac);colorbar;
title('Log_{10}(P_{surf}:P_{deep}) Pacific','fontsize',16);

figure(6);colormap(purp_map4);
subplot(221);contourf(log10(tauvec),Reffvec,log10(PsurfAtl));colorbar;
title('Surface Atlantic P')
subplot(222);contourf(log10(tauvec),Reffvec,log10(PsurfPac));colorbar;
title('Surface Pacific P')
subplot(223);contourf(log10(tauvec),Reffvec,log10(PdeepAtl));colorbar;
title('Deep Atlantic P')
subplot(224);contourf(log10(tauvec),Reffvec,log10(PdeepPac));colorbar;
title('Deep Pacific P')

figure(7);colormap(purp_map4);
contourf(log10(tauvec),Reffvec,log10(Antarctic));colorbar;
title('Antarctic P');

% Report total processing execution time

sprintf('%s%g\n','Total Elasped Time: ',etime(clock,t_init_ps06p1))