function varargout=fill_surface_beckman(cruiseid,staa,stab)
% FILL_SURFACE_BECKMAN -  fill the CTD data up to the surface
%    
% CTD Calibration toolbox.
%
% USAGE
%	cnv_fill=fill_surface(cruiseid,staa,stab);
%
% INPUT
%	cruiseid	cruise identification name.
%	staa,stab	first and last stations to process
%	cnv_load	from database file
%   
% OUTPUT
% 	cnv_fill	data extrapolated to the surface
%          
% DESCRIPTION: This function uses the cnv files from the cruise database
% and looks for the first valid data for salinity, oxygen and temperature and
% copies this value up to the surface (e.g. sets a surface "mixed layer").
%

% AUTHORS:Carlos Fonseca 
%         UM/CIMAS
%         Fri Dec 2 09:19:37 EDT 2005
%
% VERSION: 
%	$Revision: 1.5 $
%	$Date: 2008/04/16 17:28:51 $
%	$Id: fill_surface_beckman.m,v 1.5 2008/04/16 17:28:51 duncombe Exp duncombe $
% 
% CHANGELOG: 
%   23-Jul-2004, Version 1.0
%	* Initial version.
%   24-Mar-2006 Molly Baringer
%       * Major overhaul of this routine.  Should probably be called version 2.
%       * Cleaned up comments
%       * Added instructions for plots etc.
%       * uses cnv_load instead of parsing cnv directly.
%       * sets scales for plots based on actual values indendent of scan,
%       bin averaging etc.
%       * algorithm for padding mixed layer works for both down and upcast
%       profiles (not tested for both, but in principle that should work too!).
%       * log file kept that includes all modifications made to casts.
%       * sigma originally computed from Theta - should be from insitu Te.
%       * fixed problem with having to pad *all* station *every time*.
%  2008/04/10 CMDR
% 	* incorporate changes made in fill_surface.m
%

% PROGRAMMING NOTES:
%	* should be simple to change fill_surface.m to handle the Beckman
%	data (test for oxyv or oxyc and oxyct and process identically) 
%	and declare this module redundant.

% Set depth range for plot
% Used for determining the mean values for the axes of the plots.
% Mean values for the plots are taken about the mean values of these depth
max_depth1 = 18.0;
max_depth2 = 22.0;

%%%%%%%%%% 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_load');
if isempty(n)==1
    error(['cnv_load not found. Run load_cnv_molly ']);
    return;
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_load.station));

disp('****************************************** ')
disp('FILL SURFACE DATA')
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
% 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(aux_stations); end;
if isempty(stab), stab=max(aux_stations); 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(aux_stations==staa);
sta_max_idx = find(aux_stations==stab);       

% and loop through the known order of ind_stns
for i = ind_stns(sta_min_idx:sta_max_idx).',
    display(['Processing station ' int2str(cnv_load(i).station) ' index ' ...
	 int2str(i)])
    cnv_fill(i) = cnv_load(i);
    II = find(cnv_fill(i).pressure<=max_depth2 & ...
	 cnv_fill(i).pressure>=max_depth1);

    if (length(II) <= 0)
        display(['  FILL SURFACE WARNING: No pressure values near ', ...
		'20 dbar for Station ' int2str(cnv_fill(i).station)])
        display(['  Moving on to next station number.  Check cnv data '...
		'file manually.'])
    else
        
%%%%%%%%%%% CMDR
%%        nf=fopen([group_pref.cal_hdr_dir '\' group_pref.file_name_prefix ...
%% 		num2str(cnv_fill(i).station,'%3.3d') '.log'],'wt');
%%%%% heh-heh aah this might not work on a DOS machine which  may return "C" 
%%%%% as the first char of pwd
	p=pwd; dir_del=p(1);
        nf=fopen([group_pref.cal_hdr_dir dir_del ...
	group_pref.file_name_prefix ...
	num2str(cnv_fill(i).station,'%3.3d') '.log'],'wt');
