function varargout=deloop(cruiseid,staa,stab)
% DELOOP -  similar to LOOPEDIT on Seabird Data Processing     
%
% CTD Calibration toolbox.
%
% USAGE
%	cnv_dlp=deloop(cruiseid,staa,stab);
%
% INPUT
%	cruiseid	cruise identification name
%	staa,stab	first and last station numbers to process
%   
% OUTPUT
%	cnv_dlp		data with pressure reversals removed 
%          
% DESCRIPTION
%	Remove potential density inversions created due to ship roll

% AUTHORS:Carlos Fonseca 
%         UM/CIMAS
%         Fri Dec 2 09:19:37 EDT 2005
% 
% VERSION: 
%	$Revision: 1.12 $
%	$Date: 2008/04/17 15:13:22 $
%	$Id: deloop.m,v 1.12 2008/04/17 15:13:22 duncombe Exp duncombe $
%
% CHANGELOG: 
%   23-Jul-2004, Version 1.0
%		* Initial version.
%
%  2007/12/12 CMDR
%  	* Get cruiseid from a file containing the identifier in the data
%  	directory.
% 	* put staa and stab in the input variables so that the stupid thing
%	doesn't keep asking for stations when you're doing testing.
%  2008/01/09
%	* if the oxygen concentration is not present, calculate it!
%  2008/02/17 CMDR
%	* comments and style corrections
%  2008/03/11 CMDR
%	* provide defaults for staa and stab
%  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/17 CMDR
%	* incorporate beckman oxygen sensors
%  2012 Feb M. Baringer
%   * cleaned up some indenting
%   * fixed b, b2 problem if b was not a column vector
%

%%%%%%%%%%% CMDR read cruiseid if we can
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

    % Load data from the cruise database
    db_file = fullfile(cruise_dir,[cruiseid '_db.mat']);
    load(db_file)
    if exist(db_file)==2
        g=who('cnv_ds1');
        if isempty(g)==0
            ld_file=(['load ',db_file,' cnv_ds1']);
            eval(ld_file)
            disp('despiked')            
        else
            ld_file = (['load ',db_file,' cnv_gterp']);
            eval(ld_file)
            disp('quick interp')
            cnv_ds1=cnv_gterp;
        end;
    else
        error(['The cruise data base file was not found on the path.  ' ...
		db_file]);
    end

    if isempty(g)==0,
        n=who('cnv_ds1');
        disp('despiked')        
    else
        n = who('cnv_gterp');
        disp('quick interp')
    end;

%%%%%%%%% CMDR set secondary sensors to primary data
% trying to run without secondary sensors
% comment out the ff code
%% if isempty(n)==1
%%  error(['cnv_ds2 not found. run despike2.m ']);
%%  return;
%%  end
% 
% Add the workaround:
%
%if isempty(n),	
%    disp(['No secondary sensors. Winging it with primary ' ...
%			'sensors only.  \nIf this is not what you expect, perhaps run despike2.m']);		
%    eval(['load ',db_file,' cnv_ds1']);
%    cnv_ds2=cnv_ds1;
	if isempty(cnv_ds1),
		error(['Fatal error: no cnv_ds1 data']);
	end;
%end;
%%%%%%%%% 
seabird=false;
beckman=false;
if isfield(cnv_ds1,'oxcpri'),
	beckman=true;
elseif isfield(cnv_ds1,'oxvpri'),
	seabird=true;
end;

% 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_ds1.station));

disp('****************************************** ')
disp('DELOOP STATION')
disp('****************************************** ')
disp('The data found in cnv goes from') 
disp(['the station ',num2str(min(aux_stations)),' to the station ', ...
	num2str(max(aux_stations)),'.'] ) 
