function [mVAR,bin,error,N,newVAR1,newDEPTH]=meanprofile_v2(VAR,DEPTH,refDepth,pplot,clr);

% normalizes the profiles between 0-1, and 
% averages using nanmean with no extrapolation.

% Same as meanprofile2.m but allows to input variables that have not been
% previously interpolated or that have been interpolated to different
% grid bins. For example Argo dz=5, and EM-APEX dz=2.5. 
% This code will re-interpolate using the highest resolution.

% Computes mean profiles averaging over an adjusted depth axis relative to a
% reference depth (see e.g. Hense and Beckman (2008):
% Hense, I., and A. Beckmann (2008), Revisiting subsurface chlorophyll and phytoplankton distributions, Deep Sea Research Part I: Oceanographic Research Papers, 55(9), 1193?1199, doi:10.1016/j.dsr.2008.04.009.

% Input variables:
% VAR: matrix with profiles in each column, size [nz,np]
% DEPTH: matrix with DEPTH corresponding to each profile
% refDepth: is a vector of size [np,1] with the reference depth (e.g MLD,
% DCM depth, etc)
% plot=1 to plot, no plot otherwise

% NOTE: DEPTH should have equally spaced bins, and VAR should have
% interpolated values to a common grid

% Feb20 2016: cleaned up + skipping profiles where depth is nor monotonically increasing
% + allowing to interpolate when profiles have repeated values (if depth is not strictly monotonically 
% increasing interp1 doesn't work)

% Output variables:
% mVAR = mean profile
% bin =new adjusted depth axis
% error = standard error
% N= number of data points averaged in each bin (N>30)

% % Test:
% VAR=VAR(:,indm);%TEMP_interp(:,fll);%(:,spring); % ALL_TEMP; %
% DEPTH=DEPTH(:,indm); %DEPTH_interp(:,fll);%(:,spring); % ALL_DEPTH; %
% refDepth=MLD(indm); %MLDtt(fll);%(spring);%(:,4);%(summer,4);ALL_MLD; %EMAPEX_
% pplot=1;
% clr='k';

% Delete profiles where there is no MLD estimate
nomld=isnan(refDepth); %
VAR=VAR(:,~nomld);
DEPTH=DEPTH(:,~nomld);
refDepth=refDepth(~nomld);

% Interpolate to a common depth grid

% if length(unique(diff(DEPTH(:,1))))~=1
    [nz,np]=size(VAR);
    
    RD=range(DEPTH);
    diffDEPTH=diff(DEPTH); diffDEPTH(diffDEPTH==0)=NaN; % added to avoid dz=0 n line below
    dz=abs(nanmin(nanmean(diffDEPTH))); % 2.5 ~aprox vertical resolution for EM-APEX (minimum dz)
    DEPTH_interp=[floor(RD(1)):dz:ceil(RD(2))]';
    nz=length(DEPTH_interp(:,1));
    VAR_interp=nan(nz,np);
    
    for i=1:np
        % Interpolate data to regular depth bins
        DEPTH_interp(:,i)=[floor(RD(1)):dz:ceil(RD(2))]';
        indnan=find(isnan(DEPTH(:,i)));
        
        if isempty(indnan)==1
            VAR_interp(:,i)=interp1(DEPTH(:,i),VAR(:,i),...
                DEPTH_interp(:,i),'linear');% 'nearest','extrap'); % 
        elseif isempty(indnan)~=1
            indnan=indnan(1)-1;
            % Check if Dept is monotonically increasing
            if issorted(DEPTH(1:indnan,i))==1
                % if there are repeated depth values, then this code keeps
                % only one, otherwise interp1 doesn't work
                [D, inddepth] = unique(DEPTH(1:indnan,i),'stable');
                V=VAR(inddepth,i);
                
               % interpolate:
                VAR_interp(:,i)=interp1(D,V,... %DEPTH(1:indnan,i),VAR(1:indnan,i),...
                DEPTH_interp(:,i),'linear');
            elseif issorted(DEPTH(1:indnan,i))==0
                continue
            end
                
        end
        
    end
    
    DEPTH=DEPTH_interp;
    VAR=VAR_interp;
% end

% Matrix dimensions
[nz,np]=size(VAR);

% vertical resolution
dz=nanmean(diff(DEPTH(:,1)));

% save memory
refind=nan(np,1);
adjDEPTH=nan(nz,np);

% Construct adjusted depth matrix
for p=1:np;
    refind(p)=searchclosest_m(DEPTH(:,p),refDepth(p));
    adjDEPTH(:,p)=DEPTH(:,p)-DEPTH(refind(p),p);
end

% New adjusted depth vector
rng=range(adjDEPTH);
bin=[rng(1):dz:rng(2)]'; % sort(reshape(adjDEPTH,nz*np,1));%

nbin=length(bin);

mVAR=nan(nbin,1);
stdVAR=nan(nbin,1);
N=nan(nbin,1);

%%  Construct VAR matrix such that all columns have the same depth level 
% and then I can average over lines

newVAR=nan(nbin,np);
newDEPTH=nan(nbin,np);

for p=1:np;
    indi(p)=searchclosest_m(bin,adjDEPTH(1,p)); %find(bin==adjDEPTH(1,p));
    indf(p)=searchclosest_m(bin,adjDEPTH(end,p));
    newVAR(indi(p):indf(p),p)=VAR(:,p);
    newDEPTH(indi(p):indf(p),p)=adjDEPTH(:,p);
end

% Compute anomalies before averaging? or normalize the variables?
% such that all profiles are comparable?
newVAR1=nan(nbin,np);
 for p=1:np;
     newVAR1(:,p)=(newVAR(:,p)-nanmin(newVAR(:,p)))./(nanmax(newVAR(:,p))-nanmin(newVAR(:,p))); %scaling all profiles between 0,1
%     newVAR(:,p)=(newVAR(:,p)-nanmean(newVAR(:,p)))./nanstd(newVAR(:,p)); % standard normalization
%    newVAR(:,p)=(newVAR(:,p)-nanmean(newVAR,2))./(nanstd(nanmean(newVAR)));
%     -nanmean(newVAR(:,p)))./nanstd(newVAR(:,p),01);
%     nanmean(newVAR,2));%./nanstd(newVAR,0,2);%nanmean(newVAR(:,p));% nanmean(newVAR,2);% ./nanstd(newVAR,0,2); % nanmean(newVAR(:,p));%
 end

