function [trend_out,trend_ste_out,data_nino,time_nino]=compute_all_trend(data,nino_strength,time_ann,indm,perall,strength_criteria);

   jm=length(indm);
   data_nino=[]; time_nino=[];
   sz=length(size(data));
   if sz==2;
      for j=1:jm
          s0=nino_strength(indm(j),:);
          if strength_criteria==3;
              inda=find(time_ann>=datenum(perall(1),1,1) & time_ann<=datenum(perall(2),1,1));
          else
              inda=find(time_ann>=datenum(perall(1),1,1) & time_ann<=datenum(perall(2),1,1));
          end
          f1=squeeze(data(indm(j),inda));
          data_nino=cat(2,data_nino,f1);
          t1=time_ann(inda);
          time_nino=cat(2,time_nino,t1);
      end
      [yr,mo,da]=datevec(time_nino);
      lm = fitlm(yr,data_nino);
      trend_out=table2array(lm.Coefficients(2,1));
      trend_ste_out=table2array(lm.Coefficients(2,2));
   elseif sz==4;
      for j=1:jm
          s0=nino_strength(indm(j),:);
          if strength_criteria==3;
             inda=find(time_ann>=datenum(perall(1),1,1) & time_ann<=datenum(perall(2),1,1));
          else
             inda=find(time_ann>=datenum(perall(1),1,1) & time_ann<=datenum(perall(2),1,1));
          end
          f1=squeeze(data(indm(j),:,:,inda));
          data_nino=cat(3,data_nino,f1);
          t1=time_ann(inda);
          time_nino=cat(2,time_nino,t1);
      end
      [yr,mo,da]=datevec(time_nino);
      [am,bm,cm]=size(data_nino);
      trend_out=trend_md(data_nino,yr,3); %unit: trend per year
      trend_ste_out=nan;
   end
end


