function [trend_out,trend_ste_out,data_nino,time_nino,trend_pvalue_out]=compute_nino_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) & s0==strength_criteria-2);
          elseif strength_criteria==-1;
              inda=find(time_ann>=datenum(perall(1),1,1) & time_ann<=datenum(perall(2),1,1) & s0<=strength_criteria);
          else
              inda=find(time_ann>=datenum(perall(1),1,1) & time_ann<=datenum(perall(2),1,1) & s0>=strength_criteria);
          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));
      % 90% confidence interval
      ci = coefCI(lm,0.32); %0.1 
      trend_ste_out=(ci(2,2)-ci(2,1))/2;
      % 95% confidence interval
      %trend_ste_out=table2array(lm.Coefficients(2,2));
      trend_pvalue_out=table2array(lm.Coefficients(2,4));
   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) & s0==strength_criteria-2);
          elseif strength_criteria==-1;
             inda=find(time_ann>=datenum(perall(1),1,1) & time_ann<=datenum(perall(2),1,1) & s0<=strength_criteria);
          else
             inda=find(time_ann>=datenum(perall(1),1,1) & time_ann<=datenum(perall(2),1,1) & s0>=strength_criteria);
          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_out=trend_md(data_nino,yr,'dim',3,'omitnan'); %unit: trend per year
      trend_ste_out=nan;
      trend_pvalue_out=nan;
   end
end


