function [mVAR,bin,error,N]=meanprofile_ne(VAR,TEMP,DEPTH,refDepth,plt,clr);

% 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 Chl-a profiles in each column, size [nz,np]
% TEMP: matrix with Temperature profiles in each column, size [nz,np],
%       used to compute the error in the averaging method
% DEPTH: matrix with DEPTH corresponding to each profile, size [nz,np],
% refDepth: is a vector of size [np,1] with the reference depth (e.g MLD,
% DCMD depth, etc)
% plot=1 to plot error with dashed lines
% plot=2 to plot error with shaded colors
% nz is the number of depth levels (lines)
% np, is teh number of profiles (columns)

% OUTPUT variables:
% mVAR = mean profile, [nz,1]
% bin =new adjusted depth axis, [nz,1]
% error = error, [nz,1]
% N = number of data points averaged in each bin (N>30, otherwise not computing a mea)
%
% Non-native MATLAB functions required:
% searchclosest.m
% boundedline.m
%
% %% To Test the function (eg with interpolated data to a common grid):
% VAR=TEMP_interp; %CHLA_interp;
% TEMP=TEMP_interp;
% DEPTH=DEPTH_interp;
% refDepth=EMAPEX_MLDtt(:,4);
% plt=2;
% clr='k';

% %% To Test the function (eg with data that needs to be interpolated):
% VAR=EMAPEX_TEMP; %CHLA_interp;
% TEMP=EMAPEX_TEMP;
% DEPTH=EMAPEX_DEPTH;
% refDepth=EMAPEX_MLDtt;
% plt=2;
% clr='k';

% Delete profiles where there is no MLD estimate
nomld=isnan(refDepth); %
VAR=VAR(:,~nomld);
TEMP=TEMP(:,~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);
    dz=2.5; % ~aprox vertical resolution for EM-APEX
    DEPTH_interp=[floor(RD(1)):dz:ceil(RD(2))]';
    nz=length(DEPTH_interp(:,1));
    VAR_interp=nan(nz,np);
    TEMP_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');
            TEMP_interp(:,i)=interp1(DEPTH(:,i),TEMP(:,i),...
                DEPTH_interp(:,i),'linear');
        elseif isempty(indnan)~=1
            indnan=indnan(1)-1;
            
            VAR_interp(:,i)=interp1(DEPTH(1:indnan,i),VAR(1:indnan,i),...
                DEPTH_interp(:,i),'linear');
            TEMP_interp(:,i)=interp1(DEPTH(1:indnan,i),TEMP(1:indnan,i),...
                DEPTH_interp(:,i),'linear');
        end
        
    end
    
    DEPTH=DEPTH_interp;
    TEMP=TEMP_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);
% mldTEMP=nan(nz,np);

% Construct adjusted depth matrix
for p=1:np;
    
    [mx,indmax]=nanmax(DEPTH(:,p));
    
    if refDepth(p)<DEPTH(indmax,p)
        refind(p)=searchclosest(DEPTH(:,p),refDepth(p));
        adjDEPTH(:,p)=DEPTH(:,p)-DEPTH(refind(p),p);
        
        %Compute mean MLD temperature for each profile
        mldTEMP(p,1)=nanmean(TEMP(1:refind(p)));
    else
        continue
    end
end

% New adjusted depth vector (common Depth axis relative to the MLD)
rng=range(adjDEPTH);
bin=[rng(1):dz:rng(2)]';
nbin=length(bin); % number of depth bins

% Compute mean (or median) profiles
mVAR=nan(nbin,1);
mTEMP=nan(nbin,1);
stdVAR=nan(nbin,1);
N=nan(nbin,1);

% For each depth bin
for z=1:nbin;
    % Find the indeces of all values that need to be averaged
    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)=mean(VAR(ind)); % mean VAR profile
        stdVAR(z,1)=std(VAR(ind)); % standard deviation of VAR for each bin
        
        mTEMP(z,1)=mean(TEMP(ind)); % mean TEMP profile to evaluate error
    end
    
end

% Standard error
error_std=stdVAR./sqrt(N); % in units of VAR

errormethod=error_std.*0;
error=error_std;

%%  Compute error in avg method based on mean of temperature profiles

maxTEMP=nanmax(mTEMP);
mldTEMP=nanmean(mldTEMP); % mean MLD temperature among all profiles in the matrix
mmTEMP=nanmean(mTEMP);

