clear

nx=2160; ny=nx*13; nz=90;
siz=[nx ny];
dirGrid='../grid/';
hc=readbin([dirGrid 'hFacC.data'],[nx ny]);
IX=find(hc==1);
NX=length(IX);

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

%seaice_init_varia.F:
%           sIceLoad(i,j,bi,bj) = HEFF(i,j,bi,bj)*SEAICE_rhoIce
%     &                         + HSNOW(i,j,bi,bj)*SEAICE_rhoSnow
%	    SSH = Etan + sIceLoad*recip_rhoConst
%STDOUT.0000:
SEAICE_rhoIce = 9.100000000000000E+02;
SEAICE_rhoSnow = 3.300000000000000E+02;
rhoConst = 1.027500000000000E+03;
recip_rhoConst = 1./rhoConst;

flds={'Eta', 'SIheff', 'SIhsnow'};
eta=zeros([nx ny]); heff=zeros([nx ny]); hsnow=zeros([nx ny]);
ssh=zeros([nx ny]);

%output
%yearly: 1 map
ssh_yr_mean=zeros([nx ny]);
ssh_yr_sq_hr=zeros([nx ny]);
ssh_yr_sq_dy=zeros([nx ny]);
%monthly: 12 map
ssh_mn_mean=zeros([nx ny]);
ssh_mn_sq_hr=zeros([nx ny]);
ssh_mn_sq_dy=zeros([nx ny]);
%daily: 365 map
ssh_dy_mean=zeros([nx ny]);
%for tides file [NX TX]
fnout='TIDE_SSH.bin';

%time
t0 = datenum(2020,1,19,21,0,0);        deltaT = 45;
ts1=(datenum(2020,3,1)-t0)*86400/deltaT+3600/deltaT;
ts2=(datenum(2021,3,1)-t0)*86400/deltaT;
ts=ts1:3600/deltaT:ts2; %length(ts)==365*24

yr=2020;
k=0;
for mn=3:14 %2020/03 to 2021/02
	ssh_mn_mean=zeros([nx ny]);
	ssh_mn_sq_hr=zeros([nx ny]);
	ssh_mn_sq_dy=zeros([nx ny]);

	days=datenum(yr,mn+1,1)-datenum(yr,mn,1);
	for dy=1:days
			ssh_dy_mean=zeros([nx ny]);
		for hr=1:24
			k=k+1;
			t=( datenum(yr,mn,dy,hr,0,0)-t0 )*86400/deltaT;

			f=1;
			fld=flds{f}; fn=[pp fld '/' fld '.' myint2str(t,10) '.data'];
			eta = readbin(fn,[nx ny]);
			f=2;
			fld=flds{f}; fn=[pp fld '/' fld '.' myint2str(t,10) '.data'];
			heff = readbin(fn,[nx ny]);
			f=3;
			fld=flds{f}; fn=[pp fld '/' fld '.' myint2str(t,10) '.data'];
			hsnow = readbin(fn,[nx ny]);

			ssh = eta + (heff*SEAICE_rhoIce + hsnow*SEAICE_rhoSnow)*recip_rhoConst;
        		writebin(fnout,ssh(IX), 1,'real*4',k -1);

			ssh_dy_mean = ssh_dy_mean + ssh;
			ssh_yr_mean = ssh_yr_mean + ssh;
			ssh_yr_sq_hr = ssh_yr_sq_hr + ssh.^2;
			ssh_mn_sq_hr = ssh_mn_sq_hr + ssh.^2;
		end %hr			
		%daily mean:
			ssh_dy_mean = ssh_dy_mean/24;
		f_day=datestr(t0 + t*deltaT/86400-1,'yyyymmdd');
		disp(f_day)
		writebin(['ssh_' f_day '.bin'],ssh_dy_mean)

		ssh_yr_sq_dy = ssh_yr_sq_dy + ssh_dy_mean.^2;
		ssh_mn_sq_dy = ssh_mn_sq_dy + ssh_dy_mean.^2;
		ssh_mn_mean  = ssh_mn_mean  + ssh_dy_mean;
	end %dy		
	%monthly mean:
		ssh_mn_mean = ssh_mn_mean/days;
		ssh_mn_sq_dy = ssh_mn_sq_dy/days;
		ssh_mn_sq_hr = ssh_mn_sq_hr/(days*24);
	f_mon=datestr(t0 + t*deltaT/86400-1,'yyyymm');
	writebin(['ssh_' f_mon '.bin'],ssh_mn_mean, 1,'real*4',1 -1)
	writebin(['ssh_' f_mon '.bin'],ssh_mn_sq_hr,1,'real*4',2 -1)
	writebin(['ssh_' f_mon '.bin'],ssh_mn_sq_dy,1,'real*4',3 -1)
end %mn

%yearly mean
ssh_yr_mean = ssh_yr_mean/length(ts);
ssh_yr_sq_hr = ssh_yr_sq_hr/length(ts);
ssh_yr_sq_dy = ssh_yr_sq_dy/365;
%save
writebin('ssh_year.bin',ssh_yr_mean, 1,'real*4',1 -1);
writebin('ssh_year.bin',ssh_yr_sq_hr,1,'real*4',2 -1);
writebin('ssh_year.bin',ssh_yr_sq_dy,1,'real*4',3 -1);