clear all; close all;

   addpath ./utils

   ys=2003;
   ye1=2007;

   multi_year='TRUE';
%  multi_year='FALS';

%  OBS1='';
   OBS1='CARS';
%  OBS2='';
   OBS2='LEVITUS';
%  DIFF='';
   DIFF='TRUE';
   
   temp='thetao';

   exp(1).config='peru12_nemo';
   exp(1).grid='roms12nemo';
   exp(1).names={'p12n_bdy_tr12n_quik' 'p12n_bdy_soda2' 'p12n_bdy_soda2_B' 'p12n_bdy_tr12n_quik_ND' 'p12n_bdy_tr12n_quik_blp1' };
   exp(1).time='time_average_1mo';

   exp(2).config='trop12_nemo/peru';
   exp(2).grid='roms12nemo';
   exp(2).names={'tr12_quik' 'tr12_erainocrt'};
   exp(2).time='time_average_5d';

   exp(3).config='trop025peru12_nemo/peru';
   exp(3).grid='roms12nemo';
   exp(3).names={'tr025p12_quik'};
   exp(3).time='time_average_5d';

   exp(4).config='peru12_nemo';
   exp(4).grid='roms12nemo';
%  exp(4).names={'erai70b' 'erai74c'};
%  exp(4).names={'erai70a' 'erai70b' 'erai70c' 'erai72c' 'erai73c' 'erai73c2' 'erai74c' 'erai75c'};
   exp(4).names={'erai70c' 'erai72c' 'erai73c' 'erai73c2' 'erai75c'};
   exp(4).time='time_average_5d';

   exp_dim=size(exp);
   expREF=2;
   numREF=1;

   if strcmp(OBS1,'CARS')
      MfiletCARS='/loceanfs/pulsation/cholod/data/CARS/cars_mon_clim_MEAN_PERU12_time5.nc';
      ctlMtCARS=rnt_timectl({MfiletCARS}, 'time_counter');
   end
   if strcmp(OBS2,'LEVITUS')
      MfiletLEVITUS='/loceanfs/pulsation/cholod/data/LEVITUS/levitus_annual_mean_PERU12_time5.nc';
      ctlMtLEVITUS=rnt_timectl({MfiletLEVITUS}, 'time_counter');
   end

   if strcmp(multi_year,'TRUE')   % ruse pour réaliser soit une boucle de ys à ye1 soit un seul passage pour ys-ye1
      ye=ys;
   else
      ye=ye1;
   end

   for year=ys:1:ye  
      if strcmp(multi_year,'TRUE')
         period=[num2str(ys) '-' num2str(ye1)];
         rep=[num2str(ys) '0101_' num2str(ye1) '1231'];
      else
         period=num2str(year);
         rep=[period '0101_' period '1231'];
      end

%     for ii = 1:exp_dim(2)
      for ii = 1:1
         config=exp(ii).config;
         grid=exp(ii).grid;
         names=exp(ii).names;
         time=exp(ii).time;
   
         names_dim=size(names);
         grd=rnt_gridload_xa(grid); 
   