% % Option 1)
% % multiplicative error, by the inverse of the % of change between Temp at
% % one level and maximum temp
% errormethod=1./(100*(maxTEMP-mTEMP)./mTEMP); % more conservative
% errormethod(errormethod==Inf)=NaN;
% error=error_std.*(1+errormethod); % %% Need to multiply the error scale factor by the
%                                    % mean profile, otherwise the error units make
%                                    % no sense
%
% % cd /Users/Mau/Documents/ChlaProfiles/error
% % figure(1)
% % plot(errormethod,bin); % hold on; plot(errormethod2,'r')
% % xlabel('error_m','fontsize',20)
% % axis ij
% % print('Error_Opt1_errorm_maxTrel2mT.eps','-depsc2');
% % close all
%
% error(error==Inf)=NaN;

% %%  Option 2)
% % Additive error
% % errormethod=((maxTEMP-mTEMP)./mTEMP); % doesn't work, error is minimal at the MLD
% % errormethod=1./(abs(maxTEMP-mTEMP)./maxTEMP); % doesn't work
% % errormethod=abs(mTEMP-mmTEMP)./mmTEMP; % less conservative
% errormethod=abs(mTEMP-mldTEMP)./mldTEMP; % more conservative
%
% cd /Users/Mau/Documents/ChlaProfiles/error
% figure(1)
% plot(errormethod,bin); % hold on; plot(errormethod2,'r')
% xlabel('\epsilon','fontsize',20)
% axis ij
% print('Error_Opt2_epsilon_rel2mldT.eps','-depsc2');
% close all
%
% error=error_std+errormethod.*mVAR; % %% Need to multiply the error scale factor by the
%                                    % mean profile, otherwise the error units make
%                                    % no sense
% % error=error_std+(1+errormethod).*error_std; % %% Need to multiply the error scale factor by the
% %                                    % mean profile, otherwise the error units make
% %                                    % no sense
%

% % Option 3)
% errormethod=(maxTEMP-mTEMP).^-1;
% error=sqrt(error_std.^2+errormethod.^2);
%
% cd /Users/Mau/Documents/ChlaProfiles/error
% figure(1)
% plot(errormethod,bin); % hold on; plot(errormethod2,'r')
% xlabel('error_m','fontsize',20)
% axis ij
% print('Errormethod_SarahsOpt1_inverse.eps','-depsc2');
% close all
%
% error(error==Inf)=NaN;

% % Option 4)
% % Additive error
% % errormethod=((maxTEMP-mTEMP)./mTEMP); % doesn't work, error is minimal at the MLD
% % errormethod=1./(abs(maxTEMP-mTEMP)./maxTEMP); % doesn't work
% % errormethod=abs(mTEMP-mmTEMP)./mmTEMP; % less conservative
% errormethod=abs(mTEMP-mldTEMP)./mldTEMP; % more conservative
%
% % cd /Users/Mau/Documents/ChlaProfiles/error
% % figure(1)
% % plot(errormethod,bin); % hold on; plot(errormethod2,'r')
% % xlabel('\epsilon','fontsize',20)
% % axis ij
% % print('Error_SarahsOpt2_epsilon_rel2mldT.eps','-depsc2');
% % close all
%
% error=sqrt(error_std.^2+(errormethod.*mVAR).^2); % %% Need to multiply the error scale factor by the
%                                    % mean profile, otherwise the error units make
%                                    % no sense
% % error=error_std+(1+errormethod).*error_std; % %% Need to multiply the error scale factor by the
% %                                    % mean profile, otherwise the error units make
% %                                    % no sense

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

errormethod=errormethod(~isnan(error));
error_std=error_std(~isnan(error));
error=error(~isnan(error));

%% Figure

if plt==1;
    % Plot mean profile with dotted error lines
    plot(mVAR,bin,'r');
    hold on
    plot(mVAR-error,bin,'r--')
    hold on
    plot(mVAR+error,bin,'r--')
    axis ij
elseif plt==2;
    % Plot mean profiles with shaded error
    boundedline(mVAR,bin,error,clr,'orientation','horiz')
    axis ij
end

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

maximize;

% xlabel('Chla','fontsize',20)
% print('MeanChla_Error_Opt1_errorm__maxTrel2mT.eps','-depsc2');

% xlabel('Temperature','fontsize',20)
% print('MeanTemp_Error_Opt1_errorm__maxTrel2mT.eps','-depsc2');

% xlabel('Chl-a','fontsize',20)
% print('MeanChla_Error_Opt2_epsilon_rel2mldT.eps','-depsc2');

% xlabel('Temperature','fontsize',20)
% print('MeanTemp_Error_Opt2_epsilon_rel2mldT.eps','-depsc2');