function [submax,submaxchla,submaxdepth]=findsubmax_tp(ELEP_CHLA,ELEP_DEPTH,...
   SURFCHLA,threshold,method,M); %

% EMAPEX Chl-a data shows a lot of spiking and or enhanced Chl-a layers. For many
% profiles it is hard to define one deep Chl-a maxima.

% Doing a running mean of profiles first (before computing the
% maxima) will give a more conservative estimate of number of profiles with
% subsurface maxima.

% Threshold is a Chl-a value such that a submax is consider when 
% ((maxima-smoothEMAPEX_CHLA(1,i))>0.1)

% Modified to assign DCMD to the base of the homogeneous layer if Chl-a
% profile is nearly homogeneous. Criteria?

% Modified to use surface Chla (mean over the first optical depth, instead
% of the closest value to the surface (i.e.(maxima-smooth_CHLA(1,i))>threshold))

% modified to test for sensitivity to the criteria used, i.e. threshold
% criteria vs a relative threshold value
%% Smooth data 
CHLA=ELEP_CHLA; % CHLA_interp; %
DEPTH=ELEP_DEPTH; %DEPTH_interp; %
[nd,np]=size(CHLA);

if strcmp(method,'runmean');
    % M is half the window size 2*M+1
    
    int_CHLA=nan(nd,np);
    smooth_CHLA=nan(nd,np);
    
    for i=1:np
        
        int_CHLA(:,i)=naninterp(CHLA(:,i),[],[]);
        
        indnan=find(isnan(CHLA(:,i)));
        
        if ~isempty(indnan)
            indnan=indnan(1);
            smooth_CHLA(1:indnan-1,i)=runmean(int_CHLA(1:indnan-1,i),M);
        else
            smooth_CHLA(:,i)=runmean(int_CHLA(:,i),M);
        end
    end
else % if 'runmean' is not specify use raw matrices
    smooth_CHLA=CHLA;
end

% Extrapolate values to the top (and bottom)
for i=1:np
    if isnan(CHLA(1,i))
        CHLA(:,i)=naninterp(CHLA(:,i),'nearest',1);
    end
end
SURFCHLA=CHLA(1,:);
% This works best, otherwise if there are NaNs at the top, even SURFCHLA10m could
% be NaN and therefore profile is skipped but may show subsurface maxima

%%    
    % Save Memory
    submax=nan(np,1);
    submaxdepth=nan(np,1);
    submaxchla=nan(np,1);
    
    for i=1:np
        
        % Find the largest maxima in the profile (skipping the surface)
        [maxima, indmax]=nanmax(smooth_CHLA(:,i));
        
        if indmax==1 %  maxima is at the surface || isempty(maxima)==1
            submax(i)=NaN;
%       elseif ((maxima-SURFCHLA(i))>threshold) %% figures saved as _t0.1.eps %((maxima-smooth_CHLA(1,i))>threshold)
       elseif ((maxima-SURFCHLA(i))>threshold.*SURFCHLA(i)) % figures saved as _tp0.1.eps 
            submax(i)=1;
            submaxdepth(i)=DEPTH(indmax(1),i);
            submaxchla(i)=smooth_CHLA(indmax,i);
%         elseif isnan(smooth_CHLA(2,i))~=1 ; %isempty(maxima)==1 %% profile is nearly homogeneous && abs(diff(maximachla(maximadepth<EMAPEX_MLD(i)+20)))<=0.3 
%             
%             
%                 % reassign maxima to base of the mixed layer
%                 % indbml=find(abs(diff(exp(EMAPEX_CHLA(:,i)))));
%                 for j = 1:length(smooth_CHLA(:,i));
%                     % calculate the mixed layer using a threshold value
%                     if abs(smooth_CHLA(2,i)-smooth_CHLA(j,i))>0.1
%                         indbml=j;% Depth(j,i);
%                         break
%                     end
%                 end
%                 submax(i)=0; % profile is homogeneous
%                 submaxchla(i)=CHLA(indbml(1),i);
%                 submaxdepth(i)=DEPTH(indbml(1),i);
% 
% %         else
% %             submax(i)=NaN;
        end
        
        clear maxima indmax
        
%         %% Plot raw profile, smoothed profile and location of subsurface
%         % maxima assigned
%       
%         clf
%         
%         plot(CHLA(:,i),DEPTH(:,i),'--','linewidth',2)
%         hold on        
%         plot(smooth_CHLA(:,i),DEPTH(:,i),'linewidth',2)
%         hold on
%         plot(submaxchla(i),submaxdepth(i),'o','markerfacecolor','r','markersize',10)
%         hold on
%         
% %         % search for MLD
% %         indMLD(i)=searchclosest(EMAPEX_DEPTH(:,i),EMAPEX_MLD(i));
% %         chlaatMLD(i)=smoothEMAPEX_CHLA(indMLD(i),i);
% %         
% %         plot(chlaatMLD(i),EMAPEX_MLD(i),'o','markerfacecolor',[0.2 0.8 0.2],'markersize',10)
%         axis('ij')
%         ylim([0 200]);
%         xlim([0 6]);
%         set(gca,'fontsize',16);
%         ylabel('Depth [m] ','fontsize',16);
%         xlabel('Chl-a [\mug/l] ','fontsize',16);
%         title(['Profile number = ', num2str(i)])
%         % maximize;
%         %print(['ChlaProfile_EMAPEX_', num2str(i),'.eps'],'-depsc2')
%          
%         drawnow
%         pause(0.1)
        
    end