%        for jj = 1:names_dim(2);
         for jj = 1:1
            name=names{jj};
            name_title=regexprep(name,'_','\\_');

            dirfile=['/loceanfs/pulsation/cholod/' config '/outputs/' name '/' rep '/']; 
            Mfileu=[dirfile name '_1s_' rep '_grid_U_3D.nc'] 
            Mfilev=[dirfile name '_1s_' rep '_grid_V_3D.nc'] 
            Mfilet=[dirfile name '_1s_' rep '_grid_T_T3D.nc']  
            ctlMu=rnt_timectl({Mfileu}, time);
            ctlMv=rnt_timectl({Mfilev}, time);
            ctlMt=rnt_timectl({Mfilet}, time);

            if strcmp(DIFF,'TRUE')
               configREF=exp(expREF).config;
               gridREF=exp(expREF).grid;
               nameREF=exp(expREF).names{numREF};
               timeREF=exp(expREF).time;
               dirfileREF=['/loceanfs/pulsation/cholod/' configREF '/outputs/' nameREF '/' rep '/']; 
               MfileuREF=[dirfileREF nameREF '_1s_' rep '_grid_U_3D.nc'] 
               MfilevREF=[dirfileREF nameREF '_1s_' rep '_grid_V_3D.nc'] 
               MfiletREF=[dirfileREF nameREF '_1s_' rep '_grid_T_T3D.nc']  
               ctlMuREF=rnt_timectl({MfileuREF}, timeREF);
               ctlMvREF=rnt_timectl({MfilevREF}, timeREF);
               ctlMtREF=rnt_timectl({MfiletREF}, timeREF);
               name_titleREF=regexprep([name ' - ' nameREF],'_','\\_');
               fignameREF=['peru_' name '-' nameREF '_' period '_'];
            end

            season=[5]; 
            figname=['peru_' name '_' period '_'];
            vname='r2v';  prename='along_'; dist=3e5;
            latpm=3;
            latsec=-10;
            %%latsec=[-10 -25 -31 ];
   
            nc=netcdf(Mfilet,'r');
            zl=-1*nc{'deptht'}(:);
            close(nc);
            NN=length(zl);

            for k=1:length(latsec)
               latmin=latsec(k)-latpm; latmax=latsec(k)+latpm;
               [Mvar,Mmask,Mh,distx]=get_section_meanlongshore_nemo(grd,ctlMu,ctlMv,ctlMt,vname,season,dist,latmin,latmax,NN);
               [MTemp,Mmask,Mh,distx]=get_section_meanlongshore_nemo(grd,ctlMu,ctlMv,ctlMt,temp,season,dist,latmin,latmax,NN);
               if strcmp(OBS1,'CARS')
                  [MTempCARS,Mmask,Mh,distx]=get_section_meanlongshore_nemo(grd,ctlMu,ctlMv,ctlMtCARS,temp,season,dist,latmin,latmax,NN);
               end
               if strcmp(OBS2,'LEVITUS')
                  [MTempLEVITUS,Mmask,Mh,distx]=get_section_meanlongshore_nemo(grd,ctlMu,ctlMv,ctlMtLEVITUS,temp,season,dist,latmin,latmax,NN);
               end
               if strcmp(DIFF,'TRUE')
                  [MvarREF,Mmask,Mh,distx]=get_section_meanlongshore_nemo(grd,ctlMuREF,ctlMvREF,ctlMtREF,vname,season,dist,latmin,latmax,NN);
                  [MTempREF,Mmask,Mh,distx]=get_section_meanlongshore_nemo(grd,ctlMuREF,ctlMvREF,ctlMtREF,temp,season,dist,latmin,latmax,NN);
               end
      

      
% plot currents profondeur 0-300m / CARS
               Tmin=8; Tmax=22; Tinc=1;
               figure(1);
               ninterv=30;
               contourf(distx/1000,zl,Mvar,ninterv); shading flat
               hold on
               caxis([-0.1 0.1]); 
               axis([-250 0 -300 0]);
               plot(distx/1000,-Mh,'k');
%              colormap(jet(ninterv));
               colormap([cold;flipud(hot)]);
               hy=ylabel('Depth [m]'); set(hy,'FontS',12); hx=xlabel('Distance [km]'); set(hx,'FontS',12);
               set(gca,'FontS',12); set(gca,'PlotBoxAspectRatio',[0.7 1 1]); wysiwyg; box on;
               [ctt,htt]=contour(distx/1000,zl,MTemp,[Tmin:Tinc:Tmax]); set(htt,'LineColor',[0.0 0.0 0.0]);
               clabel(ctt,htt,'FontS',12,'LabelSpacing',250,'Color',[0.0 0.0 0.0]);
               hcc=colorbar; set(hcc,'FontS',12);
               title(['Current(m/s) / T(C)- ' name_title ' (' period ')'], 'FontSize',12);
            
               if strcmp(OBS1,'CARS')
                  [cttCARS,httCARS]=contour(distx/1000,zl,MTempCARS,[Tmin:Tinc:Tmax]); set(httCARS,'LineColor',[0.0 0.0 1.0]);
                  clabel(cttCARS,httCARS,'FontS',12,'LabelSpacing',290,'Color',[0.0 0.0 1.0]);
                  figname2=[prename figname 'seas' num2str(season) '_lat' num2str(latsec(k)) '_currents_0-300m+CARS'];
               else
                  figname2=[prename figname 'seas' num2str(season) '_lat' num2str(latsec(k)) '_currents_0-300m'];
               end
               print ('-dpng',figname2);

% plot currents profondeur 0-300m / REF
            if strcmp(DIFF,'TRUE')
               Tmin=8; Tmax=22; Tinc=1;
               figure(10);
               ninterv=30;
               contourf(distx/1000,zl,Mvar-MvarREF,ninterv); shading flat
               hold on
               caxis([-0.1 0.1]); 
               axis([-250 0 -300 0]);
               plot(distx/1000,-Mh,'k');
