
% combine CTD and XBT data into a single dataset of temperature,
% salinity and oxygen.
clear all; close all;

% load XBT data
load xbt
% drop zonal section from 15N, 38W to 20.34N, 23.51W
q=find(lon<-23.5);
lon(q)=[]; lat(q)=[]; temp(:,q)=[];
xtemp=temp; xlon=lon; xlat=lat;
% calculate pressure
xp=NaN*ones(length(z),length(xlat));
for i=1:length(xlat);
    xp(:,i)=sw_pres(z,xlat(i));
end;
% allocate for estimated salinity, oxygen
xsal=NaN*ones(size(xtemp));
xoxy=NaN*ones(size(xtemp));

% load CTD data (see "convert_ctd.m" in this directory:)
load 'C:\Documents and Settings\lumpkin.AOML\My Documents\data\pne06\ctd\pne06\ctddata';
temp=fillnans(temp,0);
sal=fillnans(sal,0);
oxy=fillnans(oxy,0);
ptmp=sw_ptmp(sal,temp,p,0);
sigthet=sw_pden(sal,ptmp,p,0)-1000;

% estimate xsal (xbt salinity)
xDistMin=NaN*ones(length(xlat),1);
xsal=NaN*ones(size(xtemp));
xsigthet=xsal;
xoxy=xsal;
for i=1:length(xlon);   % loop through XBT
    % find closest CTD casts
    del=NaN*ones(length(lon),1);
    for j=1:length(lon);
        del(j)=sw_dist([lat(j) xlat(i)],[lon(j) xlon(i)],'km');
    end;
    xDistMin(i)=min(del);
    % don't bother if there isn't a close CTD
    q=find(del<200);
    if (length(q))
        del=del(q);
        wts=1./(del.^2);
        ttemp=temp(:,q);
        tsal=sal(:,q);
        % loop through pressure bins
        for j=1:size(xp,1);
            if (isfinite(xtemp(j,i)))
                % identify "close" pressures (within 50 dbar,
                % assuming resolution of p, xp is )
                pind=max([1 j-25]):min([length(xp),j+25]);
                % translate to CTD pressure index
                pind=find(p>=min(xp(pind,i)) & p<=max(xp(pind,i)));
                % do a polyfit of T/S for each close cast
                xsfit=NaN*ones(size(tsal,2),1);
                for k=1:length(xsfit);
                    [pf,sf]=polyfit(ttemp(pind,k),tsal(pind,k),1);
                    xsfit(k)=polyval(pf,xtemp(j,i));
                end;
                % take average, weighting by 1/distance
                qq=find(isfinite(xsfit));
                xsal(j,i)=sum(wts(qq).*xsfit(qq))./sum(wts(qq));
            end; % if isfinite(xtemp) at this depth
        end; % loop through pressure bins
        
        % plot results of salinity estimate
        %clf;
        %for k=1:length(xsfit);
        %    plot(tsal(:,k),ttemp(:,k),'color',[.5 .5 1]);
        %    hold on;
        %end;
        %plot(xsal(:,i),xtemp(:,i),'k');
        %hold off;
        
    end; % if there's a CTD cast close enough
end; % loop through XBT drops

% calculate sigma-theta for XBT
xptmp=sw_ptmp(xsal,xtemp,xp,0);
xsigthet=sw_pden(xsal,xptmp,xp,0)-1000;

% calculate oxygen for XBT, assuming isopycnal flow
%% FUTURE VERSION: just linearly interpolate for mixed layer.
for i=1:length(xlon);   % loop through XBT
    % find closest CTD casts
    del=NaN*ones(length(lon),1);
    for j=1:length(lon);
        del(j)=sw_dist([lat(j) xlat(i)],[lon(j) xlon(i)],'km');
    end;
    xDistMin(i)=min(del);
    % don't bother if there isn't a close CTD
    q=find(del<200);
    if (length(q))
        del=del(q);
        wts=1./(del.^2);
        tsigthet=sigthet(:,q);
        toxy=oxy(:,q);
        % loop through pressure bins
        for j=1:size(xp,1);
            if (isfinite(xsigthet(j,i)))
                % find closest sigthet values to xsigthet
                xsfit=NaN*ones(size(toxy,2),1);
                for k=1:length(xsfit);
                    dd=(tsigthet(:,k)-xsigthet(j,i)).^2;
                    dd=find(dd==min(dd));
                    xsfit(k)=toxy(dd,k);
                end;
                % take average, weighting by 1/distance
                qq=find(isfinite(xsfit));
                xoxy(j,i)=sum(wts(qq).*xsfit(qq))./sum(wts(qq));
            end; % if isfinite(xtemp) at this depth
        end; % loop through pressure bins
        
        
    end; % if there's a CTD cast close enough
end; % loop through XBT drops

% interpolate xtemp, xsal, xoxy onto p
xtemp2=NaN*ones(length(p),length(xlat));
xsal2=xtemp2;
xoxy2=xtemp2;
for i=1:length(xlat);
    xtemp2(:,i)=interp1(xp(:,i),xtemp(:,i),p,'linear');
    xsal2(:,i)=interp1(xp(:,i),xsal(:,i),p,'linear');
    xoxy2(:,i)=interp1(xp(:,i),xoxy(:,i),p,'linear');
end;
xtemp=xtemp2; clear xtemp2;
xsal=xsal2; clear xsal2;
xoxy=xoxy2; clear xoxy2;
clear xp i j k xsigthet xptmp dd del pf pind ptmp q qq
clear sigthet tsigthet toxy tsal ttemp wts 
clear xsfit z sf

% merge
lat=[lat;xlat];
lon=[lon;xlon];
oxy=[oxy xoxy];
sal=[sal xsal];
temp=[temp xtemp];
[lat,ind]=sort(lat);
lon=lon(ind);
oxy=oxy(:,ind);
sal=sal(:,ind);
temp=temp(:,ind);

% save results
save ctd_xbt_merged lon lat temp sal oxy p;

% plot results
plotmerged