function [freq_ElNino freq_LaNina] = func_nino( t, sst, i_type, i_plot );
%==========================================================================
%
% Version 1.00
%
%   Compute NINO index from SST time series
%
% This program requires 
%       zerodown.m: computing zero-down crossing wave height from data
%
%  Input
%   t       time series by month
%   sst     sst data
%   i_type  1 - NINO3.4 (+0.4C)
%           2 - NINO3 and 4 (+0.5C)
%   i_plot  1 - Plot for test
%           Othre - no
%
%  Output 
%   freq_ElNino     frequency of El Nino per year
%   freq_LaNina     frequency of La Nina per year
%
%==========================================================================
%
% Example to run this program
%
%    [freq_ElNino freq_LaNina] = func_nino( t, sst, 2, 1 )
%
%==========================================================================
%
% Terms:
%
%       Distributed under the terms of the terms of the BSD License
%
% Copyright:
%
%       Nobuhito Mori
%           Disaster Prevention Research Institue
%           Kyoto University, JAPAN
%           mori@oceanwave.jp
%
%==========================================================================
%
%  Update 2010/09/10 first version
%
%==========================================================================

% % for test
% fname = 'zos_test.mat';
% load(fname)
% t   = 12*time_year;
% sst = SST_nino34;
% %sst = SST_nino3;
% i_type = 1;
% i_plot = 1;

%
% --- preprocess
%

nt = length(t);

% detrend
sst_detrend = detrend(sst,'linear');

% 5 month moving avaerage
a=1; 
n=5;
b = ones(1,n)/n;
nino = filter(b,a, sst_detrend);

%
% --- compute NINO index
%

if i_type == 1
    % NINO3.4 index +0.4
    temp_thresh = 0.4;
else
    % NINO3.4 index exceeds +0.5C
    temp_thresh = 0.5;
end

% number of El Nino by NINO3.4
[h ac at p nw] = zerodown( nino, t );
tc = cumsum(p);
I_elnino = find( ac>=+temp_thresh & p>=6.0 );
I_lanina = find( at<=-temp_thresh & p>=6.0 );

number_ElNino = length(I_elnino);
number_LaNina = length(I_lanina);
freq_ElNino   = length(I_elnino)/(t(nt)/12);
freq_LaNina   = length(I_lanina)/(t(nt)/12);

%
% --- plot
%

if i_plot==1

figure(1);clf
plot(t, nino, 'k:');
hold on
%plot(time_month, SST_nino34_detrend, 'g:');
plot(tc, ac, 'r-');
plot(tc, at, 'b-');
plot(tc(I_elnino), ac(I_elnino), 'ro', 'LineWidth', 2);
plot(tc(I_lanina), at(I_lanina), 'bo', 'LineWidth', 2);
hold off
xlabel('month');
ylabel('NINO Index');

end