%              colormap(redblue(ninterv));
               colormap([cold;flipud(hot)]);
               hy=ylabel('Depth [m]'); set(hy,'FontS',12); hx=xlabel('Distance [km]'); set(hx,'FontS',12);
               set(gca,'FontS',12); set(gca,'PlotBoxAspectRatio',[0.7 1 1]); wysiwyg; box on;
               [ctt,htt]=contour(distx/1000,zl,MTemp,[Tmin:Tinc:Tmax]); set(htt,'LineColor',[0.0 0.0 0.0]);
               clabel(ctt,htt,'FontS',12,'LabelSpacing',250,'Color',[0.0 0.0 0.0]);
               hcc=colorbar; set(hcc,'FontS',12);
               title(['Current(m/s) / T(C)- ' name_titleREF ' (' period ')'], 'FontSize',12);
            
               [cttREF,httREF]=contour(distx/1000,zl,MTempREF,[Tmin:Tinc:Tmax]); set(httREF,'LineColor',[0.0 0.0 1.0]);
               clabel(cttREF,httREF,'FontS',12,'LabelSpacing',290,'Color',[0.0 0.0 1.0]);
               figname2=[prename fignameREF 'seas' num2str(season) '_lat' num2str(latsec(k)) '_currents_0-300m'];
               print ('-dpng',figname2);
            end


% plot currents 75 niveaux / LEVITUS
               layers=[1:75];
               Tmin=0; Tmax=22; Tinc=1;
               figure(2);
               ninterv=30;
               contourf(distx/1000,layers,Mvar,ninterv); shading flat
               set (gca,'YDir','reverse')
               hold on
               caxis([-0.1 0.1]); 
               plot(distx/1000,-Mh,'k');
%              colormap(jet(ninterv));
               colormap([cold;flipud(hot)]);
               hy=ylabel('Layer'); set(hy,'FontS',12); hx=xlabel('Distance [km]'); set(hx,'FontS',12);
               set(gca,'FontS',12); set(gca,'PlotBoxAspectRatio',[0.7 1 1]); wysiwyg; box on;
               [ctt,htt]=contour(distx/1000,layers,MTemp,[Tmin:Tinc:Tmax]); set(htt,'LineColor',[0.0 0.0 0.0]);
               clabel(ctt,htt,'FontS',12,'LabelSpacing',250,'Color',[0.0 0.0 0.0]);
               hcc=colorbar; set(hcc,'FontS',12);
               title(['Current(m/s) / T(C)- ' name_title ' (' period ')'], 'FontSize',12);
            
               if strcmp(OBS2,'LEVITUS')
                  [cttLEVITUS,httLEVITUS]=contour(distx/1000,layers,MTempLEVITUS,[Tmin:Tinc:Tmax]); set(httLEVITUS,'LineColor',[0.0 0.0 1.0]);
                  clabel(cttLEVITUS,httLEVITUS,'FontS',12,'LabelSpacing',290,'Color',[0.0 0.0 1.0]);
                  figname2=[prename figname 'seas' num2str(season) '_lat' num2str(latsec(k)) '_currents_L75+LEVITUS'];
               else
                  figname2=[prename figname 'seas' num2str(season) '_lat' num2str(latsec(k)) '_currents_L75'];
               end
               print ('-dpng',figname2);

% plot currents 75 niveaux / REF
            if strcmp(DIFF,'TRUE')
               layers=[1:75];
               Tmin=0; Tmax=22; Tinc=1;
               figure(20);
               ninterv=30;
               contourf(distx/1000,layers,Mvar-MvarREF,ninterv); shading flat
               set (gca,'YDir','reverse')
               hold on
               caxis([-0.1 0.1]); 
               plot(distx/1000,-Mh,'k');
%              colormap(jet(ninterv));
               colormap([cold;flipud(hot)]);
               hy=ylabel('Layer'); set(hy,'FontS',12); hx=xlabel('Distance [km]'); set(hx,'FontS',12);
               set(gca,'FontS',12); set(gca,'PlotBoxAspectRatio',[0.7 1 1]); wysiwyg; box on;
               [ctt,htt]=contour(distx/1000,layers,MTemp,[Tmin:Tinc:Tmax]); set(htt,'LineColor',[0.0 0.0 0.0]);
               clabel(ctt,htt,'FontS',12,'LabelSpacing',250,'Color',[0.0 0.0 0.0]);
               hcc=colorbar; set(hcc,'FontS',12);
               title(['Current(m/s) / T(C)- ' name_titleREF ' (' period ')'], 'FontSize',12);
            
               [cttREF,httREF]=contour(distx/1000,layers,MTempREF,[Tmin:Tinc:Tmax]); set(httREF,'LineColor',[0.0 0.0 1.0]);
               clabel(cttREF,httREF,'FontS',12,'LabelSpacing',290,'Color',[0.0 0.0 1.0]);
               figname2=[prename fignameREF 'seas' num2str(season) '_lat' num2str(latsec(k)) '_currents_L75'];
               print ('-dpng',figname2);
            end

