% pHprofile.m                                       by:  Edward T Peltzer, MBARI
%                                                   revised:  2006 Feb 15.
%
% M-file to plot pH profiles calculated from TCO2, TAlk data on two pH scales.
%
% Data is from the WOCE CO2 Program for P17N Station 10.
%
% pH was calculated using CO2SYS.EXE by Ernie Lewis and Doug Wallace at BNL.
% Mehrbach's apparant CO2 constants for a seawater matrix were used throughout.
% Calculated pH includes corrections for total P and total Si.
%
% Data file = HydroData.xls  /  Worksheet = pHandTau
%
%   Column   Parameter        Units
%
%     D      CTD pressure     dbar  
%     F      pH (NBS)         pH on the NBS / NIST scale.
%     H      pH (TOT)         pH on the Total H+ scale.


% Clear environment space and set orientation

clear all
clf

orient tall


% Load data and assign data vectors

LAT = 38.23192;

prs = xlsread('HydroData','pHandTau','D2:D29');
pHnbs = xlsread('HydroData','pHandTau','F2:F29');
pHtot = xlsread('HydroData','pHandTau','H2:H29');

dep = DEPTH(prs,LAT);


% Fit lines to data

depX = [0:25:3950];

snbs = ppval(spline(dep,pHnbs),depX);
stot = ppval(spline(dep,pHtot),depX);


% Define plotting region and make plot

left = 1.25/8.5;
bottom = 2.5/11;
width = 6/8.5;
height = 6/11;

subplot('position',[left bottom width height])


% Make plot

plot(pHnbs,dep,'bo',snbs,depX,'LineWidth',1,'Color',[0 0 1])

hold on

plot(pHtot,dep,'ro',stot,depX,'LineWidth',1,'Color',[1 0 0])


% Define axis and set tick-label-size

axis('ij');
axis([7.4 8.4 -100 4000]);

h=gca;
set(h,'FontSize',16);
set(h,'FontWeight','normal');
set(h,'LineWidth',1);


% Manually define and control X-ticks and X-tick labels

xtl = [7.5:0.1:8.3];

N = length(xtl);

for I = 1:N
    if rem(I,2) == 0
        if xtl(I) == 8
    	    XTL(I,:) = ['8.0'];
        else
            XTL(I,:) = [num2str(xtl(I))];
        end
    end
end

set(h,'XTick',xtl);
set(h,'XTickLabel',XTL);


% Manually define and control Y-ticks and Y-tick labels

ytl = [0:250:4000];

N = length(ytl);

for I = 1:N
    if rem(I-1,4) == 0
        if ytl(I) == 0
            YTL(I,:) = ['    0'];
        elseif ytl(I) < 1000
        	YTL(I,:) = ['  ' num2str(ytl(I))];
        else
            YTL(I,:) = [' ' num2str(ytl(I))];
        end
    end    
end

set(h,'YTick',ytl);
set(h,'YTickLabel',YTL);


% Label plot -- X-Axis

xlabel(['pH']);
h=get(gca,'xlabel');
set(h,'FontSize',20);
set(h,'FontWeight','normal');


% Label plot -- Y-Axis

ylabel('Depth (m)');
h=get(gca,'ylabel');
set(h,'FontSize',20);
set(h,'FontWeight','normal');


% Label plot -- Title

title(['WOCE P17N Sta10 {\itin situ} data']);
h=get(gca,'title');
set(h,'FontSize',20);
set(h,'FontWeight','normal');


% Label plot -- Legends

leg1 = ['-o- = NBS scale'];
h = text(8,1000,leg1);
set(h,'Color',[0,0,1]);
set(h,'FontSize',16);
set(h,'FontWeight','normal');

leg2 = ['-o- = Total H^{+}'];
h = text(8,1250,leg2);
set(h,'Color',[1,0,0]);
set(h,'FontSize',16);
set(h,'FontWeight','normal');

hold off