%%%%%%%%%%%

% plot primary and secondary salinity profiles 0-20 db
%%%%%%%%%%%% CMDR do not overwrite old figures
%        figure(1);
	fig1=figure;
%%%%%%%%%%%%
%%%%%%%%%%%% CMDR
            plot(cnv_fill(i).pressure,cnv_fill(i).salpri,'b.-');
%%%%%%%%%%%% 
%           set the axes based on the value at 20 m
%% this turns out to be a really bad idea if there is great 
%% variability in the surface layer
        yaxis_mean_sal = meanmiss( ...
	cnv_fill(i).salpri(find(cnv_fill(i).pressure<=max_depth2 & ...
	cnv_fill(i).pressure>=max_depth1)));
            axis([0 20 yaxis_mean_sal-0.5 yaxis_mean_sal+0.5]);
            xlabel('Pressure');ylabel('Salinity');
        title(['Station ' int2str(cnv_fill(i).station) ...
		':  Click on first good depth to pad surface layer']);
            grid on;
            hold on;
%%%%%%%%%%%% CMDR
            if ~isempty(cnv_fill(i).salsec), 
		cnv_fill(i).salsec
		plot(cnv_fill(i).pressure,cnv_fill(i).salsec,'r.-'); 
	    end;
%%%%%%%%%%%% 
            hold off;

% Select mixed layer values
% click on db of first good salinity data record
        figure(fig1);
            [sx,y]=ginput(1); 
            sx = round(sx);
% sx now contains the (closest) deepest pressure level the user has
% determined to be a good value.  y is discarded for a more precise
% calculation.  

% copy salinities back to surface
% II is index of "mixed layer" values
% Maximum pressure II value is used for padding
% Note:  unlike the earlier algorithm - this padding will work for upcasts
% *and* downcasts.  It should actually work for casts with up *and*
% downcasts, but that's not really what we want to be dealing with.
%%%%%%%%%%% CMDR
%%%%%%%%%%% take secondary sensors away from the primary sensors.
%%%%%%%%%%% We'll be able to have one test to see if the data exist or not.
%
% First pad salinity then recompute all the extra related parameters
% recompute theta (which is a function of Salinity and is needed for 
% the next step) 
% recompute conductivity (so that conductivity would give the same S)
% recompute sigma
        II_sal = find (cnv_fill(i).pressure <= sx);
    
        cnv_fill(i).salpri(II_sal) = cnv_fill(i).salpri(...
	   find(cnv_fill(i).pressure ==max(cnv_fill(i).pressure(II_sal))));
        cnv_fill(i).thetapri(II_sal)=sw_ptmp(cnv_fill(i).salpri(II_sal), ...
		cnv_fill(i).temppri(II_sal), cnv_fill(i).pressure(II_sal),0);
        cnv_fill(i).condpri(II_sal)=sw_c3515* ...
		sw_cndr(cnv_fill(i).salpri(II_sal), ...
		cnv_fill(i).temppri(II_sal),cnv_fill(i).pressure(II_sal))/10;
        cnv_fill(i).sigmapri(II_sal)=sw_pden(cnv_fill(i).salpri(II_sal), ...
		cnv_fill(i).temppri(II_sal), ...
		cnv_fill(i).pressure(II_sal),0)-1000.;

	if ~isempty(cnv_fill(i).salsec), 
        cnv_fill(i).salsec(II_sal) = cnv_fill(i).salsec(...
	    find(cnv_fill(i).pressure ==max(cnv_fill(i).pressure(II_sal))));
        cnv_fill(i).thetasec(II_sal)=sw_ptmp(cnv_fill(i).salsec(II_sal), ...
	    cnv_fill(i).tempsec(II_sal), cnv_fill(i).pressure(II_sal),0);
        cnv_fill(i).condsec(II_sal)=sw_c3515*sw_cndr(...
	    cnv_fill(i).salsec(II_sal),cnv_fill(i).tempsec(II_sal),...
	    cnv_fill(i).pressure(II_sal))/10;
        cnv_fill(i).sigmasec(II_sal)=sw_pden(cnv_fill(i).salsec(II_sal), ...
	    cnv_fill(i).tempsec(II_sal), cnv_fill(i).pressure(II_sal),0)-1000.;
	end;
