function [mVAR,bin,error,N]=meanprofile(VAR,DEPTH,refDepth,plot);

% 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

% 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=TEMP_interp;
% DEPTH=DEPTH_interp;
% refDepth=EMAPEX_MLDtt;

% Delete profles where there is no mld
nomld=isnan(refDepth); % 8 profiles
VAR=VAR(:,~nomld);
DEPTH=DEPTH(:,~nomld);
refDepth=refDepth(~nomld);

% 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(DEPTH(:,p),refDepth(p));
    adjDEPTH(:,p)=DEPTH(:,p)-DEPTH(refind(p),p);
    
    mldVAR(p,1)=nanmean(VAR(1:refind(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);

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)>10 % only compute a mean value if there is more than 30 data points
        mVAR(z,1)=nanmedian(VAR(ind)); % mean profile
        stdVAR(z,1)=nanstd(VAR(ind)); % standar deviation
    end
    
end

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

% % load errormethod based on temperature
% cd /Users/Mau/Documents/ChlaProfiles/emapex
% er=load('AVGmethod_error.mat');
% % Interpolate error method to new depht bins (in case they're different)
% errormethod=interp1(er.bin,er.errormethod,bin,'linear',NaN');
% 
% % error normalized by % of change between Temp at one level and maximum
% % temp
% error=error./errormethod;
% error(error==Inf)=NaN;

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

%%

if plot==1;

% plot(mTEMP,bins,'r');
% hold on
% plot(mTEMP-error,bins,'r--')
% hold on
% plot(mTEMP+error,bins,'r--')
boundedline(mVAR,bin,error,'r','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



