function [anom,clim]=anomaly_ts(Z)
%substract climatology from an array
% expects a field Z(t) with t in months 
% By default, 1 is the first January month.
% MAKE SURE IT IS THE CASE in your data
% inputs :
%            Z (t) 
% outputs : anom(t), clim(12),
%  
%  
%  Julien Emile-Geay, 07/12/2006
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nt=length(Z);
clim=zeros(12,1);
anom=zeros(nt,1);
step=[0:12:nt-1]; % jump by 12 months as a time
for month=1:12
    if ( month <= mod(nt,12) ) 
        add_1=12;  
    else add_1=0;
    end
    T_ind=month+step+add_1; %add the last month for the incomplete year
    clim(month)= nmean(Z(T_ind)); 
    anom(T_ind)=Z(T_ind) - clim(month);
end