%%%%%%%%%%%%%%%%
 
% plot primary and secondary theta for comparison, similar scale
%%%%%%%%%%%% CMDR do not overwrite old figures
%        figure(2);
	fig2=figure;
%%%%%%%%%%%%
%%%%%%%%%%%% CMDR
            plot(cnv_fill(i).pressure,cnv_fill(i).thetapri,'y.-');
% this does not work with a non-mixed surface layer! 
%	    idx=find(cnv_fill(i).pressure<=max_depth2 & ...
%		cnv_fill(i).pressure>=max_depth1);
%	    yaxis_mean_th = meanmiss(cnv_fill(i).thetapri(idx));
%           axis([0 20 ...
%	yaxis_mean_th-0.5/sw_aonb(yaxis_mean_sal,yaxis_mean_th,15) ...
%	yaxis_mean_th+0.5/sw_aonb(yaxis_mean_sal,yaxis_mean_th,15)]);
%% replace it with
	    idx=find(cnv_fill(i).pressure<=max_depth2);
	axis([0 20 ...
		round(min(cnv_fill(i).thetapri(idx)))-0.5 ...
		round(max(cnv_fill(i).thetapri(idx)))+0.5]);
            xlabel('Pressure');ylabel('Theta');
            title(['Station ' int2str(cnv_fill(i).station) ...
	':  Click on first good depth to pad surface layer'])
           grid on;
            hold on;
%%%%%%%%%%%% CMDR
            if ~isempty(cnv_fill(i).thetasec),
		 plot(cnv_fill(i).pressure,cnv_fill(i).thetasec,'c.-');
	    end;
%%%%%%%%%%%%
            hold off;

% click on db of first good theta data record

%%%%%%%%%%%% CMDR
%   	figure(2)
        figure(fig2);
            [tx,y]=ginput(1); 
        tx = round(tx);

% tx now contains the deepest pressure level the user has determined to be
% a good value.  y is discarded for a more precise calculation.

% copy theta back to surface
% II is index of "mixed layer" values
% Maximum pressure II value is used for padding
% Note:  unlike the earlier algorithm - this padding will work for upcasts
% *and* downcasts.  It should actually work for casts with up *and*
% downcasts, but that's not really what we want to be dealing with.
%%%%%%%%%%% CMDR
%%%%%%%%%%% take secondary sensors away from the primary sensors.
%%%%%%%%%%% We'll be able to have one test to see if the data exist or not.
% Now have pad theta and recompute all the extra related parameters
% recompute temperature (from theta where Theta was padded)
% recompute conductivity (so that given the new T, the same S will be
% computed consistently)
% recompute sigma

        II_te = find (cnv_fill(i).pressure <= tx);
        cnv_fill(i).thetapri(II_te) = cnv_fill(i).thetapri( ...
		find(cnv_fill(i).pressure ==max(cnv_fill(i).pressure(II_te))));
        cnv_fill(i).temppri(II_te)=sw_temp(cnv_fill(i).salpri(II_te), ...
		 cnv_fill(i).thetapri(II_te), cnv_fill(i).pressure(II_te),0);
        cnv_fill(i).condpri(II_te)=sw_c3515*sw_cndr(...
		cnv_fill(i).salpri(II_te),cnv_fill(i).temppri(II_te),...
		cnv_fill(i).pressure(II_te))/10;
        cnv_fill(i).sigmapri(II_te)=sw_pden(cnv_fill(i).salpri(II_te), ...
		cnv_fill(i).temppri(II_te), ...
		 cnv_fill(i).pressure(II_te),0)-1000.;

	if ~isempty(cnv_fill(i).thetasec),
        cnv_fill(i).thetasec(II_te) = cnv_fill(i).thetasec( ...
		find(cnv_fill(i).pressure ==max(cnv_fill(i).pressure(II_te))));
        cnv_fill(i).tempsec(II_te)=sw_temp(cnv_fill(i).salsec(II_te), ...
		 cnv_fill(i).thetasec(II_te), ...
		 cnv_fill(i).pressure(II_te),0);
        cnv_fill(i).condsec(II_te)=sw_c3515*sw_cndr( ...
		cnv_fill(i).salsec(II_te),cnv_fill(i).tempsec(II_te), ...
		cnv_fill(i).pressure(II_te))/10;
        cnv_fill(i).sigmasec(II_te)=sw_pden(cnv_fill(i).salsec(II_te), ...
		cnv_fill(i).tempsec(II_te),  ...
		cnv_fill(i).pressure(II_te),0)-1000.;
	end;
