clear

fn='raw/dt_global_allsat_phy_l4_20120101_20210726.nc';
disp(fn)

lon = ncread(fn,'longitude');    
lat = ncread(fn,'latitude');    
lon_bnds = ncread(fn,'lon_bnds');
lat_bnds = ncread(fn,'lat_bnds');
%disp(isequal(lat',mean(lat_bnds,1)))
%disp(isequal(lon',mean(lon_bnds,1)))
sla = ncread(fn,'sla'); %meter

%llc2160
nx=2160; ny=nx*13; nz=90;
pp='/nobackupp11/dmenemen/DYAMOND/c1440_llc2160/mit_output/STATS/';
hFacC=readbin('../../../grid/hFacC.data',[nx ny]);
xc=readbin('../../../grid/XC.data',[nx ny]);
yc=readbin('../../../grid/YC.data',[nx ny]);
ii=find(hFacC>0);
NX=length(ii);
ssh=zeros([nx ny]);
ssh24=zeros([NX 24]); sshhr=zeros([NX 1]);

pout='25km/';
yr2012=2012;
yr2020=2020;

% interpolate to AVISO grid
ssh2 = nan(length(lon),length(lat),366*24,'single');

if ~exist('SSH_hourly_raw.bin')
%interp here:
	for i=1:length(lon)
		mydisp(i)
		ix=find( xc>=lon_bnds(1,i) & xc<=lon_bnds(2,i) & hFacC>0 );
		if length(ix)>0
    		for j=1:length(lat)
		        iy=find( yc(ix)>=lat_bnds(1,j) & yc(ix)<=lat_bnds(2,j) );
		        if length(iy)>0
       		        idx{i,j}=ix(iy);
	        	end
        	end
	    	end
	end
%


for dy=0:365
tic

%llc2160 here: ssh
	ssh=zeros([nx ny]);
	if dy==0
		fn='TIDE_SSH_20200229.bin';
		disp([fn '@ ' myint2str(dy,3)])
		ssh24 = readbin(fn,[NX 24]);
	else		
		fn='TIDE_SSH.bin';
		disp([fn '@ ' myint2str(dy,3)])
		ssh24 = readbin(fn,[NX 24],1,'real*4',dy -1);
	end
	for hr=1:24
	ssh(ii)=ssh24(:,hr);

%interp here:
	for i=1:length(lon)
 	for j=1:length(lat)
       	        ixiy=idx{i,j};
		if length(ixiy)>0
			ssh2(i,j,dy*24+hr)=mean(ssh(ixiy));
		end
       	end
	end

	end %hr
toc
end %dy

%to save raw
writebin('SSH_hourly_raw.bin',ssh2)
else
for dy=0:365
mydisp(dy)
dd=dy*24+(1:24);
ssh2(:,:,dd)=readbin('SSH_hourly_raw.bin',[length(lon),length(lat) 24],1,'real*4',dy);
end
end


%TIDES start
t0=datenum(2020,2,29,1,0,0);
for i=1:length(lon)
for j=1:length(lat)
disp([i j])
	xin = squeeze(double(ssh2(i,j,:)));
	if ~isnan(xin(1))
	lat1=lat(j);
	[~, xout]=t_tide(xin,'interval',1, 'start',t0, 'latitude',lat1, 'output','none');
%DE-tiding
	ssh2(i,j,:) = xin-xout;
	end
end
end
%TIDES end

%to save de-tided
fn='SSH_hourly_detided.bin';
writebin(fn,ssh2)
sshyear=mean(ssh2,3);

for dy=1:365
tic

sla2 = nan(length(lon),length(lat));

        f_day=datestr(datenum(yr2020,3,dy),'yyyymmdd');
        disp(f_day)
	fin=['raw/dt_global_allsat_phy_l4_' f_day '_20210726.nc'];
%in case NRT!!	
	if ~exist(fin)
        f_da2=datestr(datenum(yr2020,3,dy)+6,'yyyymmdd');
	fin=['raw/nrt_global_allsat_phy_l4_' f_day '_' f_da2 '.nc'];
	end
	sla = ncread(fin,'sla'); %meter

	dd=( (dy-1)*24+1:dy*24 )+11; %12pm to 11am
	sshday = mean(ssh2(:,:,dd),3);
	sla2 = sshday - sshyear; %SLA


	fn=[pout 'AVISO_LLC_SLA_' f_day '.mat'];
        sla_llc=single(sla2);
        sla_aviso=single(sla);
        save(fn, 'sla_aviso','sla_llc');
toc
end %dy