function varargout=despike1(cruiseid,staa,stab)
% DESPIKE1 -  remove spikes in the primary sensors from CTD Data    
%
% CTD Calibration toolbox.
%
% USAGE
%	cnv_ds1=despike1(cruiseid,staa,stab);
%
% INPUT
%	cruiseid	cruise identification name
%   
% OUTPUT
%	cnv_ds1  
%          
% DESCRIPTION: This function uses the cnv files from the cruise database
% and remove the spikes related to the primary sensors
%

% AUTHORS:Carlos Fonseca 
%         UM/CIMAS
%         Fri Dec 2 09:19:37 EDT 2005
% 
% VERSION: 
%	$Revision: 1.13 $
%	$Date: 2008/04/17 00:10:38 $
%	$Id: despike1.m,v 1.13 2008/04/17 00:10:38 duncombe Exp duncombe $
%
% CHANGELOG: 
%   23-Jul-2004, Version 1.0
%		* Initial version.
%   24-Mar-2006 Molly Baringer
%       * Fixed comments.  Reads cnv_fill and incrementally saves to
%       cnv_ds1
%
%  2007/12/12 CMDR
%  	* Get cruiseid from a file containing the identifier in the data
%  	directory.
%  2007/12/21 CMDR
%	* test for distinct x-values in interpolation
%  2007/12/31
%	* plot dots as well as lines, makes choosing the desired 
% 	pressures easier 
%  2008/03/11
%	* enter station range in the argument list.
%	* provide defaults if staa or stab are empty.
%	* provide default for cruiseid if its empty.
%  2008/03/12 CMDR
%	* interp dies if the points for interpolation are chosen in an 
%	unexpected order - fixed
%  2008/03/27 CMDR
%	* fixed the assumption that the dir command would provide a nicely
%	sorted array of station numbers. Now sort station numbers and loop 
%	through the station number index.
%  2008/04/16 CMDR
%	* incorporate beckman sensor. 
%	* fiddle around a bit with dynamic field names
%

% PROGRAMMING NOTES
% I don't know if using dynamic fieldnames to cope with the voltage and 
% current differences in the types of oxygen sensor is such a good idea. 
% Certainly it helps to cut down on the testing if fieldnames are present,
% but if this kind of scheme had been in place when I was trying to figure
% out what these programs were doing and trying to get them to work with 
% missing secondary sensors etc, I don't know if I would ever have 
% sorted it out. So from clarity POV it might be better to stick to the 
% endless repetitious testing for fieldnames.

%%%%%%%%%%% CMDR
if nargin==0, cruiseid=getcruiseid; end;
if isempty(cruiseid), cruiseid=getcruiseid; end;
%%%%%%%%%%%
% Retrieve the necessary preferences for this cruise
group = group_name(cruiseid);

% cruise_dir is the main calibration data directory 
if ispref(group,'cal_data_dir'),
	cruise_dir=getpref(group,'cal_data_dir');
else
	error(['The cal_data_dir preference was not set for cruise: ' cruiseid '.  Run register_cruise.m']);
end;

group_pref = getpref(group);

% Load data from the cruise database
db_file = fullfile(cruise_dir,[cruiseid '_db.mat']);
if exist(db_file)==2
    ld_file=(['load ',db_file]);
    eval(ld_file)
else
	error(['The cruise data base file was not found on the path.  ' db_file]);
end

n=who('cnv_fill');
if isempty(n),
    error(['cnv_fill not found. run fill_surface.m ']);
    return;
end;

%% find which field name we must use for oxygen. This kludgy mess is because
%% I cannot think of a more elegant way to find out which fieldname is
%% present. There should be a oneliner way to do it but I can't find it.
if isfield(cnv_fill,'oxvpri'), fldoxvpri='oxvpri'; end;
if isfield(cnv_fill,'oxcpri'), fldoxvpri='oxcpri'; end;
%%  might be more understandable....
% fields=fieldnames(cnv_fill);
% [s,f,t]=regexp([fields{:}],'(oxcpri)|(oxvpri)');
% fldoxvpri=fields(s:f);

% we need to be very sure which order the stations are being accessed, so 
% set up aux_stations so that we do not need to rely on any arbitrary order
% that the variable structure might have been created in.
[aux_stations,ind_stns]=sort(vertcat(cnv_fill.station));

