function varargout=fill_surface(cruiseid)
% FILL SURFACE -  Fill the CTD data up to the surface
%    
% CTD Calibration toolbox.
%
% USAGE:
% cnv_fill=fill_surface(cruiseid);
%
%INPUT: cruiseid: Cruise identification name.
%   
% OUTPUT: .  
%          
% 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
% 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*.
%   Feb 2012 M. Baringer
%       * changed to automatically fill missing surface values with last
%       "good" value automatically.
%   Feb 2012 J. Hooper
%       * changed cnv_load(1) -> cnv_load(i) in the II = find();  The
%       fill surface routine was only working on the first station and
%       leaving NaN's in the rest!

% 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 ofr these depth
max_depth1 = 18.0;
max_depth2 = 22.0;

% 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

aux_stations=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('****************************************** ')
%
staa=input('Input first station: ');
stab=input('Input last station: ');

% 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.
sta_min_idx = find([cnv_load.station]'==staa);
sta_max_idx = find([cnv_load.station]'==stab);       

for i = sta_min_idx: sta_max_idx
    cnv_fill(i) = cnv_load(i);
        
    nf=fopen([group_pref.calibration_dir '/log' '/' group_pref.file_name_prefix num2str(cnv_fill(i).station,'%3.3d') '.log'],'wt');


% 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.
    %JH - find surface NaN values only
   
    II = find(isnan(cnv_load(i).temppri) & cnv_load(i).pressure <=20);%JH changed cnv_load(1)->cnv_load(i) replaced '1' with the indice 'i'
   
    cnv_fill(i).temppri(II) = cnv_fill(i).temppri(max(II)+1);
    cnv_fill(i).tempsec(II) = cnv_fill(i).tempsec(max(II)+1);

    cnv_fill(i).salpri(II) = cnv_fill(i).salpri(max(II)+1);
    cnv_fill(i).salsec(II) = cnv_fill(i).salsec(max(II)+1);

    cnv_fill(i).oxmllpri(II) = cnv_fill(i).oxmllpri(max(II)+1);
    cnv_fill(i).oxmllsec(II) = cnv_fill(i).oxmllsec(max(II)+1);     

    cnv_fill(i).oxvpri(II) = cnv_fill(i).oxvpri(max(II)+1);
    cnv_fill(i).oxvsec(II) = cnv_fill(i).oxvsec(max(II)+1);
  
    if isfield(cnv_fill,'oxmmpri'),
        cnv_fill(i).oxmmpri(II) = cnv_fill(i).oxmmpri(max(II)+1);
    end;
    
    if isfield(cnv_fill,'oxmmsec'),
        cnv_fill(i).oxmmsec(II) = cnv_fill(i).oxmmsec(max(II)+1);
    end;
   

    
% Now you have padded salinity and temperature
% You now have to recompute all the extra related parameters
% recompute theta (which is a function of Salinity and is needed for the next step) 
    cnv_fill(i).thetapri(II)=sw_ptmp(cnv_fill(i).salpri(II), cnv_fill(i).temppri(II), cnv_fill(i).pressure(II),0);
    cnv_fill(i).thetasec(II)=sw_ptmp(cnv_fill(i).salsec(II), cnv_fill(i).tempsec(II), cnv_fill(i).pressure(II),0);
% recompute conductivity (so that conductivity would give the same S)
    cnv_fill(i).condpri(II)=sw_c3515*sw_cndr(cnv_fill(i).salpri(II),cnv_fill(i).temppri(II),cnv_fill(i).pressure(II))/10;
    cnv_fill(i).condsec(II)=sw_c3515*sw_cndr(cnv_fill(i).salsec(II),cnv_fill(i).tempsec(II),cnv_fill(i).pressure(II))/10;
% recompute sigma
    cnv_fill(i).sigmapri(II)=sw_pden(cnv_fill(i).salpri(II), cnv_fill(i).temppri(II), cnv_fill(i).pressure(II),0)-1000.;
    cnv_fill(i).sigmasec(II)=sw_pden(cnv_fill(i).salsec(II), cnv_fill(i).tempsec(II), cnv_fill(i).pressure(II),0)-1000.;

  
% 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,II),II)) =     2222222;
    cnv_fill(i).quality(II) = cnv_fill(i).quality(II) + 5005500;
    cnv_fill(i).quality(II) = cnv_fill(i).quality(II) +    550000;
    cnv_fill(i).quality(II) = cnv_fill(i).quality(II) +        55;
    cnv_fill(i).quality1oxy(II) = cnv_fill(i).quality1oxy(II) + 4;
    cnv_fill(i).quality2oxy(II) = cnv_fill(i).quality2oxy(II) + 4;
    cnv_fill(i).quality1cond(II) = cnv_fill(i).quality1cond(II) + 4;
    cnv_fill(i).quality2cond(II) = cnv_fill(i).quality2cond(II) + 4;
    cnv_fill(i).quality1cond(II) = cnv_fill(i).quality1tem(II) + 4;
    cnv_fill(i).quality2cond(II) = cnv_fill(i).quality2tem(II) + 4;
% 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))))) ' from pr = ' num2str(max(cnv_fill(i).pressure(II))) '\n']);
    fprintf(nf, ['  Theta    = ' num2str(cnv_fill(i).thetapri(find(cnv_fill(i).pressure ==max(cnv_fill(i).pressure(II))))) ' from pr = ' num2str(max(cnv_fill(i).pressure(II))) '\n']);
    fprintf(nf, ['  Oxygen   = ' num2str(cnv_fill(i).oxmllpri(find(cnv_fill(i).pressure ==max(cnv_fill(i).pressure(II))))) ' from pr = ' num2str(max(cnv_fill(i).pressure(II))) '\n']);

    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))))) ' from pr = ' num2str(max(cnv_fill(i).pressure(II))) '\n']);
    fprintf(nf, ['  Theta    = ' num2str(cnv_fill(i).thetapri(find(cnv_fill(i).pressure ==max(cnv_fill(i).pressure(II))))) ' from pr = ' num2str(max(cnv_fill(i).pressure(II))) '\n']);
    fprintf(nf, ['  Oxygen   = ' num2str(cnv_fill(i).oxmllpri(find(cnv_fill(i).pressure ==max(cnv_fill(i).pressure(II))))) ' from pr = ' num2str(max(cnv_fill(i).pressure(II))) '\n']);

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

      
 
      
      
      
    
     
     
     

  