clear

%nsidc
% read northern data grid coordinates
nnx=304; nny=448;
nlat=readbin('psn25lats.dat',[nnx nny],1,'int32')/100000;
nlon=readbin('psn25lons.dat',[nnx nny],1,'int32')/100000;
ix=find(nlon>180);
nlon(ix)=nlon(ix)-360;
% read southern data grid coordinates
snx=316; sny=332;
slat=readbin('pss25lats.dat',[snx sny],1,'int32')/100000;
slon=readbin('pss25lons.dat',[snx sny],1,'int32')/100000;
ix=find(slon>180);
slon(ix)=slon(ix)-360;

nlon2=zeros([nnx+1 nny]); nlat2=zeros([nnx nny+1]);
slon2=zeros([snx+1 sny]); slat2=zeros([snx sny+1]);
nlon2(1:nnx,:)=nlon; nlon2(end,:)=nlon(end,:);
slon2(1:snx,:)=slon; slon2(end,:)=slon(end,:);
nlat2(:,1:nny)=nlat; nlat2(:,end)=nlat(:,end);
slat2(:,1:sny)=slat; slat2(:,end)=slat(:,end);
dnlon=nlon2(2:end,:)-nlon2(1:end-1,:);
dslon=slon2(2:end,:)-slon2(1:end-1,:);
dnlat=nlat2(:,2:end)-nlat2(:,1:end-1);
dslat=slat2(:,2:end)-slat2(:,1:end-1);

%{
domain='n';
fn=['raw/bt_20200301_f17_v3.1_' domain '.bin'];
nic=readbin(fn,[nnx nny],1,'int16',0,'l')/1000;
nic(nic==1.2)=nan; nic(nic==1.1)=1;
domain='s';
fn=['raw/bt_20200301_f17_v3.1_' domain '.bin'];
sic=readbin(fn,[snx sny],1,'int16',0,'l')/1000;
sic(sic==1.2)=nan;
subplot(121)
m_proj('stereo','lat',90,'lon',0,'rad',40)
m_pcolor(nlon,nlat,nic),shading flat,colorbar
m_grid('xtick',-180:30:180,'ytick',40:10:80)
subplot(122)
m_proj('stereo','lat',-90,'lon',0,'rad',40)
m_pcolor(slon,slat,sic),shading flat,colorbar
m_grid('xtick',-180:30:180,'ytick',-80:10:-40,'Xaxislocation','top','Yaxislocation','middle')
%}	   



%llc2160
nx=2160; ny=nx*13; nz=90;
pp='/nobackupp11/dmenemen/DYAMOND/c1440_llc2160/mit_output/';
hFacC=readbin('../../../grid/hFacC.data',[nx ny]);
xc=readbin('../../../grid/XC.data',[nx ny]);
yc=readbin('../../../grid/YC.data',[nx ny]);
ssh=zeros([nx ny]); ssh24=zeros([nx ny 24]); ssh5=zeros([nx ny 5]); 
fld='SIarea';

%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

% interpolate to NSIDC grid
slan = nan*nlon; slas = nan*slon;
nic_llc=zeros([size(nlon) 365],'single');
sic_llc=zeros([size(slon) 365],'single');
nic_nsidc=zeros([size(nlon) 365],'single');
sic_nsidc=zeros([size(slon) 365],'single');

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

for dy=1:365
tic

slan = nan*nlon; slas = nan*slon;

        f_day=datestr(datenum(yr2012,3,dy),'yyyymmdd');
        disp(f_day)

	domain='n';
        fn=['raw/bt_' f_day '_f17_v3.1_' domain '.bin'];
	nic=readbin(fn,[nnx nny],1,'int16',0,'l')/1000;
	nic(nic==1.2)=nan; nic(nic==1.1)=1;

	domain='s';
        fn=['raw/bt_' f_day '_f17_v3.1_' domain '.bin'];
	sic=readbin(fn,[snx sny],1,'int16',0,'l')/1000;
	sic(sic==1.2)=nan;

%llc2160 here: ssh
	for hr=1:24
            t = ( datenum(yr2020,3,dy,hr,0,0)-t0 )*86400/deltaT;
	    fn = [pp fld '/' fld '.' myint2str(t,10) '.data'];
	    ssh5 = readbin(fn,[nx ny 5]);
	    ssh24(:,:,hr) = sum(ssh5,3);
	end    
	ssh = mean(ssh24,3);

%interp here:
	for i=1:nnx
		mydisp(i)
   	for j=1:nny
		if dy==1
		ix=find( xc>=nlon(i,j)-.1 & xc<=nlon(i,j)+.1 & ...
		         yc>=nlat(i,j)-.1 & yc<=nlat(i,j)+.1 & ...
			 hFacC>0 );
		nix{i,j}=ix;
		else
		ix=nix{i,j};
		end
		if length(ix)>0
       		        slan(i,j)=mean(ssh(ix));
        	end
    	end
	end

	for i=1:snx
		mydisp(i)
   	for j=1:sny
		if dy==1
		ix=find( xc>=slon(i,j)-.1 & xc<=slon(i,j)+.1 & ...
		         yc>=slat(i,j)-.1 & yc<=slat(i,j)+.1 & ...
			 hFacC>0 );
		six{i,j}=ix;			 
		else
		ix=six{i,j};
		end
		if length(ix)>0
       		        slas(i,j)=mean(ssh(ix));
        	end
    	end
	end

        nic_llc(:,:,dy)=single(slan);
        nic_nsidc(:,:,dy)=single(nic);
        sic_llc(:,:,dy)=single(slas);
        sic_nsidc(:,:,dy)=single(sic);
toc
end %dy
	fn=[pout 'NSIDC_LLC_AREA2.mat'];
        save(fn, 'nic_llc','nic_nsidc','sic_llc','sic_nsidc');