clear

nx=2160; ny=nx*13; nz=90;
nn=6; ss=nz/nn; %6 segments


pp='/nobackupp11/dmenemen/DYAMOND/c1440_llc2160/mit_output/';

pout='../W/';
fld='W';

yr=2020;
%monthly	
disp('monthly mean	')
mn=max(randperm(14,1),3);
f_mon=datestr(datenum(yr,mn,1),'yyyymm');
disp(f_mon)
	k=randperm(nz,1);
	k1=floor(k/ss)+1; k2=mod(k,ss);
	if k2==0; k1=k1-1;k2=ss; end
	ssh_mn_mean1=readbin([pout fld '_' f_mon '_k' myint2str(k1,2) '.bin'],[nx ny], 1,'real*4',(1-1)*ss+k2 -1);
	ssh_mn_mean2=readbin([pout fld '_' f_mon '.bin'],                     [nx ny], 1,'real*4',(1-1)*nz+k  -1);
	disp([k isequal(ssh_mn_mean1,ssh_mn_mean2)])
%monthly mean squared based on hourly output
disp('monthly mean squared based on hourly output')
mn=max(randperm(14,1),3);
f_mon=datestr(datenum(yr,mn,1),'yyyymm');
disp(f_mon)
	k=randperm(nz,1);
	k1=floor(k/ss)+1; k2=mod(k,ss);
	if k2==0; k1=k1-1;k2=ss; end
	ssh_mn_sq_hr1=readbin([pout fld '_' f_mon '_k' myint2str(k1,2) '.bin'],[nx ny], 1,'real*4',(2-1)*ss+k2 -1);
	ssh_mn_sq_hr2=readbin([pout fld '_' f_mon '.bin'],                     [nx ny], 1,'real*4',(2-1)*nz+k  -1);
	disp([k isequal(ssh_mn_sq_hr1,ssh_mn_sq_hr2)])
%monthly mean squared based on daily-mean output
disp('monthly mean squared based on daily-mean output')
mn=max(randperm(14,1),3);
f_mon=datestr(datenum(yr,mn,1),'yyyymm');
disp(f_mon)
	k=randperm(nz,1);
	k1=floor(k/ss)+1; k2=mod(k,ss);
	if k2==0; k1=k1-1;k2=ss; end
	ssh_mn_sq_dy1=readbin([pout fld '_' f_mon '_k' myint2str(k1,2) '.bin'],[nx ny], 1,'real*4',(3-1)*ss+k2 -1);
	ssh_mn_sq_dy2=readbin([pout fld '_' f_mon '.bin'],                     [nx ny], 1,'real*4',(3-1)*nz+k  -1);
	disp([k isequal(ssh_mn_sq_dy1,ssh_mn_sq_dy2)])

%yearly	
disp('yearly mean	')
	k=randperm(nz,1);
	k1=floor(k/ss)+1; k2=mod(k,ss);
	if k2==0; k1=k1-1;k2=ss; end
	ssh_mn_mean1=readbin([pout fld '_year_k' myint2str(k1,2) '.bin'],[nx ny], 1,'real*4',(1-1)*ss+k2 -1);
	ssh_mn_mean2=readbin([pout fld '_year.bin'],                     [nx ny], 1,'real*4',(1-1)*nz+k  -1);
	disp([k isequal(ssh_mn_mean1,ssh_mn_mean2)])
%yearly mean squared based on hourly output
disp('yearly mean squared based on hourly output')
	k=randperm(nz,1);
	k1=floor(k/ss)+1; k2=mod(k,ss);
	if k2==0; k1=k1-1;k2=ss; end
	ssh_mn_sq_hr1=readbin([pout fld '_year_k' myint2str(k1,2) '.bin'],[nx ny], 1,'real*4',(2-1)*ss+k2 -1);
	ssh_mn_sq_hr2=readbin([pout fld '_year.bin'],                     [nx ny], 1,'real*4',(2-1)*nz+k  -1);
	disp([k isequal(ssh_mn_sq_hr1,ssh_mn_sq_hr2)])
%yearly mean squared based on daily-mean output
disp('yearly mean squared based on daily-mean output')
	k=randperm(nz,1);
	k1=floor(k/ss)+1; k2=mod(k,ss);
	if k2==0; k1=k1-1;k2=ss; end
	ssh_mn_sq_dy1=readbin([pout fld '_year_k' myint2str(k1,2) '.bin'],[nx ny], 1,'real*4',(3-1)*ss+k2 -1);
	ssh_mn_sq_dy2=readbin([pout fld '_year.bin'],                     [nx ny], 1,'real*4',(3-1)*nz+k  -1);
	disp([k isequal(ssh_mn_sq_dy1,ssh_mn_sq_dy2)])

%daily
disp('daily')
for i=1:20
dy=randperm(365,1);
f_day=datestr(datenum(yr,3,dy),'yyyymmdd');
disp(f_day)
	k=randperm(nz,1);
	k1=floor(k/ss)+1; k2=mod(k,ss);
	if k2==0; k1=k1-1;k2=ss; end
	ssh_dy_mean1=readbin([pout fld '_' f_day '_k' myint2str(k1,2) '.bin'],[nx ny], 1,'real*4',k2 -1);
	ssh_dy_mean2=readbin([pout fld '_' f_day '.bin'],                     [nx ny], 1,'real*4',k  -1);
	disp([i k isequal(ssh_dy_mean1,ssh_dy_mean2)])
end