%exploring SST differences (with and without eddies)

%some ini
 addpath('funct')
if exist('../data/processed/both/Dwpow.mat')~=2
  T=[]; n=0;
 D=dir('/Volumes/scratch/baedw/data/mod/noed/out/forcing/*.nc');

 for i=1:length(D)
  i
  %load data
   wind_power_u=ncread(['/Volumes/scratch/baedw/data/mod/ref/out/forcing/' D(i).name],'wind_power_u');
   wind_power_v=ncread(['/Volumes/scratch/baedw/data/mod/ref/out/forcing/' D(i).name],'wind_power_v');
   ref=wind_power_u+wind_power_v; ref=permute(ref,[3 2 1]);
   wind_power_u=ncread(['/Volumes/scratch/baedw/data/mod/noed/out/forcing/' D(i).name],'wind_power_u');
   wind_power_v=ncread(['/Volumes/scratch/baedw/data/mod/noed/out/forcing/' D(i).name],'wind_power_v');
   noed=wind_power_u+wind_power_v; noed=permute(noed,[3 2 1]);
  %axis information
   xt_ocean=ncread(['/Volumes/scratch/baedw/data/mod/noed/out/forcing/' D(i).name],'xt_ocean');
   yt_ocean=ncread(['/Volumes/scratch/baedw/data/mod/noed/out/forcing/' D(i).name],'yt_ocean');
   time=ncread(['/Volumes/scratch/baedw/data/mod/noed/out/forcing/' D(i).name],'time');
   t = get_calendar_date((time./24)  + get_julian_day(1987,1,1,0));
   T = [T; t];
   x=xt_ocean; clear xt_ocean;
   y=yt_ocean; clear yt_ocean;
   I=find(ref==min(ref(:))); ref(I)=NaN; noed(I)=NaN;
  %calculations
   dmld = ref-noed;;
   for l=1:length(t)
    n=n+1;;
    tmp = dmld(l,:,:);   DMLDa(n)   = meannan(tmp(:));
    tmp = ref(l,:,:);     REFa(n)   = meannan(tmp(:));
    tmp = noed(l,:,:);    NOEDa(n)  = meannan(tmp(:));
   end
  %some late ini
    if i==1
     DMLD = zeros(12,length(y),length(x)); REF=DMLD; NOED=DMLD;
    end
   %more calculations
    mo=t(1,2);
    DMLD(mo,:,:)  = DMLD(mo,:,:)  + mean(dmld(:,:,:));
    REF(mo,:,:)   = REF(mo,:,:)   + mean(ref(:,:,:));
    NOED(mo,:,:)  = NOED(mo,:,:)  + mean(noed(:,:,:));
 end

   %final calulations
    DMLD=DMLD./(length(D)./12);
    REF=REF./(length(D)./12);
    NOED=NOED./(length(D)./12);

save ../data/processed/both/Dwpow.mat DMLD REF NOED DMLDa REFa NOEDa x y t
else
load ../data/processed/both/Dwpow.mat
load ../data/mod/grid/area.mat area_u xu_ocean yu_ocean
end

x=xu_ocean;
y=yu_ocean;

 Pref =DMLD+NaN;
 Pnoed=DMLD+NaN;
 dP   = DMLD+NaN;

%converting W ==> W/m^2
 for i=1:12
  dP(i,:,:)=squee(DMLD(i,:,:))./area_u;
  Pref(i,:,:)=squee(REF(i,:,:))./area_u;
  Pnoed(i,:,:)=squee(NOED(i,:,:))./area_u;
 end

%keyboard
  draw_diffp(x,y,mean(dP(:,:,:)*1000),[-4:0.1:4],' ',[12.5 max(x(:))],[min(y(:)) max(y(:))],2)
  print -depsc -r300 ../ps/dWPOW_annu.eps
  print -dtiff -r500 ../ps/dWPOW_annu.tiff

%display('plot')
%  figure(1)
%   subplot(2,1,1); set(gca,'box','on','linewidth',1.2,'tickdir','out','fontsize',12, ...
%                           'xgrid','on','ygrid','on','nextplot','add','xlim',[1 length(DMLDa)])
%   plot(DMLDa,'k','linewidth',1.5)
%   xlabel('days')
%   ylabel('\Delta Power [W/gridbox]')
%   print -depsc ../ps/dWPOW_ancy.eps
%
%  draw_diffsst(x,y,DMLD(7,:,:),[-20000:2000:20000],'July',[12.5 max(x(:))],[min(y(:)) max(y(:))],0)
%  print -depsc -r300 ../ps/dWPOW_jul.eps
%
%
%  draw_diffsst(x,y,DMLD(3,:,:),[-20000:2000:20000],'November',[12.5 max(x(:))],[min(y(:)) max(y(:))],0)
%  print -depsc -r300 ../ps/dWPOW_nov.eps



