function [ZEU]=zeu(Chla,Depth,method);

[nz, np]=size(Chla);

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


if method==1
    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
        
        % In MM2001 equations are good for Zeu<180m
        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)
        
    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
        
        if ZEU(i)>182
            ZEU(i)=182;
        end
        % Note it gives similar values of ZEu as in Greg Mitchell's code,
        % although slithgly shallower (diff of <10m)
    end
elseif method==3
    
    for i=1:np
        
        % If using Ctot derived from surface Chl-a (technically Csat or Cpd in
        % mg/m^3)
        Csat(i)=nanmean(Chla(1:3,i));

        if Csat(i)<=0.3
            ZEU(i)=63.442*(Csat(i))^(-0.176);
        elseif Csat(i)>0.3
            ZEU(i)=48.717*(Csat(i))^(-0.417);
        end
    end
end


% Notes:
