function [annual,seasonal]=GISS_annual(temperatures);



% Calculates the GISS seasonal and annual means

% by Jean S (9/2007)

% Assumes full year input starting in Jan



years=length(temperatures)/12;

temp=reshape(temperatures,12,years)';

seasonal=ones(years,4)*NaN;

annual=ones(years,1)*NaN;



% calculate the reference means and anomalies

ref_mean=nanmean(temp);

anomal=temp-repmat(ref_mean,years,1);



% DJF

seasonal(2:end,1)=mean([temp(1:end-1,12) temp(2:end,1:2)],2);

season=[NaN anomal(1,1:2); anomal(1:end-1,12) anomal(2:end,1:2)];

I=sum(isnan(season),2)==1;

seasonal(I,1)=nanmean(season(I,:)')'+repmat(mean(ref_mean([1:2 12])),sum(I),1);



% MAM

seasonal(:,2)=mean(temp(:,3:5),2);

season=anomal(:,3:5);

I=sum(isnan(season),2)==1;

seasonal(I,2)=nanmean(season(I,:)')'+repmat(mean(ref_mean(3:5)),sum(I),1);



% JJA

seasonal(:,3)=mean(temp(:,6:8),2); 

season=anomal(:,6:8);

I=sum(isnan(season),2)==1;

seasonal(I,3)=nanmean(season(I,:)')'+repmat(mean(ref_mean(6:8)),sum(I),1);



% SON

seasonal(:,4)=mean(temp(:,9:11),2); 

season=anomal(:,9:11);

I=sum(isnan(season),2)==1;

seasonal(I,4)=nanmean(season(I,:)')'+repmat(mean(ref_mean(9:11)),sum(I),1);



% Annual

annual=mean(seasonal,2);

ref_mean_season=nanmean(seasonal);

season=seasonal-repmat(ref_mean_season,years,1);

I=sum(isnan(season),2)==1;

annual(I)=nanmean(season(I,:)')'+repmat(mean(ref_mean_season),sum(I),1);