clear all; close all;

   filezos='/loceanfs/pulsation/cholod/data/AVISO/dt_upd_global_merged_msla_h_qd_7d_20030101_20080102_PERU.nc';

   namerun_title='AVISO';
   namerun='AVISO';

   label1='computed with 7d file and  2005-2007  mean';

%%% Calcul des indices journaliers pour les années plottées 
   ys=2005;         
   ye=2007;         
   yy=[ys:1:ye];

   xday='7 days';
   ds=((yy(1)-ys)*52)+1;
   de=(yy(end)-ys+1)*52;
   dur_plot=[num2str(ys) '-' num2str(ye)];      

%%% Lecture zos
   disp(['lecture zos dans ' filezos ' ...']);
   nz=netcdf(filezos,'r'); 
   ntime=length(nz('time'));
   lat=nz{'NbLatitudes'}(:);
   lon=nz{'NbLongitudes'}(:);
   zos_alltime=nz{'Grid_0001'}(:,:,:);
   close(nz);

   zos_alltime=zos_alltime./100;  % aviso en centimetres
   mzos=squeeze(mean(zos_alltime(:,:,:)));
   tmask=mzos;
   tmask(tmask<999)=1;
   tmask(tmask>999)=0;
   mzos=mzos.*tmask;

%break
%pcolor(squeeze(zos_alltime(1,:,:)).*tmask); shading flat;colorbar;

%%% Init
   ke=[];    % declare que vecteur vide
   time=[];
   sgzos=[];
   count=zeros(5,1);
   mvz=zeros(5,size(tmask,1),size(tmask,2));

%%% Calcul var ssh
   display('jour: 1  ');
   ll=0;
   for ii=ds:1:de
      ll=ll+1;
      if mod(ll,10)==0; display(['jour: ' num2str(ll) ' sur ' num2str(de-ds)]); end
      zos=squeeze(zos_alltime(ii,:,:)).*tmask;

      day=((ii/52)-floor(ii/52))*52;
      if day<=13 
%        mzos=squeeze(nm{'zos'}(1,:,:)); 
         seas=1;
      elseif day<=26  & day>=14
%        mzos=squeeze(nm{'zos'}(2,:,:));
         seas=2;
      elseif day<=39  & day>=27
%        mzos=squeeze(nm{'zos'}(3,:,:));
         seas=3;
      elseif day<=52  & day>=40
%        mzos=squeeze(nm{'zos'}(4,:,:));
         seas=4;
      end

      vzos=(zos-mzos).^2;
      [gzx,gzy]=gradient(zos-mzos);
      mgz=(gzx.^2+gzy.^2).*tmask;
      vzos=vzos.*tmask;       %  . pour tout calcul matriciel

      rvzos=vzos(17:134,30:end); rtmask=tmask(17:134,30:end); rsmt=sum(sum(rtmask));  
      rmgz=mgz(17:134,30:end);

      ke(ll)=(sum(sum(rvzos))./rsmt);
      sgzos(ll)=(sum(sum(rmgz))./rsmt);
      time(ll)=0.5+(ll-1)*1;
      mvz(seas,:,:)=squeeze(mvz(seas,:,:))+vzos;
      count(seas)=count(seas)+1;
      mvz(5,:,:)=squeeze(mvz(5,:,:))+vzos;
      count(5)=count(5)+1;
   end


   for jj=1:5
      mvz(jj,:,:)=mvz(jj,:,:)/count(jj);
   end

%%% save datas
%  save(filemat,'ke','sgzos','time','mvz','lon','lat');



   opengl neverselect

%%%  Plot 1
   figure;
   plot(time,sqrt(ke));
   grid on;
   hx=xlabel(xday); set(hx,'FontS',14);
   set(gca,'FontS',14); set (gca,'ylim',[0.00 0.08]);
   title(['\sigma\_ssh (m) - ' namerun_title ' (' dur_plot ')'],'FontSize',18);
   text(50,0.005,label1, 'HorizontalAlignment','left','FontSize',14) 
   figname1=['var_ssh_' namerun '_' dur_plot ]
   print ('-dpng',figname1) 


%%%  Plot 2
   figure;
   plot(time,sqrt(sgzos));
   grid on;
   hx=xlabel(xday); set(hx,'FontS',14);
   set(gca,'FontS',14); set (gca,'ylim',[0.000 0.015]);  
   title(['sqrt(grad ssh) serie - ' namerun_title ' (' dur_plot ')'], 'FontSize',18);
   text(50,0.001,label1, 'HorizontalAlignment','left','FontSize',14) 
   figname1=['grad_ssh_' namerun '_' dur_plot ]
   print ('-dpng',figname1) 

%%%  Plot 3
   figure;

   contourf(lon(30:end),lat(17:134),squeeze(sqrt(mvz(5,17:134,30:end))),40);

   shading flat;   %  enlève les lignes en chaque point
   set(gca,'FontS',14); set(gca,'PlotBoxAspectRatio',[0.7 1 1]); box on;
   hcc=colorbar; set(hcc,'FontS',14);
   caxis ([0.00 0.07])
   set(gca,'XTick',[265:5:290])                               
   set(gca,'XTickLabel',['-95';'-90';'-85';'-80';'-75';'-70'])
   title(['\sigma\_ssh (m) - ' namerun_title ' (' dur_plot ')'],'FontSize',18);
   text(280,-30,label1, 'HorizontalAlignment','center','FontSize',14) 
   figname1=['map_var_ssh_' namerun '_' dur_plot]
   print ('-dpng',figname1) 


quit