function [nino_strength]=compute_strength(sst_nino,rain_nino);

   [mm,lm]=size(sst_nino);
   if lm<=10;
       disp(['observational data nino strength']);
       nino_strength=zeros(mm,lm);
       yr_mod=[1986 1987 1991 1994 2002 2009 2015];yr_ext=[1997];
       [yr,mo,da]=datevec(rain_nino);
       for i=1:length(yr_mod);
           ind=find(yr==yr_mod(i));
           nino_strength(:,ind)=1;
       end
       ind=find(yr==yr_ext);
       nino_strength(:,ind)=2;
   else
       for m=1:mm
           nino_strength0=zeros(1,lm);
           sst_nino0=sst_nino(m,:);
           rain_nino0=rain_nino(m,:);
   
           %moderate El Nino is defined in historical period, without considering extreme nino
           sst_nino1=sst_nino0;
           %sst_nino1(rain_nino0>=5)=nan;
           if lm==37 | lm==36;  
              sst_nino_std=nanstd(sst_nino1(1:lm));
           elseif lm==200;
              %sst_nino_std=nanstd(sst_nino1(82:119));
              sst_nino_std=nanstd(sst_nino1(1:100));
           elseif lm==251;
              sst_nino_std=nanstd(sst_nino1(51:150));
           else
              error('time step of sst is wrong');
           end
           %if lm>=200; sst_nino_std=nanstd(sst_nino1(1:lm/2)); end
           coef_std=1;
           nino_strength0(sst_nino0>=coef_std*sst_nino_std)=1;
           nino_strength0(sst_nino0<=-1*coef_std*sst_nino_std)=-1;
           nino_strength0(sst_nino0<=-1.75*coef_std*sst_nino_std)=-2;
           nino_strength0(rain_nino0>=5)=2;
           %nino_strength0(rain_nino0>=5 & sst_nino0>=coef_std*sst_nino_std)=2;
           nino_strength(m,:)=nino_strength0;
       end
   end
end