% plot diff temperature profondeur 0-300m / CARS
               figure(3);
               minmax=3;
               division=0.25;
               [ctt,htt]=contourf(distx/1000,zl,MTemp-MTempCARS,[-minmax:division:minmax]);
               hold on
               contour(distx/1000,zl,MTemp-MTempCARS,[-minmax:division:minmax]);
               [ctt2,htt2]=contour(distx/1000,zl,MTemp,[Tmin:Tinc:Tmax]); set(htt2,'LineColor',[0.0 0.0 0.0]);
               clabel(ctt2,htt2,'FontS',14,'LabelSpacing',250,'Color',[0.0 0.0 0.0]);
               set(htt,'LineColor',[0.0 0.0 0.0]); 
               axis([-250 0 -300 0]);
               caxis([-minmax minmax]); 
               ncolor=minmax*2/division;
               cc=colormap(redblue(ncolor+2));
               borne=minmax/division;
               cc1=cc(1:borne,:);
               cc2=cc(borne+3:end,:);
               ccc=vertcat(cc1,cc2);
               colormap(ccc);
               plot(distx/1000,-Mh,'k');
               hy=ylabel('Depth [m]'); set(hy,'FontS',14); hx=xlabel('Distance [km]'); set(hx,'FontS',14);
               set(gca,'FontS',14); set(gca,'PlotBoxAspectRatio',[0.7 1 1]); wysiwyg; box on;
               hcc=colorbar; set(hcc,'FontS',14);
               title(['T (C)- ' name_title '-CARS (' period ')'], 'FontSize',12);
            
               if strcmp(OBS1,'CARS')
                  [cttCARS,httCARS]=contour(distx/1000,zl,MTempCARS,[Tmin:Tinc:Tmax]); set(httCARS,'LineColor',[0.0 0.0 1.0]);
                  clabel(cttCARS,httCARS,'FontS',14,'LabelSpacing',290,'Color',[0.0 0.0 1.0]);
                  figname2=[prename figname 'seas' num2str(season) '_lat' num2str(latsec(k)) '_temperature_0-300m+CARS'];
               else
                  figname2=[prename figname 'seas' num2str(season) '_lat' num2str(latsec(k)) '_temperature_0-300m'];
               end
               print ('-dpng',figname2);
         

% plot diff temperature profondeur 0-300m / REF
            if strcmp(DIFF,'TRUE')
               figure(30);
               minmax=3;
               division=0.25;
               [ctt,htt]=contourf(distx/1000,zl,MTemp-MTempREF,[-minmax:division:minmax]);
               hold on
               contour(distx/1000,zl,MTemp-MTempREF,[-minmax:division:minmax]);
               [ctt2,htt2]=contour(distx/1000,zl,MTemp,[Tmin:Tinc:Tmax]); set(htt2,'LineColor',[0.0 0.0 0.0]);
               clabel(ctt2,htt2,'FontS',14,'LabelSpacing',250,'Color',[0.0 0.0 0.0]);
               set(htt,'LineColor',[0.0 0.0 0.0]); 
               axis([-250 0 -300 0]);
               caxis([-minmax minmax]); 
               ncolor=minmax*2/division;
               cc=colormap(redblue(ncolor+2));
               borne=minmax/division;
               cc1=cc(1:borne,:);
               cc2=cc(borne+3:end,:);
               ccc=vertcat(cc1,cc2);
               colormap(ccc);
               plot(distx/1000,-Mh,'k');
               hy=ylabel('Depth [m]'); set(hy,'FontS',14); hx=xlabel('Distance [km]'); set(hx,'FontS',14);
               set(gca,'FontS',14); set(gca,'PlotBoxAspectRatio',[0.7 1 1]); wysiwyg; box on;
               hcc=colorbar; set(hcc,'FontS',14);
               title(['T (C)- ' name_titleREF ' (' period ')'], 'FontSize',12);
            
               [cttREF,httREF]=contour(distx/1000,zl,MTempREF,[Tmin:Tinc:Tmax]); set(httREF,'LineColor',[0.0 0.0 1.0]);
               clabel(cttREF,httREF,'FontS',14,'LabelSpacing',290,'Color',[0.0 0.0 1.0]);
               figname2=[prename fignameREF 'seas' num2str(season) '_lat' num2str(latsec(k)) '_temperature_0-300m'];
               print ('-dpng',figname2);
            end
         
