function [dataa,datam]=compute_anomaly(data,time,tc2,tc1);

% process monthly data
% tc2 is large and tc1 is small, e.g., cut off between 13 month and 15 years;
% tc2=12*15 and tc1=13 
    %create bandpass digit filter
    fname_bp='bp.mat';create_bp=false;
    fstat=dir(fname_bp);
    if length(fstat)==0;
       create_bp=true;
    else
       fage=(datenum(clock)-datenum(fstat.date))*24*3600;
       if fage>1800;
          create_bp=true;
       end
    end
    if ~create_bp; load(fname_bp); end

    num=0;
    if length(size(data))==3;
       [jm,im,lm]=size(data);
       dataa=nan(jm,im,lm);datam=nan(jm,im,lm);
        for j=1:jm
            if mod(j,100)==0; disp(['       bandpass filter ',num2str(j),'/',num2str(jm)]); end
            for i=1:im
                d0=squeeze(data(j,i,:));
                if length(d0(~isnan(d0)))==lm;
                   %d0_end_points=mean(d0(end-119:end))*rand(120,1); d0=cat(1,d0,d0_end_points);
                   if num==0 & create_bp
                      num=num+1; 
                      [d1,bp]=bandpass(d0,[1/tc2 1/tc1],1,'Steepness',0.95); %compute bandpass object bp, only run one time
                      [d2,hp]=highpass(d0,1/13,1,'Steepness',0.98);
                      save(fname_bp,'bp','hp') 
                   end
                   d_ia=filtfilt(bp,d0);
                   d_season=filtfilt(hp,d0);
                   %d0_m=smooth(d0-d_ia-d_season,13); dataa(j,i,1:lm)=d_ia(1:end-120); datam(j,i,1:lm)=d0_m(1:end-120);
                   dataa(j,i,1:lm)=d_ia; datam(j,i,1:lm)=smooth(d0-d_ia-d_season,13);
                end
            end
        end
    else
       error('data size is wrong')
    end
end









