clear all; close all;
 
   replot='FALS';
%  replot='TRUE';

%%% periode du fichier pluriannuel
   ys_fich=2003;         
   ye_fich=2007;         
%%% periode du plot temporel
   ys_plot=2003;         
   ye_plot=2007;         

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

   exp(1).config='peru12_nemo';
   exp(1).names={'erai70b' 'erai74c' 'p12n_bdy_tr12n_quik' 'p12n_bdy_tr12n_quik_ND' 'p12n_bdy_tr12n_quik_blp1' 'p12n_bdy_soda2' 'p12n_bdy_soda2_B' 'p12n_bdy_soda2_B_rfact' 'p12n_bdy_soda2_B_steps'};
   exp(1).day_freq='5d';
   exp(1).grid_type='T_2D';

   exp(2).config='trop12_nemo/peru';
   exp(2).names={'tr12_quik'};
   exp(2).day_freq='5d';
   exp(2).grid_type='T_2D';

   exp(3).config='trop025peru12_nemo/peru';
   exp(3).names={'tr025p12_quik'};
   exp(3).day_freq='5d';
   exp(3).grid_type='T_2D';

   exp(4).config='/loceanfs/pulsation/cholod/data/AVISO/';
   exp(4).names={'AVISO_PERU12'};
   exp(4).day_freq='7d';
   exp(4).grid_type='T_2D';

%  day_freq='5d';
%  grid_type='T_2D';

%  les 2 masks sont identiques sauf sur les caraibes (masquées ou pas) => masquées pour 
%  p12n_bdy_soda2_B_rfact' et  'p12n_bdy_soda2_B_steps' pour lesquels les caraibes ne sont pas masquées     
%  gridfile='/loceanfs/pulsation/cholod/peru12_nemo/inputs/grd/mesh_mask_peru12.L75_from_trop12.nc';
   gridfile='/loceanfs/pulsation/cholod/peru12_nemo/inputs/grd/mesh_mask_peru12.L75_bathy_Stephane.nc';
   exp_dim=size(exp);

%  for ii = 1:exp_dim(2)
   for ii = 1:1 
      day_freq=exp(ii).day_freq;
      grid_type=exp(ii).grid_type;
      config=exp(ii).config;
      names=exp(ii).names;
      names_dim=size(names);
      disp([num2str(ii) ' : ' gridfile ' : ' config ' : ' day_freq ' : ' grid_type]);

%%% Calcul des indices journaliers pour les années plottées 
      yy=[ys_plot:1:ye_plot];
      if day_freq=='5d'
         ds=((yy(1)-ys_fich)*73)+1;
         de=(yy(end)-ys_fich+1)*73;
         xday='5 days';
      elseif day_freq=='1d'
         ds=((yy(1)-ys_fich)*365)+1;
         de=(yy(end)-ys_fich+1)*365;
         xday='1 day';
      elseif day_freq=='7d'
         ds=((yy(1)-ys_fich)*52)+1;
         de=(yy(end)-ys_fich+1)*52;
         xday='7 days';
      end

%     for jj = 1:names_dim(2);
      for jj = 7:9
         name=names{jj};
         disp([num2str(jj) ' : ' name]);
         name_title=regexprep(name,'_','\\_');
         rep_fich=[num2str(ys_fich) '0101_' num2str(ye_fich) '1231'];
         rep_plot=[num2str(ys_plot) '0101_' num2str(ye_plot) '1231'];
         dur_plot=[num2str(ys_plot) '-' num2str(ye_plot)];
         if strcmp(name,'AVISO_PERU12')
            dir=config;
            labelA=['computed with 7days file and ' dur_plot ' mean'];
         else
            dir=['/loceanfs/pulsation/cholod/' config '/outputs/' name '/' rep_fich '/'];
            labelA=['computed with ' day_freq ' file and ' dur_plot ' mean'];
         end

         if replot == 'TRUE'
            filemat=[dir 'var_ssh_' name '_' rep_plot '.mat'];
            load(filemat);
         else
            if strcmp(name,'AVISO_PERU12')
               filezos=[dir name '_' dur_plot '_7d.nc']
            else
               filezos=[dir name '_' day_freq '_' rep_fich '_grid_' grid_type '.nc'];
            end

%%% Mask et coord de la grille
            ng=netcdf(gridfile,'r');
            tmask=squeeze(ng{'tmaskutil'}(1,:,:));
            lat=ng{'nav_lat'}(:);
            lon=ng{'nav_lon'}(:);
            close(ng);

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

%%% Lecture zos
            disp(['lecture zos dans ' filezos ' ...']);
            nz=netcdf(filezos,'r'); 
            ntime=length(nz('time_counter'));
            zos_alltime=nz{'zos'}(:,:,1:361);
            if strcmp(name,'AVISO_PERU12')
               zos_alltime=0.01*zos_alltime;
            end
            mzos=squeeze(mean(zos_alltime(ds:de,:,:)));
            close(nz);

