clear

nx=2160; ny=nx*13; nz=90;

%fn='../SALT/OLD/Salt_202003.bin';
%fn='../SALT/OLD/Salt_year.bin';
%fn='../THETA/Theta_202003.bin';
%fn='../THETA/Theta_year.bin';
%fn='../U/U_202003.bin';
%fn='../U/U_year.bin';
%fn='../V/V_202003.bin';
%fn='../V/V_year.bin';
fn='../W/W_202003.bin';
%fn='../W/W_year.bin';

%SALTanom:
%fn='../SALT/Salt_202003.bin';
%fn='../SALT/Salt_year.bin';

disp(fn)
for k=1:5:nz
hc=readbin('../../grid/hFacC.data',[nx ny], 1,'real*4', k-1);
%hc=readbin('../../grid/hFacW.data',[nx ny], 1,'real*4', k-1);
%hc=readbin('../../grid/hFacS.data',[nx ny], 1,'real*4', k-1);
ii=find(hc>0);

ssh_mn_mean=readbin(fn,[nx ny], 1,'real*4',(1-1)*nz+ k-1);
ssh_mn_sq_hr=readbin(fn,[nx ny],1,'real*4',(2-1)*nz+ k-1);
ssh_mn_sq_dy=readbin(fn,[nx ny],1,'real*4',(3-1)*nz+ k-1);
std_hr = zeros([nx ny]);
std_dy = zeros([nx ny]);

std_hr = ssh_mn_sq_hr - ssh_mn_mean.^2;
std_dy = ssh_mn_sq_dy - ssh_mn_mean.^2;
%SALTanom:
%std_hr(ii) = ssh_mn_sq_hr(ii) - (ssh_mn_mean(ii)-35).^2;
%std_dy(ii) = ssh_mn_sq_dy(ii) - (ssh_mn_mean(ii)-35).^2;

ihr=find(std_hr<0);     idy=find(std_dy<0);
lih = ismember(ihr,ii); lid = ismember(idy,ii);

%disp([k length(ii) length(ihr) sum(lih) length(idy) sum(lid)])
%disp(sprintf('%d %d %d %d %d %d %d\n',[k length(ii) length(ihr) sum(lih) length(idy) sum(lid)]))
fprintf('%d %f %f\n',[k sum(lih)/length(ii) sum(lid)/length(ii)])

end