%read and cut data

%define input
 fils = dir('../data/processed/ref/O*')

%read data
 for k=1:length(fils)
  fil = ['../data/processed/ref/' fils(k).name];
  display(['loading ' fil])
  if k==1
   x=squeeze(ncread(fil,'XT_OCEAN'));
   y=squeeze(ncread(fil,'YT_OCEAN'));
   t=squeeze(ncread(fil,'TAXIS'));
   M=single(squeeze(ncread(fil,'MOD_SST')));
   O=single(squeeze(ncread(fil,'OBS_SST')));
  else
   L=length(squee(M(1,1,:)));
  
   tmp=single(squeeze(ncread(fil,'MOD_SST')));
   Ls=length(squee(M(1,1,:)))+1; Le=length(squee(tmp(1,1,:)))+Ls-1;
   M(:,:,Ls:Le)=tmp;

   tmp=single(squeeze(ncread(fil,'OBS_SST')));
   O(:,:,Ls:Le)=tmp;
   
   tmp=squeeze(ncread(fil,'TAXIS'));
   t(Ls:Le)=tmp;
  end
 end

%convert time axis
 for i=1:length(t);
  tmp=get_calendar_date(t(i)./24+get_julian_day(1901,1,15,0));
  yr(i)=tmp(1); mo(i)=tmp(2); da(i)=tmp(3);
 end

%cut
 xlims=[25 30]; ylims=[59.4 60.6];
 xlims=[0 100]; ylims=[0 100];

 I=find(x>=xlims(1) & x<=xlims(2)); J=find(y>=ylims(1) & y<=ylims(2));

 M=M(I,J,:); O=O(I,J,:); x=x(I); y=y(J);

%look
 flat_pcolor(x,y,squee(double(M(:,:,2)))')

%exclude timeslices 
 Iob = [237:241 470]; %apparently deficient observations
 Imo = [153:157 418]; %apparently deficient model snaps

 M(:,:,[Iob Imo])=[];
 O(:,:,[Iob Imo])=[]; 
    yr([Iob Imo])=[];
    mo([Iob Imo])=[];
    da([Iob Imo])=[];
     t([Iob Imo])=[];

%save
 save ../data/processed/ref/BSHObsModSSTdat.mat M O x y yr mo da t