%%% Calcul var ssh
            display(['extraction dans le fichier ' num2str(ys_fich) '-' num2str(ye_fich) ' des ' num2str(de-ds) ' pdts de ' num2str(ds) ' a ' num2str(de) ]);
%  display('pdt: 1  ');
            display(['pdt:  ' num2str(ds)]);
            ll=0;
            for ii=ds:1:de
               ll=ll+1;
%              if mod(ll,10)==0; display(['jour: ' num2str(ll) ' sur ' num2str(de-ds)]); end
               if mod(ll,10)==0; display(['pdt: ' num2str(ii) ' sur ' num2str(de)]); end
               zos=squeeze(zos_alltime(ii,:,:));

               if day_freq=='5d'
                  day=((ii/73)-floor(ii/73))*73;
                  if day<=18 
%                    mzos=squeeze(nm{'zos'}(1,:,:)); 
                     seas=1;
                  elseif day<=36  & day>=19
%                    mzos=squeeze(nm{'zos'}(2,:,:));
                     seas=2;
                  elseif day<=54  & day>=37
%                    mzos=squeeze(nm{'zos'}(3,:,:));
                     seas=3;
                  elseif day<=73  & day>=55
%                    mzos=squeeze(nm{'zos'}(4,:,:));
                     seas=4;
                  end
               elseif day_freq=='1d'
                  day=((ii/365)-floor(ii/365))*365;
                  if day<=90 
%                   mzos=squeeze(nm{'zos'}(1,:,:)); 
                    seas=1;
                  elseif day<=180  & day>=91
%                   mzos=squeeze(nm{'zos'}(2,:,:));
                    seas=2;
                  elseif day<=270  & day>=181
%                   mzos=squeeze(nm{'zos'}(3,:,:));
                    seas=3;
                  elseif day<=365  & day>=271
%                   mzos=squeeze(nm{'zos'}(4,:,:));
                    seas=4;
                  end
               elseif day_freq=='7d'
                  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
               end

               vzos=(zos-mzos).^2;
               [gzx,gzy]=gradient(zos-mzos);
               mgz=(gzx.^2+gzy.^2).*tmask;
               vzos=vzos.*tmask;       %  . pour tout calcul matriciel
      
%              if resolution=='012'
%                 rvzos=vzos(50:400,70:end); rtmask=tmask(50:400,70:end); rsmt=sum(sum(rtmask));  
%                 rmgz=mgz(50:400,70:end);
                  rvzos=vzos(:,:); rtmask=tmask(:,:); rsmt=sum(sum(rtmask));  
                  rmgz=mgz(:,:);
%              elseif resolution=='025'
%                 rvzos=vzos(17:134,23:end); rtmask=tmask(17:134,23:end); rsmt=sum(sum(rtmask));  
%                 rmgz=mgz(17:134,23:end);
%              elseif resolution=='075'
%                 rvzos=vzos(5:45,8:end); rtmask=tmask(5:45,8:end); rsmt=sum(sum(rtmask));  
%                 rmgz=mgz(5:45,8:end);
%              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
            filemat=[dir 'var_ssh_' name '_' rep_plot '.mat'];
            save(filemat,'ke','sgzos','time','mvz','lon','lat');
         end

         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); set (gca,'xlim',[ds de]);
         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 '_ZOOM'];
         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); set (gca,'xlim',[ds de]);
         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 '_ZOOM'];
         print ('-dpng',fignameA); 

%%%  Plot 3
         figure;
%        pcolor(lon(:,1:361),lat(:,1:361),100*squeeze(sqrt(mvz(5,:,1:361))));
         pcolor(lon(50:400,70:end),lat(50:400,70:end),100*squeeze(sqrt(mvz(5,50:400,70:end))));
%        if resolution=='012'
%           contourf(lon(50:400,70:end),lat(50:400,70:end),squeeze(sqrt(mvz(5,50:400,70:end))),40);
%           contourf(lon(50:400,70:end),lat(50:400,70:end),squeeze(sqrt(mvz(5,50:400,70:end))),40);
%        elseif resolution=='025'
%           contourf(lon(17:134,23:end),lat(17:134,23:end),squeeze(sqrt(mvz(5,17:134,23:end))),40);
%        elseif resolution=='075'
%           contourf(lon(5:45,8:end),lat(5:45,8:end),squeeze(sqrt(mvz(5,5:45,8:end))),40);
%        end
         shading flat;   %  enlève les lignes en chaque point
         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 (m) ' name_title ' (' dur_plot ')'],'FontSize',18);
         text(xA3,yA3,labelA, 'HorizontalAlignment','center','FontSize',14);
         fignameA=['map_var_ssh_' name '_' dur_plot '_' day_freq '_ZOOM'];
         print ('-dpng',fignameA); 

     end
   end

%quit