% % With Extrapolation above the surface and below the profile's max depth
% % to fill matrix?
% for p=1:np;
%      if nansum(~isnan(newVAR1(:,p)))~=0
%       newVAR1(:,p)=naninterp(newVAR1(:,p),'nearest'); % extrapolates
%       newDEPTH1(:,p)=naninterp(newDEPTH(:,p),'linear');
%      end
% end

% Test mean average to see if it is the same as the method used before
mVAR=nanmean(newVAR1,2); %!!! if using nanmean(newVAR,2); then averaging over lines with NaNs and it gives a subsurface max
stdVAR=nanstd(newVAR1,0,2); % nanstd(newVAR,0,2); that is an artifact of the avg method.. 
N=nansum(double(~isnan(newVAR1)),2);


% Note, it does give the same result, so the code below is OK

% %% Will extrapolate end values up to the top/end of matrix 
% % (as suggested by Sarah) but this assumption only works for temperature
% 
% new_VAR=nan(nbin,np);
% for p=1:np
% %     new_VAR(:,p)=interp1(bin,newVAR(:,p),bin,'nearest','extrap');
%     indi=~isnan(newVAR(:,p));
%     indi=find(indi==1);
%     
%     new_VAR(1:indi(1),p)=ones(indi(1),1)*newVAR(indi(1),p);
%     new_VAR(indi(1):indi(end),p)=newVAR(indi(1):indi(end),p);
%     new_VAR(indi(end):end,p)=ones(nbin-indi(end)+1,1)*newVAR(indi(end),p);
% end
% 
% % Conpute mean average
% mVAR=nanmean(new_VAR,2);
% stdVAR=nanstd(new_VAR,0,2);
% N=nansum(double(~isnan(new_VAR)),2);

% %% First attempt
% 
% for z=1:nbin;
%     ind=find(adjDEPTH==bin(z));
%     N(z,1)=nansum(~isnan(VAR(ind))); % number of data points averaged in each bin
%     
%     if N(z,1)>30 % only compute a mean value if there is more than 30 data points
%         mVAR(z,1)=nanmean(VAR(ind)); % mean profile
%         stdVAR(z,1)=nanstd(VAR(ind)); % standar deviation
%     end
%     
% end

%%  Standard error
error=stdVAR./sqrt(N);

% Delete NaN in mean profiles
mVAR=mVAR(~isnan(error));
stdVAR=stdVAR(~isnan(error));
N=N(~isnan(error));
bin=bin(~isnan(error));
error=error(~isnan(error));

%%

if pplot==1;

% plot(mTEMP,bins,'r');
% hold on
% plot(mTEMP-error,bins,'r--')
% hold on
% plot(mTEMP+error,bins,'r--')
boundedline(mVAR,bin,error,clr,'orientation','horiz')
axis ij

% Change fontsize
set(gca,'FontSize',20); %,'fontWeight','bold')
set(findall(gcf,'type','text'),'FontSize',20); %,'fontWeight','bold')
set(0,'defaultlinelinewidth',1)
 
maximize;

end