%%%%%%%%%%%%%%%%%

% plot primary and secondary ctd oxygen
%%%%%%%%%%% CMDR
%        figure(3);
        fig3=figure;
%%%%%%%%%%% CMDR
            plot(cnv_fill(i).pressure,cnv_fill(i).oxcpri,'r.-');
            yaxis_mean_ox = meanmiss(cnv_fill(i).oxcpri( ...
		find(cnv_fill(i).pressure<=max_depth2 &  ...
		cnv_fill(i).pressure>=max_depth1)));
            axis([0 20 yaxis_mean_ox-1 yaxis_mean_ox+1]);
            xlabel('Pressure');ylabel('Oxygen');
            title(['Station ' int2str(cnv_fill(i).station)  ...
		':  Click on first good depth to pad surface layer'])
            grid on;
            hold on;
%%%%%%%%%% CMDR
            if ~isempty(cnv_fill(i).oxcsec), 
            	plot(cnv_fill(i).pressure,cnv_fill(i).oxcsec,'g.-');
	    end;
            hold off;

% click on db of first good oxygen data record
%%%%%%%%%%%% CMDR
%       figure(3);
        figure(fig3);
            [ux,y]=ginput(1); 
        ux = round(ux);

% copy values back to surface
% II is index of "mixed layer" values
% Maximum pressure II value is used for padding
% Note:  unlike the earlier algorithm - this padding will work for upcasts
% *and* downcasts.  It should actually work for casts with up *and*
% downcasts, but that's not really what we want to be dealing with.
        II_ox = find (cnv_fill(i).pressure <= ux);
%%%%%%%%%CMDR

        if ~isempty(cnv_fill(i).oxmllpri),
        	cnv_fill(i).oxmllpri(II_ox) =  ...
			cnv_fill(i).oxmllpri(...
			find(cnv_fill(i).pressure ==max(...
			cnv_fill(i).pressure(II_ox))));
	end;

        if ~isempty(cnv_fill(i).oxcpri),
            cnv_fill(i).oxcpri(II_ox) = cnv_fill(i).oxcpri(...
	    find(cnv_fill(i).pressure ==max(cnv_fill(i).pressure(II_ox))));
	end;
        if ~isempty(cnv_fill(i).oxctpri),
            cnv_fill(i).oxctpri(II_ox) = cnv_fill(i).oxctpri(...
	    find(cnv_fill(i).pressure ==max(cnv_fill(i).pressure(II_ox))));
	end;

        if ~isempty(cnv_fill(i).oxmllsec),
        	cnv_fill(i).oxmllsec(II_ox) = cnv_fill(i).oxmllsec(...
			find(cnv_fill(i).pressure ==max( ...
			cnv_fill(i).pressure(II_ox))));
	end;

        if ~isempty(cnv_fill(i).oxcsec),
        	cnv_fill(i).oxcsec(II_ox) = cnv_fill(i).oxcsec(...
			find(cnv_fill(i).pressure ==max(...
			cnv_fill(i).pressure(II_ox))));
	end;
        if ~isempty(cnv_fill(i).oxctsec),
        	cnv_fill(i).oxctsec(II_ox) = cnv_fill(i).oxctsec(...
			find(cnv_fill(i).pressure ==max(...
			cnv_fill(i).pressure(II_ox))));
	end;

