lag correlation

clear

% read the ssh file
fid1=fopen('zeta_year.grd','r')
zeta1=fread(fid1,'float32');
fclose(fid1);
zeta=reshape(zeta1,601,182,50);

for i=1:601
for j=1:182
for t=1:50
if zeta(i,j,t)>1e+10
zeta(i,j,t)=NaN;
end
end
end
end

% select the 4 point
lat=[5:0.25:50.25]';
lon=[110:0.25:260]';

% Point 1
latmn = max(size(lat(lat<=43.5)));
latmx = max(size(lat(lat<=43.5)));
lonmn = max(size(lon(lon<=140.5)));
lonmx = max(size(lon(lon<=140.5)));
p1=zeta(lonmn:lonmx,latmn:latmx,:);
clear latmn latmx lonmn lonmx

% Point 2
latmn = max(size(lat(lat<=37.25)));
latmx = max(size(lat(lat<=37.25)));
lonmn = max(size(lon(lon<=137)));
lonmx = max(size(lon(lon<=137)));
p2=zeta(lonmn:lonmx,latmn:latmx,:);
clear latmn latmx lonmn lonmx

% Point 3
latmn = max(size(lat(lat<=35)));
latmx = max(size(lat(lat<=35)));
lonmn = max(size(lon(lon<=132)));
lonmx = max(size(lon(lon<=132)));
p3=zeta(lonmn:lonmx,latmn:latmx,:);
clear latmn latmx lonmn lonmx

% Point 4
latmn = max(size(lat(lat<=32.5)));
latmx = max(size(lat(lat<=32.5)));
lonmn = max(size(lon(lon<=131.75)));
lonmx = max(size(lon(lon<=131.75)));
p4=zeta(lonmn:lonmx,latmn:latmx,:);
clear latmn latmx lonmn lonmx

% domain mean
for t=1:50
romsdam(t)=nanmean(nanmean(zeta(:,:,t),1),2);
end

% Mean values of 4 points
p4m=(p1+p2+p3+p4)/4;
%p4m=squeeze((p1+p2+p3+p4)/4)-romsdam';

%for i=1:601
%for j=1:182
%for t=1:50

%zeta(i,j,t)=zeta(i,j,t)-romsdam(t);
%end
%end
%end

%p4m=p3;
% Compute 4 years lag-correlation

for lag=40
%lag=0;
N=length(p4m);
nstr=1; % since 1971 year
% standard deviation of p4m time series
sd_p4m=std(p4m(nstr+lag:N),1);

% mean of p4m time series in the responding period (namely, nstr:N-lag)
mean_p4m=mean(p4m(nstr+lag:N));

% standard deviation of ssh map
sd_ssh=std(zeta(:,:,nstr:N-lag),1,3);

% mean of ssh map in the responding period (namely, nstr+lag:N)
mean_ssh=mean(zeta(:,:,nstr:N-lag),3);

summ=0;
for i=nstr:N-lag
summ=summ+(p4m(i+lag)-mean_p4m)*(zeta(:,:,i)-mean_ssh);
end

r(:,:,5-lag)=summ/(N-lag-nstr+1)./(sd_p4m.*sd_ssh);
end

clear lag
lag=1;
% standard deviation of p4m time series
sd_p4m=std(p4m(nstr:N-lag),1);

% mean of p4m time series in the responding period (namely, nstr:N-lag)
mean_p4m=mean(p4m(nstr:N-lag));

% standard deviation of ssh map
sd_ssh=std(zeta(:,:,nstr+lag:N),1,3);

% mean of ssh map in the responding period (namely, nstr+lag:N)
mean_ssh=mean(zeta(:,:,nstr+lag:N),3);

summ=0;
for i=nstr:N-lag
summ=summ+(p4m(i)-mean_p4m)*(zeta(:,:,i+lag)-mean_ssh);
end

r(:,:,6)=summ/(N-lag-nstr+1)./(sd_p4m.*sd_ssh);
addpath(genpath('/home/liu'))
wrtgrd1(r,'lagcorr2original.grd')