namerun='erai51c';
%namerun='erai61c';

kk=0;

% periode du plot temporel
 yy=[1991:1:1994];
%yy=[1991:1:1994];
%yy=1994;

% periode du fichier pluriannuel
 y0=1990;         
%y0=1990;         
 yf=1994;         
%yf=1994;         
ds=((yy(1)-y0)*365)+1 
de=(yy(end)-y0+1)*365 

%rep='<REP>';
%rep=[num2str(y0) '0101_' num2str(yy(end)) '1231'];
 rep=[num2str(y0) '0101_' num2str(yf) '1231'];

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

gridfile=['/net/adonis/usr/adonis/phybiocar/cholod/peru12_nemo/mesh_mask_PERU12_L46.nc'] 
ng=netcdf(gridfile,'r');
tmask=squeeze(ng{'tmaskutil'}(1,:,:));
lat=ng{'nav_lat'}(:);
lon=ng{'nav_lon'}(:);
smt=sum(sum(tmask));  % nombre de pooints d'ocean (1)
close(ng);

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

dirfile=['/net/adonis/usr/adonis/phybiocar/cholod/peru12_nemo/outputs/' namerun '/' rep '/'] 
filezos=[dirfile 'peru12_1d_' rep '_grid_T.' namerun '.nc'] 
nz=netcdf(filezos,'r'); 


display('jour: 1  ');
for ii=ds:1:de
%display(['ii = '  num2str(ii)]);

    kk=kk+1;
    if mod(kk,50)==0; display(['jour: ' num2str(kk) ]); end
%display(['i  ' num2str(i)]);
    zos=squeeze(nz{'zos'}(ii,:,:));

    if ii<=(de-45) 
      mzos=squeeze(mean(nz{'zos'}(ii-45:ii+45,:,:),1));
    else
      mzos=squeeze(mean(nz{'zos'}(ii-45:de,:,:),1));
    end

    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

    vzos=(zos-mzos).^2;
    [gzx,gzy]=gradient(zos-mzos);
    mgz=(gzx.^2+gzy.^2).*tmask;
    vzos=vzos.*tmask;       %  . pour tout calcul matriciel
    rvzos=vzos(50:400,70:end); rtmask=tmask(50:400,70:end); rsmt=sum(sum(rtmask));  
    rmgz=mgz(50:400,70:end);
    ke(kk)=(sum(sum(rvzos))./rsmt);
    sgzos(kk)=(sum(sum(rmgz))./rsmt);
    time(kk)=0.5+(kk-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

  close(nz);


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

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

opengl neverselect

figure;
plot(time,ke);
grid on;
%axis([0 2000 0 0.06]);
xlabel('jours')
year_dur=[num2str(yy(1)) '-' num2str(yy(length(yy)))]
title(['Variance ssh serie (m^2)- ' namerun ' (' year_dur ')'], 'FontSize',15);
%set (gca,'ylim',[0.5e-3 2.5e-3]);
set (gca,'ylim',[0.0e-3 2.0e-3]);
figname1=['var_ssh_' namerun '_' year_dur ]
print ('-dpng',figname1) 

figure;
plot(time,sgzos);
grid on;
%axis([0 2000 0 0.06]);
xlabel('jours')
year_dur=[num2str(yy(1)) '-' num2str(yy(length(yy)))]
title(['Gradient ssh serie (m^2)- ' namerun ' (' year_dur ')'], 'FontSize',15);
%set (gca,'ylim',[0.5e-5 2e-5]);  
set (gca,'ylim',[0.0e-5 6e-5]);  
figname1=['grad_ssh_' namerun '_' year_dur ]
print ('-dpng',figname1) 

figure;
pcolor(lon(50:400,70:end),lat(50:400,70:end),squeeze(mvz(5,50:400,70:end)));
shading flat;   %  enlève les lignes en chaque point
set(gca,'PlotBoxAspectRatio',[0.7 1 1]);  % par défaut 1 et 1 => 0.6 et 1
colorbar;
set(gca,'FontS',14);
%axis([0 2000 0 0.06]);
%caxis ([0 0.8e-2])
caxis ([0 3.0e-3])
year_dur=[num2str(yy(1)) '-' num2str(yy(length(yy)))]
title(['Annual ssh Variance (m^2)- ' namerun ' (' year_dur ')'], 'FontSize',15);
figname1=['map_var_ssh_' namerun '_' year_dur ]
print ('-dpng',figname1) 

 quit