function [r,r25pct,r75pct]=meanCorr(modeldata,obsdata)
%
% Calculates spatial mean of temporal grid cell correlations from model and 
% observation data with bootstrap confidence intervals (25th and 75th percentiles).
% If the bootstrap uncertainty intervals of two correlation estimates do not 
% overlap then the correlations can be considered statistically separated, 
% with approximately 5% random probability of obtaining no overlap. 
%
% The dimensions of the model data are either [nx,ny,nt] or [nx,ny,nt,ne]. 
% The model data can optionally have an ensemble dimension ne. If present 
% it is assume to be the forth dimension and the algorithm will then account 
% for uncertainty in the ensemble-mean from the limited ensemble size. 
% If not present the algorithm will only account for uncertainty related to temporal sampling. 
% This script does not account for effects of serial correlation along the time 
% dimension. 
%
% For the observations only space and time dimensions i.e. [nx,ny,nt] are supported.
%
% Example with x,y dimensions=10 and time dimension=30:
%   nx=10; ny=10; nt=30; 
%   modeldata = rand(nx,ny,nt,'single'); 
%   obsdata = modeldata + rand(nx,ny,nt,'single'); 
%   [r,r25pct,r75pct]=meanCorr(modeldata,obsdata)
%

% check dimensions
[nx,ny,nt]=size(obsdata); % set time dimension 
if length(size(modeldata)) == 3
  ne=1; % modeldata has no ensemble dimension
elseif length(size(modeldata)) == 4 
  ne=size(modeldata,4); % modeldata has ensemble dimension, assume it is the third dimension
else 
  display('modeldata has wrong number of dimensions')
  return
end 

nboot=10000; % resample this often 
rboot=zeros(nboot,1); % bootstrap correlations
for n=1:nboot
  % resample along time dimension
  tind=floor(rand(nt,1)*nt)+1; 
  if ne>1 
    % resample along ensemble dimension 
    eind=floor(rand(ne,1)*ne)+1; 
    % re-compute ensemble mean of resampled data and resample along time dimension
    a=squeeze(nanmean(modeldata(:,:,tind,eind),4)); 
  else
    % no ensemble dimension, resample along time dimension
    a=modeldata(:,:,tind); 
  end
  % assume obs data has no ensemble dimension 
  b=obsdata(:,:,tind); 
  % make sure that missing values are the same for model and obs data 
  a(isnan(b))=nan;
  b(isnan(a))=nan;
  % subtract time mean 
  a=a-nanmean(a,3); 
  b=b-nanmean(b,3); 
  % compute correlation of resampled data in each grid cell 
  r2d=nansum(a.*b,3)./(sqrt(nansum(a.*a,3)).*sqrt(nansum(b.*b,3)));
  % compute average of correlation values over all grid cells
  rboot(n)=nanmean(r2d(:));  
end 
r25pct=prctile(rboot,25);
r75pct=prctile(rboot,75);

% compute correlation 
if ne>1 
  a=squeeze(mean(modeldata,4));
else
  a=modeldata;
end 
b=obsdata; 
a(isnan(b))=nan;
b(isnan(a))=nan;
a=a-nanmean(a,3);
b=b-nanmean(b,3);
r2d=nansum(a.*b,3)./(sqrt(nansum(a.*a,3)).*sqrt(nansum(b.*b,3)));
r=double(nanmean(r2d(:)));
