function [ZEU]=zeu(Chla,Depth,method);
% Magdalena Carranza (SIO)
% maucarranza@ucsd.edu
[nz, np]=size(Chla);

% Input:
% - Chla and Depth are matrices with profile data laid in columns.
% nz: number of depth levels
% np: number of profiles

% Computation of euphotic depth
ZEU=nan(np,1);
flg_nozeu=nan(np,1);
SUMCHLA_zeu=nan(np,1);


if method==1 % uses method from Morel and Maritorena (2001)
    for i=1:np
        
%             % If using Ctot derived from surface Chl-a (technically Csat or Cpd in
%             % mg/m^3)
%             Csat=nanmean(Chla(1:3,i));
%         
%             if Csat<=0.1
%                 Ctot=38*(Csat)^0.425;
%             elseif Csat>0.1
%                 Ctot=40.2*(Csat)^0.507;
%             end
%         
%             if Ctot>13.65 % mg/m^2, point that corresponds with Zeu=102 m
%                 ZEU(i)=912.5*Ctot^(-0.839);
%             elseif Ctot<=13.65
%                 ZEU(i)=426.3*Ctot^(-0.547);
%             end
%             
%             if ZEU(i)>180
%                 ZEU(i)=182;
%             end
%             % Note it gives similar values of ZEu as in Greg Mitchell's code,
%             % although slithgly deeper (diff of <10m)
        
        % If integrating Chl-a using the iterative process until Zeu is reached
        % (as explained in Morel and Berthon 1989)
        indnan=find(isnan(Depth(:,i)));       
        if nansum(isnan(Depth(:,i)))>0
            indnan=indnan(1);
        else
            indnan=nz+1;
        end
        
%         indnan=find(isnan(Chla(:,i)));
%         if nansum(isnan(Chla(:,i)))>0
%             indnan=indnan(1);
%         else
%             indnan=nz+1;
%         end
        
        for z=2:indnan-1;
            % Ctot (depth-integrated Chl-a should be in units of mg/m^2!!!):
            % IntChla(z)=(Depth(z,i)-Depth(1,i))*((Chla(1,i)+Chla(z,i))./2); % trapezoidal integration
            IntChla(z)=trapz(Depth(1:z,i)',Chla(1:z,i)');
            % nansum(Chla(1:z,i))*Depth(z,i); % Bug in my code? Ctot should be in units of mg/m^2
            
%             if Depth(z,i)<102 % For Zeu<102 m
%                 Ze(z)=912.5*IntChla(z)^(-0.839);
%             elseif Depth(z,i)>102 && Depth(z,i)<180% For Zeu>102 m
%                 Ze(z)=426.3*IntChla(z)^(-0.547);
%             elseif Depth(z,i)>=180
%                 Ze(z)=180;
%                 ZEU(i)=Ze(z);
%             end

            if IntChla(z)>13.65 % For Zeu<102 m
                Ze(z)=912.5*IntChla(z)^(-0.839);
            elseif IntChla(z)<13.65
                Ze(z)=426.3*IntChla(z)^(-0.547);
            end
            
            if Ze(z)<Depth(z,i)
                ZEU(i)=Ze(z-1);
                
                if ZEU(i)>=180
                ZEU(i)=182; % asymptotic value for pure water
                end
                break
            end
        end
    end
    
    
elseif method==2
    for i=1:np
        
%             % If using Ctot derived from surface Chl-a (technically Csat or Cpd in
%             % mg/m^3)
%             Csat=nanmean(Chla(1:3,i));
%         
%             if Csat<=0.1
%                 Ctot=38*(Csat)^0.425;
%             elseif Csat>0.1
%                 Ctot=40.2*(Csat)^0.507;
%             end
%         
%             if Ctot>10 % mg/m^2, point that corresponds with Zeu=102 m
%                 ZEU(i)=568.2*Ctot^(-0.746);
%             elseif Ctot<=10
%                 ZEU(i)=200*Ctot^(-0.293);
%             end
%             % Note it gives similar values of ZEu as in Greg Mitchell's code,
%             % although slithgly shallower (diff of <10m)
        
        
        % If using trapezoidal integration of Chl-a
        indnan=find(isnan(Depth(:,i)));
        if nansum(isnan(Depth(:,i)))>0
            indnan=indnan(1);
        else
            indnan=nz+1;
        end
        
%         indnan=find(isnan(Chla(:,i)));
%         if nansum(isnan(Chla(:,i)))>0
%             indnan=indnan(1);
%         else
%             indnan=nz+1;
%         end
        
        for z=2:indnan-1;
%             % IntChla(z)=nansum(Chla(1:z,i))*Depth(z,i);
%             IntChla(z)=trapz(Depth(1:z,i)',Chla(1:z,i)'); %trapezoidal integration
%             if Depth(z,i)<102 % For Zeu<102 m
%                 Ze(z)=568.2*IntChla(z)^(-0.746); % as in Morel and Berthon 89
%             elseif Depth(z,i)>102 % For Zeu>102 m
%                 Ze(z)=200.0*IntChla(z)^(-0.293); % as in Morel and Berthon 89
%             end
            
            % Using Ctot value as midpoint
            IntChla(z)=trapz(Depth(1:z,i)',Chla(1:z,i)'); %trapezoidal integration
            if IntChla(z)>10 % mg/m^2, point that corresponds with Zeu=102 m
                Ze(z)=568.2*IntChla(z)^(-0.746); % as in Morel and Berthon 89
            elseif IntChla(z)<10 % For Zeu>102 m
                Ze(z)=200.0*IntChla(z)^(-0.293); % as in Morel and Berthon 89
            end
            
            if Ze(z)<Depth(z,i)
                ZEU(i)=Ze(z-1);
                SUMCHLA_zeu(i)=IntChla(z-1);
                if ZEU(i)>182
                    ZEU(i)=182; % for pure water according to MB89
                end
                break
            end
        end
    end
    
    % In some profiles, no ZEU because Ze(z)<Depth(z,i) is not reached
    % extrapolating and recomputing
    % plot(ZEU);hold on; plot(flg_nozeu.*ZEU,'r*')
end



%%  %% ----- Compute euphotic depth ZEU (as in Morel and Maritorena, 2001) ----
%     % see also Morel and Berthon, 1989
%     for i=1:np
%         indnan=find(isnan(Depth(:,i)));
%         if nansum(isnan(Depth(:,i)))>0
%             indnan=indnan(1);
%         else
%             indnan=nz+1;
%         end
%
%         for z=1:indnan-1;
%             IntChla(z)=nansum(Chla(1:z,i)).*Depth(z,i);
%             if Depth(z,i)<102 % For Zeu<102 m
%                 Ze(z)=912.5*IntChla(z)^(-0.839);
%
%             elseif Depth(z,i)>102 % For Zeu>102 m
%                 Ze(z)=426.3*IntChla(z)^(-0.547);
%             end
%
%             if Ze(z)<Depth(z,i)
%                 ZEU(i)=Ze(z-1);
%                 SUMCHLA_zeu(i)=IntChla(z-1);
%                 break
%             end
%         end

%% References:
% - Morel, A., and S. Maritorena (2001), Bio?optical properties of oceanic waters: 
% A reappraisal, Journal of Geophysical Research: Oceans (1978--2012), 106(C4), 7163?7180, 
% doi:10.1029/2000JC000319.
% - Morel, A., and J.-F. Berthon (1989), Surface pigments, algal biomass profiles, 
% and potential production of the euphotic layer: Relationships reinvestigated in 
% view of remote-sensing applications, Limnol. Oceanogr, 34(8), 1545?1562.
