function [anom,clim]=anomaly(array)
%substract climatology from an array
% expects a field F(T,Y,X) with T in months 
% By default, 1 is the first January month.
% MAKE SURE IT IS THE CASE in your data
% inputs :
%            array (T,Y,X) 
% outputs : anom, clim
%
%  Hepta Technologies, Feb 2007
%    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 

%define grid vectors
nx=size(array,3);
ny=size(array,2);
nt=size(array,1); % in months since 1960
clim=zeros(12,ny,nx);
anom=zeros(size(array));
step_full=12*[0:1:ceil(nt/12)-1]; % jump by 12 months as a time
step_partial=step_full(1:end-1);  % will fail if you try this on a two-year record, but who
%would be cracked enough to do that anyway ? I know a good psychiatrist if you need one. 
%
if (mod(nt,12) ~=0)  % if there are incomplete years
   % treat full years
   step= step_full;
   for month=1:mod(nt,12)
      T_ind=month+step; 
      tmp=squeeze(array(T_ind,:,:)); %select only relevant months
      clim(month,:,:)=nmean(tmp,1); % arithmetic average along first dim (time)
      anom(T_ind,:,:)=array(T_ind,:,:)-repmat(clim(month,:,:),[length(T_ind),1,1]);
   end
   % treat partial years
   step= step_partial;
   for month=mod(nt,12)+1:12
      T_ind=month+step; 
      tmp=squeeze(array(T_ind,:,:)); %select only relevant months
      clim(month,:,:)=nmean(tmp,1); % arithmetic average along first dim (time)
      anom(T_ind,:,:)=array(T_ind,:,:)-repmat(clim(month,:,:),[length(T_ind),1,1]) ;
   end
else  % else if everything is good, sit back and enjoy
   for month=1:12
      T_ind=month+step_full;
      tmp=squeeze(array(T_ind,:,:)); %select only relevant months
      clim(month,:,:)=nmean(tmp,1); % arithmetic average along first dim (time)
      anom(T_ind,:,:)=array(T_ind,:,:)-repmat(clim(month,:,:),[length(T_ind),1,1]) ;
   end   
end