%exploring SST differences (with and without eddies)

if exist('../data/processed/both/sst.mat')~=2
%some ini
 addpath('funct')
 T=[]; n=0; YR=[1989:1993];
for yr = YR 
 display(['loading     ' num2str(yr) ' reference sst'])
 ncload(['/Volumes/scratch/baedw/data/processed/ref/sst_' num2str(yr) '.nc']);
 ref  = squee(temp);
  display(['loading     ' num2str(yr) ' noed sst'])
 ncload(['/Volumes/scratch/baedw/data/processed/noed/sst_' num2str(yr) '.nc']);
 noed = squee(temp); 
 if length(ref(:,1,1)) ~= length(noed(:,1,1))
  display('severe ERROR; timeaxis do ot comply')
  break
 end
   display(['calculating ' num2str(yr)])
 
 %axis information
  t = get_calendar_date((time./24)  + get_julian_day(1987,1,1,0));
  x=xt_ocean; clear xt_ocean;
  y=yt_ocean; clear yt_ocean;
  I=find(ref<-1e10); ref(I)=NaN; noed(I)=NaN;
  
 %calculate sst difference               
  dtemp = ref-noed;; 
 
 %calculate areal means (inexact since not weihghted by grid box aerea)
   for i=1:length(t)
    n=n+1;;
    tmp = dtemp(i,:,:);   DTa(n)   = meannan(tmp(:));
    tmp = ref(i,:,:);     REFa(n)  = meannan(tmp(:));
    tmp = noed(i,:,:);    NOEDa(n) = meannan(tmp(:));
   end
   T=[T; t];
 %some late ini
  if yr==YR(1)
   DT = zeros(12,length(y),length(x)); REF=DT; NOED=DT;
  end
 
 %calculate monthly means
  for k=1:12
   I=find(t(:,2)==k);
   DT(k,:,:)  =DT(k,:,:)  + mean(dtemp(I,:,:));
   REF(k,:,:) =REF(k,:,:) + mean(ref(I,:,:));
   NOED(k,:,:)=NOED(k,:,:)+ mean(noed(I,:,:));
  end
end

DT=DT./length(YR);
REF=REF./length(YR);
NOED=NOED./length(YR); 

display('save')
save ../data/processed/both/sst.mat x y T DT REF NOED DTa REFa NOEDa
else
display('loading already-processed data')
load ../data/processed/both/sst.mat
end
display('plot')
  figure(1)
   sdate=datenum(T(:,1),T(:,2),T(:,3),0,0,0);

   subplot(2,1,1); set(gca,'box','on','linewidth',1.2,'tickdir','out','fontsize',18, ...
                           'xgrid','on','ygrid','on','nextplot','add','ylim',[-0.1 0.2],'xlim', ...
                           [sdate(1) sdate(end)])
   plot(sdate,DTa,'k','linewidth',1.5)
   set(gca,'xtick',[726469      726834      727199      727564      727930      728295]);
   datetick('x',28,'keeplimits','keepticks')   
   ylabel('\Delta SST [K]')
   print -depsc ../ps/dSST_ancy.eps

% climatological monthly sst differences (REF - noED) + => eddies warm || - => eddies cool
  draw_diffsst(x,y,DT(7,:,:),[-0.5:0.05:0.5],'July',[12.5 max(x(:))],[min(y(:)) max(y(:))],2)
  print -depsc -r300 ../ps/dSST_jul.eps
  print -dtiff -r500 ../ps/dSST_jul.tiff

  draw_diffsst(x,y,DT(11,:,:),[-0.5:0.05:0.5],'November',[12.5 max(x(:))],[min(y(:)) max(y(:))],2)
  print -depsc -r300 ../ps/dSST_nov.eps
  print -dtiff -r500 ../ps/dSST_nov.tiff


