clear

fn = 'woa18_A5B7_s00_04.nc';
lon = ncread(fn,'lon');    
lat = ncread(fn,'lat');  
depth = ncread(fn,'depth');  
s_an = ncread(fn,'s_an');

%fn = 'woa18_A5B7_t00_04.nc';
%t_an = ncread(fn,'t_an');
%[Y X Z]=meshgrid(lat,lon,depth); PRES=pressure(Z,Y);
%temp=insitutemp(s_an,t_an,PRES,0*PRES);
temp=s_an;


lon_bnds = ncread(fn,'lon_bnds');
lat_bnds = ncread(fn,'lat_bnds');
depth_bnds = ncread(fn,'depth_bnds');
%disp(isequal(lat',mean(lat_bnds,1)))
%disp(isequal(lon',mean(lon_bnds,1)))


%llc2160
nx=2160; ny=nx*13; nz=90;
pp='/nobackupp11/dmenemen/DYAMOND/c1440_llc2160/mit_output/STATS/';
hFacC=readbin('../../../grid/hFacC.data',[nx ny nz]);
xc=readbin('../../../grid/XC.data',[nx ny]);
yc=readbin('../../../grid/YC.data',[nx ny]);
load('../../../grid/thk90.mat')
DPT25=dpt90; dpt=depth;

% interpolate to LEVITUS grid
temp2=nan*temp;
salt2=nan*temp;
st = nan(length(lon),length(lat));

lpth='/nobackupp11/dmenemen/DYAMOND/c1440_llc2160/mit_output/STATS/';
fin=[lpth 'SALT/Salt_year.bin'];
disp(fin)
temp1=readbin(fin,[nx ny nz]);
	
for k=1:1
	disp(k)
	hcc=hFacC(:,:,k);
	tmp=temp1(:,:,k);
	for i=1:length(lon)
		mydisp(i)
		ix=find( xc>=lon_bnds(1,i) & xc<=lon_bnds(2,i) & hcc>0 );
		if length(ix)>0
%    		for j=1:length(lat)
		jj=find(~isnan(temp(i,:,k)));
		for j=jj
		        iy=find( yc(ix)>=lat_bnds(1,j) & yc(ix)<=lat_bnds(2,j) );
		        if length(iy)>0
       		        st(i,j)=mean(tmp(ix(iy)));
	        	end
        	end
	    	end
	end
end %dpt
temp2(:,:,k)=st;
st = nan(length(lon),length(lat));

for k=2:length(dpt)
	disp(k)
	kk=closest(dpt(k),DPT25,0);
	tmp=(temp1(:,:,kk(1))*abs(dpt(k)-DPT25(kk(2)))+...
	     temp1(:,:,kk(2))*abs(dpt(k)-DPT25(kk(1))))/...
	     abs(DPT25(kk(2))-DPT25(kk(1)));
%	hcc=(hFacC(:,:,kk(1))*abs(dpt(k)-DPT25(kk(2)))+...
%	     hFacC(:,:,kk(2))*abs(dpt(k)-DPT25(kk(1))))/...
%	     abs(DPT25(kk(2))-DPT25(kk(1)));
	hcc=hFacC(:,:,kk(2));
	for i=1:length(lon)
		mydisp(i)
		ix=find( xc>=lon_bnds(1,i) & xc<=lon_bnds(2,i) & hcc>0 );
		if length(ix)>0
%    		for j=1:length(lat)
		jj=find(~isnan(temp(i,:,k)));
		for j=jj
		        iy=find( yc(ix)>=lat_bnds(1,j) & yc(ix)<=lat_bnds(2,j) );
		        if length(iy)>0
       		        st(i,j)=mean(tmp(ix(iy)));
	        	end
        	end
	    	end
	end
	temp2(:,:,k)=st;
st = nan(length(lon),length(lat));
end %dpt
salt_woa=single(s_an); salt_llc=single(temp2);
save WOA_LLC_salt.mat salt_woa salt_llc