function [O1,O2,MWmean]=onset_agro_bis(X,lseason,sw,wet,sd,dry,window)

% [O1,O2,MWmean]=onset_agro_bis(X,lseason,sw,wet,sd,dry,window)
%
% This function serves to define the onset of a rainy season from an
% "agronomic" approach, that is the onset is the first day of a sequence of
% 'sw' wet days receving at least 'wet' mm, provided that no dry spells of
% 'sd' days receiving less or equal than 'dry' mm occurs in the following
% 'window' days. 
%
% input
% 'X' : 2D matrix with seasons (usually a window that include a possible
% onset date) beneath each others in row and columns describing stations or
% grid-points. The computations are done independently on each column and
% year.
% 'lseason' : integer giving the length of each season
% 'sw' : integer giving the length of wet sequence for defining onset
% 'wet' : real giving the amount to define the onset. If 'wet' is empty,
% the amount to reach is computed as the mean rainfall for wet sequences of
% 'sw' length during the available period. йIt is computed for each station
% 'sd' : integer giving the length of dry sequence
% 'dry' : real giving the amount to define a dry sequence
% 'window' : integer giving the window after the first onset where the dry
% sequence is searched for
%
% output
% 'O1' : 2D matrix giving for each station and season (in row) the first
% day of a wet sequence receiving 'wet' mm
% 'O2' : 2D matrix giving for each station and season (in rox) the first
% day of the second wet sequence receiving 'wet' mm after a dry sequence
% 'MWmean' : 1D vector giving the threshold for the initial wet sequence
% O1 and O2 are filled with NaN when they are not defined and a sole onset
% could be defined as nanmax(O1,O2) for each season-station. for example,
% onset=reshape(nanmax([O1(:)';O2(:)']),nyear,C);
%
% routines needed
% sequence_overlap.m 
%
% Vincent Moron
% May 2007
%
% modifications
% July 2007 : the onset is set on the first wet day of the wet sequence
% instead of the first day of the wet sequance (that is usually dry !). The
% function is also fasten through the use of the function "while"
% December 2007 : the onset could be defined as the last wet days if the
% criteria of a dry sequence is never filled. These onset should be defined
% as NaN
% January 2008 : The function is faster than before by eliminating few loops. A
% possible bug for dry stations never realizing the second criteria is now
% fixed
% June 2008 : A flexible approach is added for computing the amount of
% rainfall during the wet sequence (it is the mean rainfall for wet
% sequences lasting 'sw' days for each station.

[N,C]=size(X);
nyear=N/lseason;

W=zeros(size(X));
W(find(X > 0))=ones(size(find(X > 0))); % matrix of wet days > 0 mm

if round(nyear) ~= nyear
йййerror('the number of rows should equal length of season x number of season. Please check');
end

if sw > 1;
йййswet=sequence_overlap([1:lseason]',lseason,sw);
йййswet=swet(sw:lseason,:);
йййswet=swet';
йййswet=(swet(:)*ones(1,C))+(ones((lseason-(sw-1))*sw,1)*[0:lseason:lseason*(C-1)]);
йййswet=reshape(swet,sw,C*(lseason-(sw-1)));
end

if sd > 1;
йййsdry=sequence_overlap([1:lseason]',lseason,sd);
йййsdry=sdry(sd:lseason,:);
йййsdry=sdry';
йййsdry=(sdry(:)*ones(1,C))+(ones((lseason-(sd-1))*sd,1)*[0:lseason:lseason*(C-1)]);
йййsdry=reshape(sdry,sd,C*(lseason-(sd-1)));
end

O1=NaN*ones(nyear,C);
O2=NaN*ones(nyear,C);

S=window-(sd-1);
S2=sequence_overlap([1:lseason]',lseason,S);
S2=S2(S:lseason,:)';

% computation of a station-dependent threshold based on the mean rainfall
% in sw days for wet sequence > 1 mm

Lw=lseason-(sw-1);
for i=1:nyear
йййsample=X(((i-1)*lseason)+1:i*lseason,:);
йййif sw > 1
ййййййSWmean(((i-1)*Lw)+1:Lw*i,:)=reshape(sum(sample(swet)),lseason-(sw-1),C); % cumulative rainfall over nday
йййelse
ййййййSWmean(((i-1)*Lw)+1:Lw*i,:)=sample;
йййend
end

for i=1:C;
йййMWmean(i)=mean(SWmean(find(SWmean(:,i) > 1),i));
end

% MWmean is the mean rainfall amount for wet sequence of 'sw' days > 1 mm
if isempty(wet);
йййwet=MWmean;
else
йййwet=wet*ones(1,C);
end


for i=1:nyear
йййclear sample
йййsample=X(((i-1)*lseason)+1:i*lseason,:);
йййwsample=W(((i-1)*lseason)+1:i*lseason,:);
йййif sw > 1
йййййййSW=reshape(sum(sample(swet)),lseason-(sw-1),C); % cumulative rainfall over nday
йййelse
йййййййSW=sample;
йййend
йййif sd > 1
йййййййSD=reshape(sum(sample(sdry)),lseason-(sd-1),C); % cumulative rainfall over nday_dry
йййelse
йййййййSD=sample;
йййend
ййй[nrw,ncw]=size(SW);
ййй[nrd,ncd]=size(SD);
йййfor j=1:C;
йййййййclear tab o1
йййййййtab=[sample(:,j) wsample(:,j) [SW(:,j);ones(sw-1,1)*SW(lseason-sw+1,j)] [SD(:,j);zeros(sd-1,1)]];
ййййййй[nrtab,nctab]=size(tab);
йййййййo1=find(tab(:,3) >= wet(j) & tab(:,2) == 1);
йййййййO2(i,j)=NaN;
йййййййO1(i,j)=NaN;
йййййййD=tab(:,4);
йййййййD=D(S2);
йййййййD=D';
йййййййD=[D;zeros(window-sd,S)]; ййййййй
йййййййif ~isempty(o1);
йййййййййййO1(i,j)=o1(1);
йййййййййййtab2=D(o1,:);
йййййййййййo2=o1(find(min(tab2') > dry));
йййййййййййif ~isempty(o2);
йййййййййййййййO2(i,j)=o2(1);
йййййййййййend
йййййййend
йййend
end