disp('****************************************** ')
disp('DESPIKE PRIMARY SENSORS')
disp('****************************************** ')
disp('The data found in cnv goes from') 
disp(['the station ',num2str(min(aux_stations)),' to the station ', ...
	num2str(max(aux_stations)),'.'] ) 
disp('****************************************** ')
%
%%%%%%
if ~exist('staa'), staa=input('Input first station: '); end;
if ~exist('stab'), stab=input('Input last station: '); end;
if isempty(staa), staa=min([cnv_fill.station]); end;
if isempty(stab), stab=max([cnv_fill.station]); end;
%%%%%%%

% Assumes stations are in "order" in cnv data array (which they should be
% since they were loaded with a "dir" command that puts things in order.
% And it turns out they may not be in the expected order. So aux_stations 
% has been rearranged (sorted) above so that we know the order.
sta_min_idx=find(staa==aux_stations);
sta_max_idx=find(stab==aux_stations);

%
% initial setup
% loop through the known order of ind_stns
for i=ind_stns(sta_min_idx:sta_max_idx).',
    cnv_ds1(i) = cnv_fill(i);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% JH Temp interp %%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

         %cnv_dlp=cnv_dlp;
     tt=find(cnv_ds1(i).temppri~=0 & cnv_ds1(i).temppri>=0);
 
	% sal==0 ;  all pr >= first good sal ; all pr <= last good sal
    % but only if there are good dalinity
    if ~isempty(tt),
      hh=find(abs(cnv_ds1(i).temppri==0 | cnv_ds1(i).temppri<-2) | isnan(cnv_ds1(i).temppri)   & ...                             
          cnv_ds1(i).pressure >= cnv_ds1(i).pressure(tt(1)) & ...         
          cnv_ds1(i).pressure <= cnv_ds1(i).pressure(tt(length(tt))));
      if length(hh) > 0,
	if isfield(cnv_ds1,'temppri'),
          cnv_ds1(i).temppri(hh)=interp1(cnv_ds1(i).pressure(tt), ...
		  cnv_ds1(i).temppri(tt),cnv_ds1(i).pressure(hh));
	end;
    if isfield(cnv_ds1,'thetapri'),
          cnv_ds1(i).thetapri(hh)=interp1(cnv_ds1(i).pressure(tt), ...
		  cnv_ds1(i).thetapri(tt),cnv_ds1(i).pressure(hh));
	end;
          for n=1:length(hh);
              dum=int2str(cnv_ds1(i).quality(hh(n)));
              if dum(6)=='2'
                  cnv_ds1(i).quality(hh(n))=cnv_ds1(i).quality(hh(n))+40;
              end; % if dum(6)=='2'
%%%%%%%%%%%%%%%%%%%%%%%%%%%% JH qq %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%              
              dm=int2str(cnv_ds1(i).quality1tem(hh(n)));
              if dm == '2'
                  cnv_ds1(i).quality1tem(hh(n))=...
                      cnv_ds1(i).quality1tem(hh(n))+4;
              end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% JH qq %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%              
          end; % for n=1:length(ii);
      end; % if length(ii) > 0,
    end; % if ~isempty(jj),
    clear tt hh dum dm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%% OXYGEN %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% primary oxv despiking: linearly interpolate over zeros; WOCE flag to 6
    % All good O2
    	jj=find(cnv_ds1(i).(fldoxvpri)~=0 & cnv_ds1(i).(fldoxvpri)>=0);
	% O2==0 ;  all pr >= first good O2 ; all pr <= last good O2
    % but only if there are good oxygens
    if ~isempty(jj),
      ii=find(abs(cnv_ds1(i).(fldoxvpri)==0 | cnv_ds1(i).(fldoxvpri)<0) | isnan(cnv_ds1(i).(fldoxvpri))   & ...                             
          cnv_ds1(i).pressure >= cnv_ds1(i).pressure(jj(1)) & ...         
          cnv_ds1(i).pressure <= cnv_ds1(i).pressure(jj(length(jj))));
      if length(ii) > 0,
	if isfield(cnv_ds1,'oxvpri'),
          cnv_ds1(i).oxvpri(ii)=interp1(cnv_ds1(i).pressure(jj), ...
		  cnv_ds1(i).oxvpri(jj),cnv_ds1(i).pressure(ii));          
	end;
    if isfield(cnv_ds1,'oxmllpri'),
          cnv_ds1(i).oxmllpri(ii)=interp1(cnv_ds1(i).pressure(jj), ...
		  cnv_ds1(i).oxmllpri(jj),cnv_ds1(i).pressure(ii));    
    end;
    if isfield(cnv_ds1,'oxmmpri'),
          cnv_ds1(i).oxmmpri(ii)=interp1(cnv_ds1(i).pressure(jj), ...
		  cnv_ds1(i).oxmmpri(jj),cnv_ds1(i).pressure(ii));    
    end;
          for n=1:length(ii);
              dum=int2str(cnv_ds1(i).quality(ii(n)));
              if dum(6)=='2'
                  cnv_ds1(i).quality(ii(n))=cnv_ds1(i).quality(ii(n))+40;
              end; % if dum(6)=='2'
%%%%%%%%%%%%%%%%%%%%%%%%%%%% JH qq %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
              dm=int2str(cnv_ds1(i).quality1oxy(ii(n)));
              if dm == '2'
                  cnv_ds1(i).quality1oxy(ii(n))=cnv_ds1(i).quality1oxy(ii(n))+4;
              end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% JH qq %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%              
          end; % for n=1:length(ii);
      end; % if length(ii) > 0,
    end; % if ~isempty(jj),
    
    
   
   % primary oxv despiking: use cleaner program to despike deep water;
   % WOCE flag to 6
   if 1,
    if length(cnv_ds1(i).pressure)>1000;
        ii=find(cnv_ds1(i).pressure>1000);
        factor = 1;
        order = 2;
        [pc,pm]=cleaner(cnv_ds1(i).(fldoxvpri)(ii),factor,order);
	
        jj=find(pm==0);
        for n=1:length(jj),
       
           if isfield(cnv_ds1,'oxvpri'),
              qualtest=abs(cnv_ds1(i).oxvpri(ii(jj(n)))-pc(jj(n)))>0.0004;
           end;
           
           if qualtest,
              cnv_ds1(i).(fldoxvpri)(ii(jj(n)))=pc(jj(n));
             
                dum=int2str(cnv_ds1(i).quality(ii(jj(n))));
                if dum(6)=='2',
                    cnv_ds1(i).quality(ii(jj(n)))= ...
                    cnv_ds1(i).quality(ii(jj(n)))+40;
                end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% JH qq %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                dm=int2str(cnv_ds1(i).quality1oxy(ii(jj(n))));
                if dm == '2'
                    cnv_ds1(i).quality1oxy(ii(jj(n)))= ...
                    cnv_ds1(i).quality1oxy(ii(jj(n)))+5;
                end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% JH qq %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%           
           end;
        end;
        
        clear pc pm qualtest dum dm
        [pc,pm]=cleaner(cnv_ds1(i).oxmmpri(ii),factor,order);
        
        jj=find(pm==0);
        for n=1:length(jj),
       
           if isfield(cnv_ds1,'oxmmpri'),
              qualtest=abs(cnv_ds1(i).oxmmpri(ii(jj(n)))-pc(jj(n)))>.0004;
           end;
           if qualtest,
              cnv_ds1(i).oxmmpri(ii(jj(n)))=pc(jj(n));
           end;
        end;
        
        clear pc pm qualtest dum dm
        [pc,pm]=cleaner(cnv_ds1(i).oxmllpri(ii),factor,order);
        
        jj=find(pm==0);
        for n=1:length(jj),
       
           if isfield(cnv_ds1,'oxmllpri'),
              qualtest=abs(cnv_ds1(i).oxmllpri(ii(jj(n)))-pc(jj(n)))>.0004;
           end;
           if qualtest,
              cnv_ds1(i).oxmllpri(ii(jj(n)))=pc(jj(n));
           end;
        end;
        
      end;% if length pressure
      clear ii jj dum dm
    end;% if 0 loop
   
   
   
% plot primary ctd oxygen voltage vs pressure 0-1000 db
 ii=find(cnv_ds1(i).pressure<=1000);
    %%%%%%%%%%% CMDR
    %fh=figure;
    fh=bigplot;%makes a larger plot window

    %buffer for plotting oxgyen axis 
    do = 0.03; 

    plot(cnv_ds1(i).pressure,cnv_ds1(i).(fldoxvpri),'g.-','linew',1.25);
    % we test if there are valid data for oxvpri
    if max(cnv_ds1(i).(fldoxvpri)(ii))-min(cnv_ds1(i).(fldoxvpri)(ii)) <= eps,
      axis([0 1000 min(cnv_ds1(i).(fldoxvpri)(ii))-do max(cnv_ds1(i).(fldoxvpri)(ii))+do]);
    else
      axis([0 1000 min(cnv_ds1(i).(fldoxvpri)(ii))-do max(cnv_ds1(i).(fldoxvpri)(ii))+do]);
    end;
    if (cnv_ds1(i).pressure<1000)
	if max(cnv_ds1(i).(fldoxvpri)(ii))-min(cnv_ds1(i).(fldoxvpri)(ii))<=eps,
            axis([0 max(cnv_ds1(i).pressure)+20  ...
		min(cnv_ds1(i).(fldoxvpri)(ii))-do max(cnv_ds1(i).(fldoxvpri)(ii))+do]);
        else
            axis([0 max(cnv_ds1(i).pressure)+20  ...
		min(cnv_ds1(i).(fldoxvpri)(ii))-do max(cnv_ds1(i).(fldoxvpri)(ii))+do]);
	end; 
    end
    title(['Station ' int2str(cnv_fill(i).station)  ...
		':  Primary Oxygen Voltage Despiking'])

% remove artifacts <=1000 db interactively; WOCE flag to 6

    loop1=menu('Interpolate?','yes','no');
    while loop1==1
        title(['Station ' int2str(cnv_fill(i).station)  ...
		':  Zoom in on current area of interest, type ' ...
		'''return'' when done']);
        zoom on
        keyboard
        zoom off
        title(['Station ' int2str(cnv_fill(i).station)  ...
		':  Click on end-points for interpolation']);
        [x,y]=ginput(2);
        x=round(x);
%%%%%% CMDR interp complains if the points are the same;
%%%%%% and dies if the order is reversed.
	if diff(x)<0, x=flipud(x); end;
	if diff(x)>1,
          kk=x(1)+1:x(2)-1;
          ll=[x(1),x(2)];
	if isfield( cnv_ds1,'oxvpri'),
          cnv_ds1(i).oxvpri(kk)=interp1(cnv_ds1(i).pressure(ll),  ...
		cnv_ds1(i).oxvpri(ll),cnv_ds1(i).pressure(kk));
    end;
    if isfield( cnv_ds1,'oxmllpri'),
          cnv_ds1(i).oxmllpri(kk)=interp1(cnv_ds1(i).pressure(ll),  ...
		cnv_ds1(i).oxmllpri(ll),cnv_ds1(i).pressure(kk));    
    end;
    if isfield( cnv_ds1,'oxmmpri'),
          cnv_ds1(i).oxmmpri(kk)=interp1(cnv_ds1(i).pressure(ll),  ...
		cnv_ds1(i).oxmmpri(ll),cnv_ds1(i).pressure(kk));    
    end;
	%if isfield(cnv_ds1,'oxuMkgpri'),
    %      cnv_ds1(i).oxuMkgpri(kk)=interp1(cnv_ds1(i).pressure(ll), ...
	%	cnv_ds1(i).oxuMkgpri(ll),cnv_ds1(i).pressure(kk));
	%end;
          for n=1:length(kk);
            dum=int2str(cnv_ds1(i).quality(kk(n)));
            if dum(6)=='2',
                cnv_ds1(i).quality(kk(n))=cnv_ds1(i).quality(kk(n))+40;
            end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% JH qq %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%            
            dm=int2str(cnv_ds1(i).quality1oxy(kk(n)));
            if dm == '2',
                cnv_ds1(i).quality1oxy(kk(n))=...
                    cnv_ds1(i).quality1oxy(kk(n))+5;
            end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% JH qq %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%            
          end;
%%%%%%%%%%% CMDR
figure(fh);
        plot(cnv_ds1(i).pressure,cnv_ds1(i).(fldoxvpri),'g.-','linew',1.25);
        axis([0 1000 min(cnv_ds1(i).(fldoxvpri)(ii))-do max(cnv_ds1(i).(fldoxvpri)(ii))+do]);
        if (cnv_ds1(i).pressure<2100)
            axis([0 max(cnv_ds1(i).pressure)+20 ...
		min(cnv_ds1(i).(fldoxvpri)(ii))-do max(cnv_ds1(i).(fldoxvpri)(ii))+do]);
        end;
else
	disp('Choose distinct x values');
end;
        title(['Station ' int2str(cnv_fill(i).station) ':  More Interpolating?'])
        loop1=menu('More Interpolating?','yes','no');     
    end 
    close
    clear ii kk ll dum dm do 
    
%%%######################################################################%%
% plot primary ctd oxygen voltage vs pressure 1000 - max pressure db

    ii=find(cnv_ds1(i).pressure>1000);
    if ~isempty(ii)
%%%%%%%%%%% CMDR
%fh=figure;
fh=bigplot;%makes a larger plot window
dod = 0.03;   
    plot(cnv_ds1(i).pressure,cnv_ds1(i).(fldoxvpri),'g.-','linew',1.25);
% we test if there are valid data for oxvpri
    if max(cnv_ds1(i).(fldoxvpri)(ii))-min(cnv_ds1(i).(fldoxvpri)(ii)) <= eps,
      axis([1000 max(cnv_ds1(i).pressure)+100 min(cnv_ds1(i).(fldoxvpri)(ii))-dod max(cnv_ds1(i).(fldoxvpri)(ii))+dod]);
    else
      axis([1000 max(cnv_ds1(i).pressure)+100 min(cnv_ds1(i).(fldoxvpri)(ii))-dod max(cnv_ds1(i).(fldoxvpri)(ii))+dod]);
    end;
    
    title(['Station ' int2str(cnv_fill(i).station)  ...
		':  Primary Oxygen Voltage Despiking'])

% remove artifacts <=1000 db interactively; WOCE flag to 6

    loop1=menu('Interpolate?','yes','no');
    while loop1==1
        title(['Station ' int2str(cnv_fill(i).station)  ...
		':  Zoom in on current area of interest, type ' ...
		'''return'' when done']);
        zoom on
        keyboard
        zoom off
        title(['Station ' int2str(cnv_fill(i).station)  ...
		':  Click on end-points for interpolation']);
        [x,y]=ginput(2);
        x=round(x);
%%%%%% CMDR interp complains if the points are the same;
%%%%%% and dies if the order is reversed.
	if diff(x)<0, x=flipud(x); end;
	if diff(x)>1,
          kk=x(1)+1:x(2)-1;
          ll=[x(1),x(2)];
	if isfield( cnv_ds1,'oxvpri'),
          cnv_ds1(i).oxvpri(kk)=interp1(cnv_ds1(i).pressure(ll),  ...
		cnv_ds1(i).oxvpri(ll),cnv_ds1(i).pressure(kk));    
    end;
    if isfield( cnv_ds1,'oxmllpri'),
          cnv_ds1(i).oxmllpri(kk)=interp1(cnv_ds1(i).pressure(ll),  ...
		cnv_ds1(i).oxmllpri(ll),cnv_ds1(i).pressure(kk));    
    end;
    if isfield( cnv_ds1,'oxmmpri'),
          cnv_ds1(i).oxmmpri(kk)=interp1(cnv_ds1(i).pressure(ll),  ...
		cnv_ds1(i).oxmmpri(ll),cnv_ds1(i).pressure(kk));    
    end;
	%if isfield(cnv_ds1,'oxuMkgpri'), 
    %      cnv_ds1(i).oxuMkgpri(kk)=interp1(cnv_ds1(i).pressure(ll), ...
	%	cnv_ds1(i).oxuMkgpri(ll),cnv_ds1(i).pressure(kk));
	%end;
          for n=1:length(kk);
            dum=int2str(cnv_ds1(i).quality(kk(n)));
            if dum(6)=='2',
                cnv_ds1(i).quality(kk(n))=cnv_ds1(i).quality(kk(n))+40;
            end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% JH qq %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%            
            dm=int2str(cnv_ds1(i).quality1oxy(kk(n)));
            if dm == '2',
                cnv_ds1(i).quality1oxy(kk(n))=cnv_ds1(i).quality1oxy(kk(n))+5;
            end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% JH qq %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%            
          end;
%%%%%%%%%%% CMDR
figure(fh);
        plot(cnv_ds1(i).pressure,cnv_ds1(i).(fldoxvpri),'g.-','linew',1.25);
        axis([1000 max(cnv_ds1(i).pressure)+100 min(cnv_ds1(i).(fldoxvpri)(ii))-dod max(cnv_ds1(i).(fldoxvpri)(ii))+dod]);
else
	disp('Choose distinct x values');
end;
        title(['Station ' int2str(cnv_fill(i).station) ':  More Interpolating?'])
        loop1=menu('More Interpolating?','yes','no');     
    end 
end;
    close 
    clear ii kk ll dum dm dod
%%%######################################################################%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%% SALINITY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% salinity despiking: plot full scale, interactively interpolate over spikes,
% and backcalculate cond, theta, sigma-theta.  repeat for deep water scale.
% Ammend salinity flag to 6.

% primary sal despiking: linearly interpolate over zeros; WOCE flag to 6
  % All good O2
  kk=find(cnv_ds1(i).salpri~=0 & cnv_ds1(i).salpri>=0);
	% sal==0 ;  all pr >= first good sal ; all pr <= last good sal
    % but only if there are good dalinity
    if ~isempty(kk),
      hh=find(abs(cnv_ds1(i).salpri==0 | cnv_ds1(i).salpri<0) | isnan(cnv_ds1(i).salpri)   & ...                             
          cnv_ds1(i).pressure >= cnv_ds1(i).pressure(kk(1)) & ...         
          cnv_ds1(i).pressure <= cnv_ds1(i).pressure(kk(length(kk))));
      if length(hh) > 0,
	if isfield(cnv_ds1,'salpri'),
          cnv_ds1(i).salpri(hh)=interp1(cnv_ds1(i).pressure(kk), ...
		  cnv_ds1(i).salpri(kk),cnv_ds1(i).pressure(hh));
	end;
    if isfield(cnv_ds1,'condpri'),
          cnv_ds1(i).condpri(hh)=interp1(cnv_ds1(i).pressure(kk), ...
		  cnv_ds1(i).condpri(kk),cnv_ds1(i).pressure(hh));
	end;
    
          for n=1:length(hh);
              dum=int2str(cnv_ds1(i).quality(hh(n)));
              if dum(6)=='2'
                  cnv_ds1(i).quality(hh(n))=cnv_ds1(i).quality(hh(n))+40;
              end; % if dum(6)=='2'
%%%%%%%%%%%%%%%%%%%%%%%%%%%% JH qq %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%              
              dm=int2str(cnv_ds1(i).quality1cond(hh(n)));
              if dm == '2'
                  cnv_ds1(i).quality1cond(hh(n))=...
                      cnv_ds1(i).quality1cond(hh(n))+4;
              end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% JH qq %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%              
          end; % for n=1:length(ii);
      end; % if length(ii) > 0,
    end; % if ~isempty(jj),
    clear kk hh dum dm
    
     if 0, % Not sure if I should use this on Salinity as well
% Secondary sal despiking: use cleaner program to despike deep water;
% WOCE flag to 6

    if length(cnv_ds1(i).pressure)>1000;
        ii=find(cnv_ds1(i).pressure>1000);
        factor = 1;
        order = 2;
        [pc,pm]=cleaner(cnv_ds1(i).salpri(ii),factor,order);

        jj=find(pm==0);
        for n=1:length(jj),
% Hmm. I wonder if trying to do this with dynmic fieldnames is such a good idea after all...
	if isfield(cnv_ds1,'salpri'),
	  qualtest=abs(cnv_ds1(i).salpri(ii(jj(n)))-pc(jj(n)))>0.0004; % not sure what to make this value????
	end;
        if qualtest,
                cnv_ds1(i).salpri(ii(jj(n)))=pc(jj(n));
                
%                 if isfield(cnv_ds1,'condpri')
%                     cnv_ds1(i).condpri(ii(jj(n)))=pc(jj(n));
%                 end;
                
                dum=int2str(cnv_ds1(i).quality(ii(jj(n))));
                if dum(6)=='2',
                    cnv_ds1(i).quality(ii(jj(n)))= ...
			cnv_ds1(i).quality(ii(jj(n)))+40;
                end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% JH qq %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%              
              dm=int2str(cnv_ds1(i).quality1cond(ii(jj(n))));
              if dm == '2'
                  cnv_ds1(i).quality1cond(ii(jj(n)))=...
                      cnv_ds1(i).quality1cond(ii(jj(n)))+5;
              end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% JH qq %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                
            end;
        end;
    end;
end;% if 0 loop
    

 % plot primary ctd salinity vs pressure 0-1000 db
 ii=find(cnv_ds1(i).pressure <= 1000);
    %%%%%%%%%%% CMDR
    %fh=figure;
    fh=bigplot;%makes a larger plot window
    %buffer for plotting the salinity axis
    ds = 0.05;
   
    plot(cnv_ds1(i).pressure, cnv_ds1(i).salpri,'r.-','linew',1.25);
    axis([0 1000 min(cnv_ds1(i).salpri(ii))-ds max(cnv_ds1(i).salpri(ii))+ds]);
    if (cnv_ds1(i).pressure<1000)
        axis([	0 max(cnv_ds1(i).pressure)+20  ...
		min(cnv_ds1(i).salpri(ii))-ds max(cnv_ds1(i).salpri(ii))+ds]);
    end
    title(['Station ' int2str(cnv_fill(i).station) ...
		 ':  Primary Salinity Despiking'])

% remove artifacts <=1000 db interactively; WOCE flag to 6

    loop1=menu('Interpolate?','yes','no');
    while loop1==1
        title(['Station ' int2str(cnv_fill(i).station) ...
		':  Zoom in on current area of interest, type ' ...
		'''return'' when done']);
        zoom on
        keyboard
        zoom off
        title(['Station ' int2str(cnv_fill(i).station) ...
		':  Click on end-points for interpolation']);
        [x,y]=ginput(2);
        x=round(x);
%%%%%% CMDR interp complains if the points are the same; 
%%%%%% and dies if the order is reversed.
	if diff(x)<0, x=flipud(x); end;
	if diff(x)>1,
          kk=x(1)+1:x(2)-1;
          ll=[x(1),x(2)];
          cnv_ds1(i).salpri(kk)=interp1(cnv_ds1(i).pressure(ll), ...
                cnv_ds1(i).salpri(ll),cnv_ds1(i).pressure(kk));
          %cnv_ds1(i).temppri(kk)=interp1(cnv_ds1(i).pressure(ll), ...
                %cnv_ds1(i).temppri(ll),cnv_ds1(i).pressure(kk));
          cnv_ds1(i).condpri(kk)=sw_c3515*sw_cndr(cnv_ds1(i).salpri(kk), ...
                cnv_ds1(i).temppri(kk),cnv_ds1(i).pressure(kk))/10;
          cnv_ds1(i).thetapri(kk)=sw_ptmp(cnv_ds1(i).salpri(kk), ...
                cnv_ds1(i).temppri(kk),cnv_ds1(i).pressure(kk),0);
          cnv_ds1(i).sigmapri(kk)=sw_pden(cnv_ds1(i).salpri(kk), ...
                cnv_ds1(i).thetapri(kk),cnv_ds1(i).pressure(kk),0)-1000;
          for n=1:length(kk),
            dum=int2str(cnv_ds1(i).quality(kk(n)));
            if dum(4)=='2',
                cnv_ds1(i).quality(kk(n))=cnv_ds1(i).quality(kk(n))+4000;
            end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% JH qq %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%              
              dm=int2str(cnv_ds1(i).quality1cond(kk(n)));
              if dm == '2'
                  cnv_ds1(i).quality1cond(kk(n))=...
                      cnv_ds1(i).quality1cond(kk(n))+5;
              end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% JH qq %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%            
          end;
%%%%%%%%%%% CMDR
figure(fh);
        plot(cnv_ds1(i).pressure, cnv_ds1(i).salpri, 'r.-','linew',1.25);
        axis([0 1000 min(cnv_ds1(i).salpri(ii))-ds max(cnv_ds1(i).salpri(ii))+ds]);
else
	disp('Choose distinct x values');
end;
        title(['Station ' int2str(cnv_fill(i).station)  ...
		':  More Interpolating?'])
        loop1=menu('More Interpolating?','yes','no');     
    end 
    close
    clear ii kk ll ds dum dm
   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot primary ctd salinity vs pressure 1000-6000 db

%buffer for plotting the deep salinity axis
dsd = 0.005;
    if length(cnv_ds1(i).pressure)>1000;
        ii=find(cnv_ds1(i).pressure>1000);
%%%%%%%%%%% CMDR
%fh=figure;
        fh=bigplot;%makes a larger plot window
        plot(cnv_ds1(i).pressure,cnv_ds1(i).salpri,'r.-','linew',1.25);
        axis([1000 6000 min(cnv_ds1(i).salpri(ii))-dsd max(cnv_ds1(i).salpri(ii))+dsd]);
        if (cnv_ds1(i).pressure<6000)
            axis([1000 max(cnv_ds1(i).pressure)+20  ...
		min(cnv_ds1(i).salpri(ii))-dsd max(cnv_ds1(i).salpri(ii))+dsd]);
        end
        title(['Station ' int2str(cnv_fill(i).station)  ...
		':  Primary Salinity Despiking - Deep'])

% remove artifacts >500 db interactively; WOCE flag to 6

        loop1=menu('Interpolate?','yes','no');
        while loop1==1
            title(['Station ' int2str(cnv_fill(i).station)  ...
		':  Zoom in on current area of interest, type ' ...
		'''return'' when done']);
            zoom on
            keyboard
            zoom off
            title(['Station ' int2str(cnv_fill(i).station)  ...
		':  Click on end-points for interpolation']);
            [x,y]=ginput(2);
            x=round(x);
%%%%%% CMDR interp complains if the points are the same;
%%%%%% and dies if the order is reversed.
	if diff(x)<0, x=flipud(x); end;
	if diff(x)>1,
            kk=x(1)+1:x(2)-1;
            ll=[x(1),x(2)];
         cnv_ds1(i).salpri(kk)=	...
            interp1(cnv_ds1(i).pressure(ll),	...
			cnv_ds1(i).salpri(ll),	...
			cnv_ds1(i).pressure(kk));
         %cnv_ds1(i).temppri(kk)=	...
            %interp1(cnv_ds1(i).pressure(ll),	...
			%cnv_ds1(i).temppri(ll),	...
			%cnv_ds1(i).pressure(kk));
         cnv_ds1(i).condpri(kk)=	...
            sw_c3515*sw_cndr(	...
			cnv_ds1(i).salpri(kk),	...
			cnv_ds1(i).temppri(kk),	...
			cnv_ds1(i).pressure(kk))/10;
         cnv_ds1(i).thetapri(kk)=	...
            sw_ptmp(cnv_ds1(i).salpri(kk),	...
			cnv_ds1(i).temppri(kk),	...
			cnv_ds1(i).pressure(kk),0);
         cnv_ds1(i).sigmapri(kk)=	...
            sw_pden(cnv_ds1(i).salpri(kk),	...
			cnv_ds1(i).thetapri(kk),	...
			cnv_ds1(i).pressure(kk),0)-1000;
            for n=1:length(kk);
                dum=int2str(cnv_ds1(i).quality(kk(n)));
                if dum(4)=='2'
                    cnv_ds1(i).quality(kk(n))=cnv_ds1(i).quality(kk(n))+4000;
                end
%%%%%%%%%%%%%%%%%%%%%%%%%%%% JH qq %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%              
              dm=int2str(cnv_ds1(i).quality1cond(kk(n)));
              if dm == '2'
                  cnv_ds1(i).quality1cond(kk(n))=...
                      cnv_ds1(i).quality1cond(kk(n))+5;
              end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% JH qq %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                
            end
%%%%%%%%%%% CMDR
figure(fh); 
            plot(cnv_ds1(i).pressure,cnv_ds1(i).salpri,'r.-','linew',1.25);
            axis([	500 ...
			6000 ...
			min(cnv_ds1(i).salpri(ii))-dsd ...
			max(cnv_ds1(i).salpri(ii))+dsd]);
            if (cnv_ds1(i).pressure<6000)
                axis([  0 ...
			max(cnv_ds1(i).pressure)+20 ...
			min(cnv_ds1(i).salpri(ii))-dsd ...
			max(cnv_ds1(i).salpri(ii))+dsd]);
            end
else
	disp('Choose distinct x values');
end;
            title(['Station ' int2str(cnv_fill(i).station) ...
		':  More Interpolating?'])
            loop1=menu('More Interpolating?','yes','no');     
        end 
    end
    close
%end %original end of station loop
    
end; %new end of station loop

%Save the primary deskipe summary structure.
save(db_file,'cnv_ds1','-append');

if nargout==0
	assignin('caller','ans',cnv_ds1);
else
	varargout{1}=cnv_ds1;
end;
return;


