%extract bottom temperatures

%define output
 outfil='../../data/processed/botTemp_reference.mat'
%define input
 idir='/Volumes/MOMSO/data/reference/yearly/';
 fils = dir([idir 'ocean_physic_*.nc']);


%get grid once
 x=ncread([idir fils(1).name],'xt_ocean');
 y=ncread([idir fils(1).name],'yt_ocean');
 z=ncread([idir fils(1).name],'st_ocean');
 [X,Y]=meshgrid(x,y); 

%some ini
 n=0;
 btemp = zeros(length(x),length(y),length(fils))+NaN;

for ii=1:length(fils)
 n=n+1;
 
 t=ncread([idir fils(ii).name],'temp'); 
 I=find(t<-1e10); t(I)=NaN; t(I)=NaN;
 t=double(t);

 tim=ncread([idir fils(ii).name],'time');
 timeStr(n,:)=datevec(daynoleap2datenum(tim,1900));


for i=1:length(x)
 for j=1:length(y)
  I=find(isnan(squee(t(i,j,:)))~=1);
  if length(I)>0 
   if z(I(end))<=1500
    btemp(i,j,n)=t(i,j,I(end));
   end
  end
 end
end

end

save(outfil, 'x', 'y', 'timeStr', 'btemp')

