clear all; close all;

%  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.L46_bathy_Stephane.nc';

%  resolution='012';
   resolution='012';

%  config='peru12_nemo';
   config='peru12_nemo';

%  namerun='bdy4b';
   namerun='p12n_bdy_soda2';

%  day_freq='5d';
   day_freq='5d';

%  grid_type='T_2D';
   grid_type='T_2D';

%%% periode du fichier pluriannuel
%  ys_fich=2000;         
   ys_fich=2003;         
%  ye_fich=2007;         
   ye_fich=2007;         

%%% periode du plot temporel
%  ys_plot=2003;         
   ys_plot=2003;         
%  ye_plot=2007;         
   ye_plot=2007;         


   namerun_title=regexprep(namerun,'_','\\_');

   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)];

   dir_in=['/loceanfs/pulsation/cholod/' config '/outputs/' namerun '/' rep_fich '/'];

%  if config(1:4) == 'peru'    %  nomenclature chris
%     filezos=[dir_in 'peru12_' day_freq '_' rep_fich '_grid_' grid_type '.' namerun '.nc']; 
%%    filezos=[dir_in 'peru12_' day_freq '_' rep_fich '_grid_' grid_type '_detrend.' namerun '.nc']; 
%  elseif config(1:4) == 'trop'    %  nomenclature Seb
      filezos=[dir_in namerun '_' day_freq '_' rep_fich '_grid_' grid_type '.nc']; 
%  end
%_________________

   filemat=['var_ssh_' namerun '_' rep_plot '.mat'];

   label1=['computed with ' day_freq ' file and ' dur_plot ' mean'];


%%% 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';
   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'}(:,:,:);
   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
      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);
      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
%  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 '_' day_freq]
   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 '_' day_freq]
   print ('-dpng',figname1) 

%%%  Plot 3
   figure;

   if resolution=='012'
      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;
   hcc=colorbar; set(hcc,'FontS',14);
   caxis ([0.00 0.07])
   title(['\sigma\_ssh (m) - ' namerun_title ' (' dur_plot ')'],'FontSize',18);
   if resolution=='012'
      text(-81,-29,label1, 'HorizontalAlignment','center','FontSize',14) 
   elseif resolution=='025' | resolution=='075'
      text(-81,-30,label1, 'HorizontalAlignment','center','FontSize',14) 
   end
   figname1=['map_var_ssh_' namerun '_' dur_plot '_' day_freq]
   print ('-dpng',figname1) 

quit