
%plot-oxygen
clear all
o2_params
p_ylim=[0 3500];
p_filter=9; %median and running filter over 9 points (m)

%
%Location of final CTD data
%
rootdir='./';
dirctdf=[rootdir 'Caliprofiles/'];
%
% output directory
%

dirop=[rootdir 'FilteredOP/'];

ttlstr=[];

giprof=1:length(gprof2get);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% END PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pname={'temp';'sali';'oxyg'};
pfieldname={'dtemp68i';'dsalii';'oxcali'}; %variable ID in calibrated CTD files
punits={'^\circle C';'PSU';'umol/kg'};
%
% LOAD DATA
%
nprof=length(giprof);
actd.dpres=(0:5:6000)';
ndep=length(actd.dpres);
actd.stat=floor(gprof2get(giprof)/100);
actd.cast=gprof2get(giprof)-100*actd.stat;
actd.dtemp68if=NaN*ones(ndep,nprof);
actd.dsaliif=NaN*ones(ndep,nprof);
actd.oxcalif=NaN*ones(ndep,nprof);

% FIRT SUBSAMPLE MAIN STATIONS
% for iprofindex=1:length(giprof)
%     istat=actd.stat(iprofindex)
%     icast=actd.cast(iprofindex)
%     eval(sprintf('load %s/ox_cal_%s_STN%03d_CAST%02d.mat ',...
%         dirctdf,cname,istat,icast))
%     actd.latx(iprofindex)=ctd.latx;
%     actd.lonx(iprofindex)=ctd.lonx;
%     actd.maxd(iprofindex)=max(find(~isnan(ctd.dtemp68i)));
%     %
%     % filter and subsamples every 5m
%     %
%     for iprop=1:length(pname)
%         eval(sprintf('propi=ctd.%s;',pfieldname{iprop}))
%         eval(sprintf('propi_qc=ctd.%s_qc;',pfieldname{iprop}))
%         %keeps only valid data
%         gibad=find((propi_qc~=2)&(propi_qc~=3)&(propi_qc~=6)); %only takes good/doubtful/interpolated data
%         propi(gibad)=NaN;
%         figure(10);clf;plot(propi,ctd.dpres)
%         %Median filter
%         propi=median_filter(propi,p_filter);
%         hold on; plot(propi,ctd.dpres,'k')
%         propi=filtend(propi,p_filter);
%         %Hanning filter
%         plot(propi,ctd.dpres,'r')
%         propii=interp1(ctd.dpres, propi, actd.dpres);
%         %tricky qc
%         [lia,locb]=ismember(actd.dpres,ctd.dpres);
%         propii_qc=6*ones(size(actd.dpres)); %6 by default because if we have no pressure/no data it will have been interpolated in final
%         propii_qc(find(lia))=propi_qc(locb(find(lia)));
%         %complete with flag 9 where remaining NaN (typically beginning or
%         %end of profile)
%         propii_qc(isnan(propii))=9;
%         plot(propii,actd.dpres,'mo');set(gca,'ydir','rev')
%         eval(sprintf('actd.%sf(1:length(propii),iprofindex)=propii;',pfieldname{iprop}))
%         eval(sprintf('actd.%sf_qc(1:length(propii),iprofindex)=propii_qc;',pfieldname{iprop}))
%         %ppause
%         pause(1)
%     end
% end %iprofindex

%SECOND CREATE INDIVIDUAL OUTPUT FILES
%eval(['cd ' dirctdf])
lsstr=ls([dirctdf 'ox_cal_' cname '_STN*.mat']);
ll=4+min(strfind(lsstr, '.mat'));
lsstr(ll);
lsstr=reshape(lsstr,ll,length(lsstr)/ll)';
iprofindex=0;
for ifile=1:size(lsstr,1)
    clear ctd
    disp(lsstr(ifile,:))
    eval(['load ' lsstr(ifile,:)])
    iprofindex=iprofindex+1;
    actd.latx(iprofindex)=ctd.latx;
    actd.lonx(iprofindex)=ctd.lonx;
    actd.maxd(iprofindex)=max(find(~isnan(ctd.dtemp68i)));     %
    % filter and subsamples every 5m
    %
    for iprop=1:length(pname)
        eval(sprintf('propi=ctd.%s;',pfieldname{iprop}))
        eval(sprintf('propi_qc=ctd.%s_qc;',pfieldname{iprop}))
        %keeps only valid data
        gibad=find((propi_qc~=2)&(propi_qc~=3)&(propi_qc~=6)); %only takes good/doubtful/interpolated data
        propi(gibad)=NaN;
        figure(10);clf;
        plot(propi,ctd.dpres)
        %Median filter
        propi=median_filter(propi,p_filter);
        hold on; plot(propi,ctd.dpres,'k')
        propi=filtend(propi,p_filter);
        %Hanning filter
        plot(propi,ctd.dpres,'r')
        propii=interp1(ctd.dpres, propi, actd.dpres);
        %tricky qc
        [lia,locb]=ismember(actd.dpres,ctd.dpres);
        propii_qc=6*ones(size(actd.dpres));%6 by default because if we have no pressure/no data it will have been interpolated in final
        propii_qc(find(lia))=propi_qc(locb(find(lia)));
        %complete with flag 9 where remaining NaN (typically beginning or
        %end of profile)
        propii_qc(isnan(propii))=9;
        plot(propii,actd.dpres,'mo');set(gca,'ydir','rev','ylim',p_ylim)
        drawnow
        
        ctd.dpresf=actd.dpres;
        eval(sprintf('ctd.%sf=propii;',pfieldname{iprop}))
        eval(sprintf('ctd.%sf_qc=propii_qc;',pfieldname{iprop}))
        eval(sprintf('actd.%sf(1:length(propii),iprofindex)=propii;',pfieldname{iprop}))
        eval(sprintf('actd.%sf_qc(1:length(propii),iprofindex)=propii_qc;',pfieldname{iprop}))
        %ppause
    end
    ctdf.stat=ctd.stat;
    ctdf.cast=ctd.cast;
    ctdf.dayd=ctd.dayd;
    ctdf.latx=ctd.latx;
    ctdf.lonx=ctd.lonx;
    ctdf.dpres=ctd.dpresf;
    ctdf.dtemp90if=ctd.dtemp68if/1.00024;
    ctdf.dtemp90if_qc=ctd.dtemp68if_qc;
    ctdf.dsaliif=ctd.dsaliif;
    ctdf.dsaliif_qc=ctd.dsaliif_qc;
    ctdf.oxcalif=ctd.oxcalif;
    ctdf.oxcalif_qc=ctd.oxcalif_qc;
    %now erase ctd/replace to have consistent variable names
    ctd=ctdf;
    
    [pathstr, fname, ext] = fileparts(lsstr(ifile,:));
    disp(['save ' dirop fname '.mat ctd'])
    eval(['save ' dirop fname '.mat ctd'])
    disp('****************************************')
end

actd.latx=actd.latx(:);
actd.lonx=actd.lonx(:);
actd.maxd=actd.maxd(:);
actd.dtemp90if=actd.dtemp68if/1.00024;
actd.dtemp90if_qc=actd.dtemp68if_qc;
   
save all_stations_singleprofile_filtered.mat actd p_filter
   