% plot diff temperature 75 niveaux / LEVITUS
               layers=[1:75];
               figure(4);
               minmax=3;
               division=0.25;
               [ctt,htt]=contourf(distx/1000,layers,MTemp-MTempLEVITUS,[-minmax:division:minmax]);
               set (gca,'YDir','reverse')
               hold on
               contour(distx/1000,layers,MTemp-MTempLEVITUS,[-minmax:division:minmax]);
               [ctt2,htt2]=contour(distx/1000,layers,MTemp,[Tmin:Tinc:Tmax]); set(htt2,'LineColor',[0.0 0.0 0.0]);
               clabel(ctt2,htt2,'FontS',14,'LabelSpacing',250,'Color',[0.0 0.0 0.0]);
               set(htt,'LineColor',[0.0 0.0 0.0]); 
               caxis([-minmax minmax]); 
               ncolor=minmax*2/division;
               cc=colormap(redblue(ncolor+2));
               borne=minmax/division;
               cc1=cc(1:borne,:);
               cc2=cc(borne+3:end,:);
               ccc=vertcat(cc1,cc2);
               colormap(ccc);
               plot(distx/1000,-Mh,'k');
               hy=ylabel('Layers'); set(hy,'FontS',14); hx=xlabel('Distance [km]'); set(hx,'FontS',14);
               set(gca,'FontS',14); set(gca,'PlotBoxAspectRatio',[0.7 1 1]); wysiwyg; box on;
               hcc=colorbar; set(hcc,'FontS',14);
               title(['T (C)- ' name_title '-LEVITUS (' period ')'], 'FontSize',12);

               if strcmp(OBS2,'LEVITUS')
                  [cttLEVITUS,httLEVITUS]=contour(distx/1000,layers,MTempLEVITUS,[Tmin:Tinc:Tmax]); set(httLEVITUS,'LineColor',[0.0 0.0 1.0]);
                  clabel(cttLEVITUS,httLEVITUS,'FontS',14,'LabelSpacing',290,'Color',[0.0 0.0 1.0]);
                  figname2=[prename figname 'seas' num2str(season) '_lat' num2str(latsec(k)) '_temperature_L75+LEVITUS'];
               else
                  figname2=[prename figname 'seas' num2str(season) '_lat' num2str(latsec(k)) '_temperature_L75'];
               end
               print ('-dpng',figname2);
         
% plot diff temperature 75 niveaux / REF
            if strcmp(DIFF,'TRUE')
               layers=[1:75];
               figure(40);
               minmax=3;
               division=0.25;
               [ctt,htt]=contourf(distx/1000,layers,MTemp-MTempREF,[-minmax:division:minmax]);
               set (gca,'YDir','reverse')
               hold on
               contour(distx/1000,layers,MTemp-MTempREF,[-minmax:division:minmax]);
               [ctt2,htt2]=contour(distx/1000,layers,MTemp,[Tmin:Tinc:Tmax]); set(htt2,'LineColor',[0.0 0.0 0.0]);
               clabel(ctt2,htt2,'FontS',14,'LabelSpacing',250,'Color',[0.0 0.0 0.0]);
               set(htt,'LineColor',[0.0 0.0 0.0]); 
               caxis([-minmax minmax]); 
               ncolor=minmax*2/division;
               cc=colormap(redblue(ncolor+2));
               borne=minmax/division;
               cc1=cc(1:borne,:);
               cc2=cc(borne+3:end,:);
               ccc=vertcat(cc1,cc2);
               colormap(ccc);
               plot(distx/1000,-Mh,'k');
               hy=ylabel('Layers'); set(hy,'FontS',14); hx=xlabel('Distance [km]'); set(hx,'FontS',14);
               set(gca,'FontS',14); set(gca,'PlotBoxAspectRatio',[0.7 1 1]); wysiwyg; box on;
               hcc=colorbar; set(hcc,'FontS',14);
               title(['T (C)- ' name_titleREF ' (' period ')'], 'FontSize',12);

               [cttREF,httREF]=contour(distx/1000,layers,MTempREF,[Tmin:Tinc:Tmax]); set(httREF,'LineColor',[0.0 0.0 1.0]);
               clabel(cttREF,httREF,'FontS',14,'LabelSpacing',290,'Color',[0.0 0.0 1.0]);
               figname2=[prename fignameREF 'seas' num2str(season) '_lat' num2str(latsec(k)) '_temperature_L75'];
               print ('-dpng',figname2);
            end

         
            if strcmp(DIFF,'TRUE')
               close(40);
               close(30);
               close(20);
               close(10);
            end
               close(4);
               close(3);
               close(2);
               close(1);
            end  % end of latsec loop
         end   % jj
      end   % ii
   end     %   for year