clear all; close all;

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

   name_title='AVISO';
   name='AVISO';
   day_freq='7d';

   ylimA1=[0 10]; xA1=50; yA1=0.5;                               % plot1
   ylimA2=[0 0.015]; xA2=50; yA2=0.001;                          % plot2
   caxisA3=[0 10]; xA3=-85; yA3=-33; colormapA3='default';       % plot3

   Xmin=2;
   Xmax=126;
   Ymin=2;
   Ymax=163;

   labelA=['computed with ' day_freq 'file and  2003-2007  mean'];

%%% Calcul des indices journaliers pour les années plottées 
   ys=2003;         
   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'}(:)-360;
   zos_alltime=nz{'Grid_0001'}(:,:,:);
   close(nz);

   zos_alltime=0.01*zos_alltime;  % 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,5:end); rsmt=sum(sum(rtmask));  
%     rmgz=mgz(17:134,30:end);
      rvzos=vzos(Ymin:Ymax,Xmin:Xmax); rtmask=tmask(Ymin:Ymax,Xmin:Xmax); rsmt=sum(sum(rtmask));  
      rmgz=mgz(Ymin:Ymax,Xmin:Xmax);

      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,100*(sqrt(ke)));
   grid on;
   hx=xlabel(xday); set(hx,'FontS',14);
   set(gca,'FontS',14); set (gca,'ylim',ylimA1);
   title(['\sigma\_ssh (cm) ' name_title ' (' dur_plot ')'],'FontSize',18);
   text(xA1,yA1,labelA, 'HorizontalAlignment','left','FontSize',14);
   fignameA=['plot_var_ssh_' name '_' dur_plot '_' day_freq];
   print ('-dpng',fignameA);

%  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 (cm) - ' name_title ' (' dur_plot ')'],'FontSize',18);
%  text(xA1,yA1,labelA, 'HorizontalAlignment','left','FontSize',14);
%  text(50,0.005,labelA, 'HorizontalAlignment','left','FontSize',14) ;
%  fignameA=['plot_var_ssh_' name '_' dur_plot ];
%  print ('-dpng',fignameA) ;


%%%  Plot 2
   figure;
   plot(time,sqrt(sgzos));
   grid on;
   hx=xlabel(xday); set(hx,'FontS',14);
   set(gca,'FontS',14); set (gca,'ylim',ylimA2);
   title(['sqrt(grad ssh) ' name_title ' (' dur_plot ')'], 'FontSize',18);
   text(xA2,yA2,labelA, 'HorizontalAlignment','left','FontSize',14);
   fignameA=['plot_grad_ssh_' name '_' dur_plot '_' day_freq];
   print ('-dpng',fignameA);

%  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) ' name_title ' (' dur_plot ')'], 'FontSize',18);
%  text(50,0.001,labelA, 'HorizontalAlignment','left','FontSize',14) ;
%  fignameA=['plot_grad_ssh_' name '_' dur_plot ];
%  print ('-dpng',fignameA) ;

%%%  Plot 3

   figure;
   pcolor(lon(Xmin:Xmax),lat(Ymin:Ymax),100*squeeze(sqrt(mvz(5,Ymin:Ymax,Xmin:Xmax)))); shading flat;
   set(gca,'FontS',14); set(gca,'PlotBoxAspectRatio',[0.7 1 1]); box on;
   colormap(colormapA3); hcc=colorbar; set(hcc,'FontS',14);
   caxis(caxisA3);
   title(['\sigma\_ssh (cm) ' name_title ' (' dur_plot ')'],'FontSize',18);
   text(xA3,yA3,labelA, 'HorizontalAlignment','center','FontSize',14);
   fignameA=['map_var_ssh_' name '_' dur_plot '_' day_freq];
   print ('-dpng',fignameA);

%quit