% Written by Katja Fennel (katja.fennel@dal.ca)
% September, 5 2008

% load HOT data
load -ascii HOT_CTD_potential_density_0_1000m_noheadder.txt
depth = HOT_CTD_potential_density_0_1000m_noheadder(:,2);
dens = HOT_CTD_potential_density_0_1000m_noheadder(:,3);
% remove bad data points
bad = find(dens<0);
depth(bad) = [];
dens(bad) = [];

% plot density profile
figure
plot(dens,-depth,'.')
xlabel('density [kg/m3]')
ylabel('depth [m]')

% calculate mean density profiles (note that data are in a vertical intervall of 2 m)
dz=2;
depthvec = [0:dz:1000];
for k=1:length(depthvec)
    ind = find(depth==depthvec(k));
    mean_dens(k) = mean(dens(ind));
end
% and plot
hold on
plot(mean_dens,-depthvec,'o-r')

% Target amounts of C to be sequesterd:
C2Bsequ = [0.01 0.05 0.1 0.5 1]; % amounts of C [Gt] to be sequesterd  
% one Gt of C in mol C:
oneGtonC = 8.3E13; % mol C
% same in kg C
C2Bsequ_in_kgC = C2Bsequ*oneGtonC*12/1000;
C2Bsequ_in_metrictonC = C2Bsequ*oneGtonC*12/1000000;

% Divide amount of C to be sequestered by area of NPSG to determine 
% how much C would have to be sequestered per m2:
area = 34E12; % m2
sequC = C2Bsequ*oneGtonC/area; % mol C/m2 to be sequestered

% data from table 1 in Karl&Letelier (2008)
z_upw = [100 200 300 400 500 750 1000];
Csequ_potential = [0.0152 0.002 0.0327 0.0545 0.1218 0.1418 0.1149]; 

% density of upwelled water:
dens_upw = mean_dens(round(z_upw/dz));

% loop over different source waters
for z=1:7
    sw2Bupw = sequC/Csequ_potential(z); % seawater to be upwelled in m3
    % loop over different C amounts to be sequestered
    for j=1:length(sw2Bupw)
       sw_upw = sw2Bupw(j);
       % iteration
       k=1; 
       new_dens = mean_dens;
       xaccum = 0;
       while sw_upw>0
         % density of next underlying layer is dens_crit
         dens_crit = mean_dens(k+1);
         % volume and density of current "mixed layer"
         V=dz*k;
         V=V+xaccum;
         dens_curr = new_dens(k);
         % calculate how much of the dense upwelled water can be mixed in before
         % reaching dens_crit
         x = V*(dens_crit-dens_curr)/ (dens_upw(z)-dens_crit); % x in m3
         xaccum=xaccum+x;
         % mix in
         for i=1:k
           new_dens(i) = dens_crit;
         end
         sw_upw=sw_upw-x;
         k=k+1;
       end
       fprintf('To sequester %2.3f Gt C by upwellng from %4.1f m the top %5.1f would be mixed.\n', C2Bsequ(j),z_upw(z), dz*k+xaccum)
       if (z==4 & j==5)
           % plot result for sequestration of 1 Gt C from 400 m
           figure
           plot(mean_dens,-depthvec,'-k')
           hold on
           plot(new_dens,-depthvec,'--k')
           xlabel('Potential density [kg/m3]')
           ylabel('Depth [m]')
           %title(['To sequester ',num2str(C2Bsequ(j)),'Gt of C, the NPSG would have to be homogenized down to ',num2str(dz*k+xaccum),' m']) 
           set(gcf,'PaperPosition',[0.25 2.5 6 4])
           %print(gcf,'-dpng','Figure_1xx')
       end
    end
end




       