%%%%%%%%%%%
  
% Need to adjust quality flags for the values replaced
% set up iql=2222222+pres*0;
% s   iql(j)=iql(j)+5005500;
% t  iql(j)=iql(j)+550000;
% o2   iql(j)=iql(j)+55;
% write out log file for all station modifications
        cnv_fill(i).quality(union(union(II_sal,II_te),II_ox)) =     2222222;
        cnv_fill(i).quality(II_sal) = cnv_fill(i).quality(II_sal) + 5005500;
        cnv_fill(i).quality(II_te) = cnv_fill(i).quality(II_te) +    550000;
        cnv_fill(i).quality(II_ox) = cnv_fill(i).quality(II_ox) +        55;

% write out log file that includes adjustments for this station
        fprintf(nf, ['PAD TO SURFACE at Station ' ...
		int2str(cnv_fill(i).station) ' with:\n']);
        fprintf(nf, ['  Salinity = ' ...
		num2str(cnv_fill(i).salpri(...
		find(cnv_fill(i).pressure ==max(...
		cnv_fill(i).pressure(II_sal))))) ...
		' from pr = ' ...
		num2str(max(cnv_fill(i).pressure(II_sal))) '\n']);
        fprintf(nf, ['  Theta    = ' num2str(cnv_fill(i).thetapri(...
		find(cnv_fill(i).pressure ==max(...
		cnv_fill(i).pressure(II_te))))) ' from pr = ' ...
		num2str(max(cnv_fill(i).pressure(II_te))) '\n']);
%%%%%%%% CMDR
        if ~isempty(cnv_fill(i).oxmllpri),
        	fprintf(nf, ['  Oxygen   = ' ...
		num2str(cnv_fill(i).oxmllpri(...
		find(cnv_fill(i).pressure ==max(...
		cnv_fill(i).pressure(II_ox))))) ' from pr = ' ...
		num2str(max(cnv_fill(i).pressure(II_ox))) '\n']);
	end;
%%%%%%% 

        fprintf(nf, ['PAD TO SURFACE at Station ' ...
		int2str(cnv_fill(i).station) ' with:\n']);
        fprintf(nf, ['  Salinity = ' ...
		num2str(cnv_fill(i).salpri(find(...
		cnv_fill(i).pressure ==max(cnv_fill(i).pressure(II_sal))))) ...
		' from pr = ' num2str(max(cnv_fill(i).pressure(II_sal))) '\n']);
        fprintf(nf, ['  Theta    = ' num2str(cnv_fill(i).thetapri(...
		find(cnv_fill(i).pressure ==max(...
		cnv_fill(i).pressure(II_te))))) ' from pr = ' ...
		num2str(max(cnv_fill(i).pressure(II_te))) '\n']);
%%%%%%%% CMDR
        if ~isempty(cnv_fill(i).oxmllpri),
        	fprintf(nf, ['  Oxygen   = ' num2str(...
		cnv_fill(i).oxmllpri(...
		find(cnv_fill(i).pressure ==max(...
		cnv_fill(i).pressure(II_ox))))) ' from pr = ' ...
		num2str(max(cnv_fill(i).pressure(II_ox))) '\n']);
	end;
%%%%%%% 
%%
        close(fig1); close(fig2); close(fig3);
    end         % there exist pressure values near 20 dbar
end             % proceed to next station

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

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

return;
  