disp('****************************************** ')
%
%%%%%%%%%%%% CMDR enable putting station limits in function call
% staa=input('Input first station: ');
% stab=input('Input last station: ');
if ~exist('staa'), staa=input('Input first station: '); end;
if ~exist('stab'), stab=input('Input last station: '); end;
if isempty(staa), staa=min([cnv_ds1.station]); end;
if isempty(stab), stab=max([cnv_ds1.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
 cnv_dlp = cnv_ds1;

%initial setup
% and loop through the known order of ind_stns
for i = ind_stns(sta_min_idx:sta_max_idx).',
    % Compute N^2 which is returned as b
    [b,dum,pav]=sw_bfrq(cnv_dlp(i).salpri,cnv_dlp(i).temppri,cnv_dlp(i).pressure);
    % make sure b is a column vector
    b = b(:);
    
    % Create b2 which is composed of two vectors: one is b padded 
    % with a zero at the surface and one is b padded with a zero
    % at depth (zero is above the N^2 threshhold that follows).

    b2=[[0;b],[b;0]];

    % This brackets b such that it can be tested and the first
    % and second points of a -N^2 failing the criteria (<=-1e-5) 
    % interpolated over.  Keep the filled surface data in spite
    % of any inversions.
 
    ii=find((b2(:,1)<=-1e-5|b2(:,2)<=-1e-5)&cnv_dlp(i).quality<2227777);
    jj=find((b2(:,1)>-1e-5&b2(:,2)>-1e-5)|cnv_dlp(i).quality>=2227777);

    % Linearly interpolate over inversions

    cnv_dlp(i).temppri(ii)=		...
        interp1(cnv_dlp(i).pressure(jj),		...
            cnv_dlp(i).temppri(jj),		...
            cnv_dlp(i).pressure(ii), 'linear', 'extrap');
    cnv_dlp(i).condpri(ii)=		...
        interp1(cnv_dlp(i).pressure(jj),		...
            cnv_dlp(i).condpri(jj),		...
            cnv_dlp(i).pressure(ii), 'linear', 'extrap');
    if seabird,
        cnv_dlp(i).oxvpri(ii)=		...
            interp1(cnv_dlp(i).pressure(jj),		...
            cnv_dlp(i).oxvpri(jj),		...
            cnv_dlp(i).pressure(ii), 'linear', 'extrap');
    elseif beckman,
        cnv_dlp(i).oxcpri(ii)=			...
            interp1(cnv_dlp(i).pressure(jj),		...
            cnv_dlp(i).oxcpri(jj),		...
            cnv_dlp(i).pressure(ii),'linear','extrap');
        cnv_dlp(i).oxctpri(ii)=		...
            interp1(cnv_dlp(i).pressure(jj),		...
            cnv_dlp(i).oxctpri(jj),		...
            cnv_dlp(i).pressure(ii),'linear','extrap');
    end;
    %
    %%%%%%%%%% CMDR check if oxygen data are available
    %%%%%%%%%% and calculate the oxygen concentration if it is not 
    if isempty(cnv_dlp(i).oxmllpri),
    	disp(['cnv_dlp(' num2str(i) ').oxmllpri is empty.']);
    	disp(['Oxygen concentration will be ' ...
    		'calculated using default sensor constants.']);
    	if seabird,
    	   cnv_dlp(i).oxmllpri=oxy_calc_none( ...
    		cnv_dlp(i).oxvpri, ...
    		cnv_dlp(i).salpri, ...
    		cnv_dlp(i).temppri, ...
    		cnv_dlp(i).pressure);
    	elseif beckman,
    	   error(['You will have to go back and create a module to ' ...
    		'calculate oxygen for the beckman sensor analogous ' ...
    		'to oxy_calc_none.m']);
    	   return;
        end;
    else
        cnv_dlp(i).oxmllpri(ii)=  ...
            interp1(cnv_dlp(i).pressure(jj),  ...
			cnv_dlp(i).oxmllpri(jj),  ...
			cnv_dlp(i).pressure(ii),'linear','extrap');
    end;
    %%%%%%%%%% 
    %
    % Compute N^2 which is returned as b
    %%%%%%%%%%%%% CMDR if secondary salinity not available use primary for bfq
    if isempty(cnv_dlp(i).salsec), 
        [b,dum,pav]=sw_bfrq(cnv_dlp(i).salpri,cnv_dlp(i).temppri,cnv_dlp(i).pressure);
    else
        [b,dum,pav]=sw_bfrq(cnv_dlp(i).salsec,cnv_dlp(i).tempsec,cnv_dlp(i).pressure);
    end;
    b = b(:);
    %%%%%%%%%%%%%

    % Create b2 which is composed of two vectors: one is b padded 
    % with a zero at the surface and one is b padded with a zero
    % at depth (zero is above the N^2 threshhold that follows).

    b2=[[0;b],[b;0]];

    % This brackets b such that it can be tested and the first
    % and second points of a -N^2 failing the criteria (<=-1e-5) 
    % interpolated over.  Keep the filled surface data in spite
    % of any inversions.
 
    ii=find((b2(:,1)<=-1e-5|b2(:,2)<=-1e-5)&cnv_dlp(i).quality<2227777);
    jj=find((b2(:,1)>-1e-5&b2(:,2)>-1e-5)|cnv_dlp(i).quality>=2227777);

    % Linearly interpolate over inversions
    %
    %%%%%%%%%%%%%%%%%%%%%%%% CMDR skip if secondary sensors not available
    if ~isempty(cnv_dlp(i).tempsec), 
        cnv_dlp(i).tempsec(ii)=interp1(cnv_dlp(i).pressure(jj), ...
				cnv_dlp(i).tempsec(jj), ...
				cnv_dlp(i).pressure(ii),'linear','extrap');
        cnv_dlp(i).condsec(ii)=interp1(cnv_dlp(i).pressure(jj), ...
				cnv_dlp(i).condsec(jj), ...
				cnv_dlp(i).pressure(ii),'linear','extrap');
    end  
    if seabird,
        jj=find(cnv_dlp(i).oxvsec~=0); 
    elseif  beckman,
        jj=find(cnv_dlp(i).oxcsec~=0);
    end
    if ~isempty(jj) 
      if seabird,
          cnv_dlp(i).oxvsec(ii)=interp1(cnv_dlp(i).pressure(jj), ...
				cnv_dlp(i).oxvsec(jj), ...
				cnv_dlp(i).pressure(ii),'linear','extrap');
      elseif beckman,
          cnv_dlp(i).oxcsec(ii)=interp1(cnv_dlp(i).pressure(jj), ...
				cnv_dlp(i).oxcsec(jj), ...
				cnv_dlp(i).pressure(ii),'linear','extrap');
          cnv_dlp(i).oxctsec(ii)=interp1(cnv_dlp(i).pressure(jj), ...
				cnv_dlp(i).oxctsec(jj), ...
				cnv_dlp(i).pressure(ii),'linear','extrap');
      end;
      cnv_dlp(i).oxmllsec(ii)=interp1(cnv_dlp(i).pressure(jj), ...
				cnv_dlp(i).oxmllsec(jj), ...
				cnv_dlp(i).pressure(ii),'linear','extrap');
    end;
    %%%%%%%%%%%%%%%%%%%%%%%% 

    % Set WOCE quality word to '2666' over interpolated records

    for n=1:length(ii);
        dum=int2str(cnv_dlp(i).quality(ii(n)));
        if dum(2)=='2'
            cnv_dlp(i).quality(ii(n))=cnv_dlp(i).quality(ii(n))+400000;
        end
        if dum(3)=='2'
            cnv_dlp(i).quality(ii(n))=cnv_dlp(i).quality(ii(n))+40000;
        end
        if dum(4)=='2'
            cnv_dlp(i).quality(ii(n))=cnv_dlp(i).quality(ii(n))+4000;
        end
        if dum(5)=='2'
            cnv_dlp(i).quality(ii(n))=cnv_dlp(i).quality(ii(n))+400;
        end
        if dum(6)=='2'
            cnv_dlp(i).quality(ii(n))=cnv_dlp(i).quality(ii(n))+40;
        end
        if dum(7)=='2'
            cnv_dlp(i).quality(ii(n))=cnv_dlp(i).quality(ii(n))+4;
        end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% JH qq %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
        d1ox = int2str(cnv_dlp(i).quality1oxy(ii(n)));
            if d1ox == '2'
                cnv_dlp(i).quality1oxy(ii(n)) =...
                    cnv_dlp(i).quality1oxy(ii(n)) + 4;
            end;
        d2ox = int2str(cnv_dlp(i).quality2oxy(ii(n)));
            if d2ox == '2'
                cnv_dlp(i).quality2oxy(ii(n)) =...
                    cnv_dlp(i).quality2oxy(ii(n)) + 4;
            end;
        d1co = int2str(cnv_dlp(i).quality1cond(ii(n)));
            if d1co == '2'
                cnv_dlp(i).quality1cond(ii(n)) =...
                    cnv_dlp(i).quality1cond(ii(n)) + 4;
            end;
        d2co = int2str(cnv_dlp(i).quality2cond(ii(n)));
            if d2co == '2'
                cnv_dlp(i).quality2cond(ii(n)) =...
                    cnv_dlp(i).quality2cond(ii(n)) + 4;
            end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% JH qq %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
    end

    % Recompute derived variables
    cnv_dlp(i).salpri(ii)=sw_salt(cnv_dlp(i).condpri(ii)*10/sw_c3515,cnv_dlp(i).temppri(ii),cnv_dlp(i).pressure(ii));
    cnv_dlp(i).thetapri(ii)=sw_ptmp(cnv_dlp(i).salpri(ii),cnv_dlp(i).temppri(ii),cnv_dlp(i).pressure(ii),0);
    cnv_dlp(i).sigmapri(ii)=sw_pden(cnv_dlp(i).salpri(ii),cnv_dlp(i).thetapri(ii),cnv_dlp(i).pressure(ii),0)-1000;

    %%%%%%%%%%% CMDR skip if secondary sensors not available
    if ~isempty(cnv_dlp(i).condsec),
        cnv_dlp(i).salsec(ii)=sw_salt(  ...
            cnv_dlp(i).condsec(ii)*10/sw_c3515,cnv_dlp(i).tempsec(ii),  ...
            cnv_dlp(i).pressure(ii));
        cnv_dlp(i).thetasec(ii)=sw_ptmp(  ...
            cnv_dlp(i).salsec(ii),cnv_dlp(i).tempsec(ii),  ...
            cnv_dlp(i).pressure(ii),0);
        cnv_dlp(i).sigmasec(ii)=sw_pden(  ...
            cnv_dlp(i).salsec(ii),cnv_dlp(i).thetasec(ii),  ...
            cnv_dlp(i).pressure(ii),0)-1000;
    end;
    %%%%%%%%%%% 

    disp(['delooped station ' num2str(cnv_dlp(i).station) ]);
    % setting up the correction variables
end;

%Save the conductivity summary structure.
save(db_file,'cnv_dlp','-append');

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

return;

