function combined=GISScombine(temperatures);



% combines stations ala GISS

% by Jean S (9/2007)

% assumes integers; and full years starting in Jan



% Minimum overlap for combining

min_overlap=4;



[years,stations]=size(temperatures);



% 1st station is assumed to be the "best"

fixed_stations(:,1)=temperatures(:,1);

no_fixed=1;

ref_station=temperatures(:,1);



% annuals for all stations

for k=1:stations;

    annual(:,k)=GISSannual(temperatures(:,k));

end



available_stations=2:stations;

% repeat until no new records to add

while ~isempty(available_stations),

  [ref_annual,ref_seasonal]=GISSannual(ref_station);

  % Find the annual mean matches

  overlaps=repmat(isfinite(ref_annual),1,length(available_stations)) & ... 

                              isfinite(annual(:,available_stations));

  % Choose the longest match (if equal then the last maximum)

  [overlap_years,overlaps_loc]=sort(sum(overlaps));

  if overlap_years(end)<min_overlap,

      % Not long enough

      fprintf('Station %d has too short overlap (%d years).\n',...

               available_stations,overlap_years');

      available_stations=[];

  else

      match=available_stations(overlaps_loc(end));

      match_overlap=overlaps(:,overlaps_loc(end));

      % calculate bias

      bias=mean(annual(match_overlap,match))-mean(ref_annual(match_overlap));

      % remove bias and add to the fixed stations

      no_fixed=no_fixed+1;

      fixed_stations(:,no_fixed)=temperatures(:,match)-bias;

      % calculate new reference station

      ref_station=nanmean(fixed_stations')';

      % remove from the station list

      available_stations=setdiff(available_stations,match);

  end

end

combined=round(ref_station);

