function [ds]=mdrought(di,mthd,mn)
% MDROUGHT Finds the megadroughts in a given time series, 'di'.
%    ds = mdrought(di,mthd,mn) where 'di' is the drought index,'mthd'
%    indicates the method for computing the droughts  and 'mn' is
%    an option for removing the time mean.
%
%    This function outputs a matrix with three columns. The first two
%    columns are the start and end years of the drought and the last 
%    column is the strength of the drought, computed as the sum of the 
%    drought index over the period of the drought. Droughts are sorted 
%    by the strength of the drought, beginning with the most severe.
%    Droughts that aren't at least 10 years long are thrown out.
%
%    The first method assumes that drought = negative anomalies around 
%    11 year running mean. The second method is the 2 start 2 end 
%    criteria, where droughts start with two consecutive years of dry
%    anomalies and are ended with two consecutive wet anomalies.


drghts=struct([]); % initialize structure

if mn=='y'
   dinm=di(:)-mean(di(:)); % Remove the time mean
else
   dinm=di(:);
end


if mthd==1

   dr_v=smooth(dinm,11); % smooth it at 11 years running mean
   neg_ids=find(dr_v<0); % Get the negative indices
   negpos=ones(length(di),1); % positive years will be ones
   negpos(neg_ids)=0;% set negative years as zeros

   i=1;j=1;flag=0;
   while i+1 <= length(di)
     % Find the start of a drought
     if negpos(i)==0
	 d_start=i;
	 i=i+1;
	 % Keep going to the end of a drought
	 while negpos(i)~=1
	     i=i+1;
	     if i+1 > length(di)
		 d_end=i-1;
		 drghts(j).strtend=[d_start d_end];
		 drghts(j).dsum=sum(dr_v(d_start:d_end));
		 flag=1;
		 break
	     end
	 end
	 if flag==1;break;end % get out of nested loops
	 
	 d_end=i-1;
	 drghts(j).strtend=[d_start d_end];
	 drghts(j).dsum=sum(dr_v(d_start:d_end));
	 
	 j=j+1;
     end
     
     i=i+1;
   end

   % Rank the droughts, find the most severe
   d=zeros(length(drghts),3);
   for kk=1:length(drghts)
     d(kk,1:2)=drghts(kk).strtend;
     d(kk,3)=drghts(kk).dsum;
   end

   % Throw out droughts that aren't at least 10 years long
   rmv=[];
   ii=1;
   for jj=1:length(drghts)
      if (d(jj,2)-d(jj,1)) < 10
	 rmv(ii)=jj;
	 ii=ii+1;
      end
   end
   d(rmv,:)=[];  


   % Sort the droughts according to the most severe drought density (sum of
   % length and intensity)
   ds=sortrows(d,3);



elseif mthd==2

   % Find the droughts
       
    dr_v=dinm; 
    neg_ids=find(dr_v<0);
    negpos=ones(length(di),1); % positive years will be ones
    negpos(neg_ids)=0;% set negative years as zeros
    
    i=1;j=1;flag=0;
    while i+1 <= length(di)
	% Find the start of a drought
	if negpos(i)==0 && negpos(i+1)==0
	    d_start=i;
	    i=i+2;
            % Stop if it's the end of the time series
            if i>=length(negpos);break;end
	    % Keep going to the end of a drought
	    while negpos(i)~=1 || negpos(i+1)~=1
		i=i+1;
		if i+1 > length(di)
		    d_end=i;
		    drghts(j).strtend=[d_start d_end];
		    drghts(j).dsum=sum(dr_v(d_start:d_end));
		    flag=1;
		    break
		end
	    end
	    if flag==1;break;end % get out of nested loops
	    
	    d_end=i+1;
	    drghts(j).strtend=[d_start d_end];
	    drghts(j).dsum=sum(dr_v(d_start:d_end));
	    
	    j=j+1;
	end
	
	i=i+1;
    end
    
   % Rank the droughts, find the most severe
   d=zeros(length(drghts),3);
   for kk=1:length(drghts)
     d(kk,1:2)=drghts(kk).strtend;
     d(kk,3)=drghts(kk).dsum;
   end

   % Throw out droughts that aren't at least 10 years long
   rmv=[];
   ii=1;
   for jj=1:length(drghts)
      if (d(jj,2)-d(jj,1)) < 10
	 rmv(ii)=jj;
	 ii=ii+1;
      end
   end
   d(rmv,:)=[];  


   % Sort the droughts according to the most severe drought density (sum of
   % length and intensity)
   ds=sortrows(d,3);

end
