clear all; close all;

   period='2000-2008';
   root_dir='/loceanfs/pulsation/cholod/peru12_nemo/inputs/bdy/';

   opengl neverselect

   ifile=1;

%    ZMAX=40; % couche 40 ~ 500m
%    ZMAX=75;  % pour plot complet (debug)
 for ZMAX=40:35:75  
%for ZMAX=75:75  
% data
   data(1)=struct('name','soda2','U3d','U','U2d','BTU','V3d','V','V2d','BTV','T3d','TEMP','S3d','SALT','T2d','SSH');
   data(2)=struct('name','tr12_quik','U3d','uo','U2d','uobtrope','V3d','vo','V2d','vobtrope','T3d','thetao','S3d','so','T2d','zos'); 
   data(3)=struct('name','ORCA025_LIM-T323','U3d','vozocrtx','U2d','vobtcrtx','V3d','vomecrty','V2d','vobtcrty','T3d','votemper','S3d','vosaline','T2d','sossheig'); 
   data(4)=struct('name','LEVITUS','U3d','none','U2d','none','V3d','none','V2d','none','T3d','temp','S3d','sal','T2d','none'); 
   data(5)=struct('name','CARS','U3d','none','U2d','none','V3d','none','V2d','none','T3d','TEMP','S3d','SALT','T2d','none'); 
%  data(6)=struct('name','soda2.old','U3d','U','U2d','BTU','V3d','V','V2d','BTV','T3d','TEMP','S3d','SALT','T2d','SSH');
%  data(7)=struct('name','tr12_quik_from_obc','U3d','uo','U2d','uobtrope','V3d','vo','V2d','vobtrope','T3d','thetao','S3d','so','T2d','zos');

   plot(1)=struct('section','West','is','xb1','ie','xb2','current','U3B','minmax','0.5','incr','0.1');
   plot(2)=struct('section','South','is','1','ie','xb1','current','V3B','minmax','0.2','incr','0.02');
   plot(3)=struct('section','North','is','xb2','ie','xbT-2','current','V3B','minmax','0.2','incr','0.02');

   bdy_U_u3d=struct('file_ext','U_u3d');
   bdy_U_u3d.field(1)=struct('name','U3d','out_name','U_u3d','labelA','U (m/s)','colormapA',[cold;flipud(hot)],'mini',-0.6,'delta',0.01,'maxi',0.6);
   bdy_V_u3d=struct('file_ext','V_u3d');
   bdy_V_u3d.field(1)=struct('name','V3d','out_name','V_u3d','labelA','V (m/s)','colormapA',[cold;flipud(hot)],'mini',-0.15,'delta',0.01,'maxi',0.15);
   bdy_T_tra=struct('file_ext','T_tra');
   bdy_T_tra.field(1)=struct('name','T3d','out_name','T_tra','labelA','T (Celcius)','colormapA',[cold;flipud(hot)],'mini',-3,'delta',0.1,'maxi',3);
   bdy_T_tra.field(2)=struct('name','S3d','out_name','S_tra','labelA','Salt (PSU)','colormapA',[cold;flipud(hot)],'mini',-0.8,'delta',0.1,'maxi',0.8);

   fcoord='/loceanfs/pulsation/cholod/peru12_nemo/inputs/bdy/bdy_coordinates_PERU12.nc';
   nccoord=netcdf(fcoord,'r');
   nav_lon=nccoord{'glamt'}(1,:);
   nav_lat=nccoord{'gphit'}(1,:);
   close(nccoord);

%  list = {'T_u2d' 'U_u2d' 'V_u2d' 'T_tra' 'U_u3d' 'V_u3d'}       
   list = {'T_tra' 'U_u3d' 'V_u3d'}
%  list = {'U_u3d' 'V_u3d'}
%  list = {'T_tra'}

   ndatamax=size(data);
%  for d2=1:ndatamax(2)
   for d2=5:5
% On plot d2 - d1
%  d2=1;
   d1=4;   %   TOUJOURS METTRE LES OBS en d1!!!
   if (d2 == 4 | d2 == 5) & (d1 ~= 4 & d1 ~= 5)
      disp('pour diff (run,OBS) mettre OBS dans d1 et non dans d2  => quit'); quit;
   end
      name1=data(d1).name;
      U3d1=data(d1).U3d;
      V3d1=data(d1).V3d;
      T3d1=data(d1).T3d;
      S3d1=data(d1).S3d;

      name2=data(d2).name;
      U3d2=data(d2).U3d;
      V3d2=data(d2).V3d;
      T3d2=data(d2).T3d;
      S3d2=data(d2).S3d;

      name1_title=regexprep(name1,'_','\\_');
      name2_title=regexprep(name2,'_','\\_');

      for word = list
         if (word{1} == 'U_u3d' | word{1} == 'V_u3d') & (strcmp(name1,'CARS') | strcmp(name1,'LEVITUS') | strcmp(name2,'CARS') | strcmp(name2,'LEVITUS'))
            disp('Pas de U et V pour CARS et LEVITUS => saute'); 
         else
            bdy=eval(['bdy_' word{1}]);
            file_ext=bdy.file_ext;
   
            nfieldmax=size(bdy.field);
            for ii=1:nfieldmax(2)
   %        for ii=1:1
               field_name1=eval([bdy.field(ii).name '1']);
               field_name2=eval([bdy.field(ii).name '2']);
               out_name=bdy.field(ii).out_name;
               labelA=bdy.field(ii).labelA;
               colormapA=bdy.field(ii).colormapA;
               mini=bdy.field(ii).mini;
               delta=bdy.field(ii).delta;
               maxi=bdy.field(ii).maxi;
      
               clear fieldA fieldB fieldBmean fieldC;
               if     word{1} == 'T_tra' | word{1} == 'U_u3d' | word{1} == 'V_u3d' 
                  if  strcmp(name1,'CARS') | strcmp(name1,'LEVITUS')
                     if     word{1} == 'T_tra' 
                        file1=[ root_dir name1 '/bdy' file_ext '_PERU12_' name1 '_mean.nc'] 
                     else
                        disp('Pas de U et V pour CARS et LEVITUS => quit'); quit; 
                     end
                  else
                     file1=[ root_dir name1 '/bdy' file_ext '_PERU12_' name1 '_y' period '_mean.nc'] 
                     if     word{1} == 'T_tra' 
                        fileU1=[ root_dir name1 '/bdy' 'U_u3d' '_PERU12_' name1 '_y' period '_mean.nc'] 
                        fileV1=[ root_dir name1 '/bdy' 'V_u3d' '_PERU12_' name1 '_y' period '_mean.nc'] 
                        ncu1=netcdf(fileU1,'r');
                        U3A1=ncu1{U3d1}(:,:,1,:);
                        U3B1=squeeze(U3A1(:,:));
                        close(ncu1);
                        ncv1=netcdf(fileV1,'r');
                        V3A1=ncv1{V3d1}(:,:,1,:);
                        V3B1=squeeze(V3A1(:,:));
                        close(ncv1);
                     end
                  end
   
                  if  strcmp(name2,'CARS') | strcmp(name2,'LEVITUS')
                     if     word{1} == 'T_tra' 
                        file2=[ root_dir name2 '/bdy' file_ext '_PERU12_' name2 '_mean.nc'] 
                     else
                        disp('Pas de U et V pour CARS et LEVITUS => quit'); quit; 
                     end
                  else
                     file2=[ root_dir name2 '/bdy' file_ext '_PERU12_' name2 '_y' period '_mean.nc'] 
                     if     word{1} == 'T_tra' 
                        fileU2=[ root_dir name2 '/bdy' 'U_u3d' '_PERU12_' name2 '_y' period '_mean.nc'] 
                        fileV2=[ root_dir name2 '/bdy' 'V_u3d' '_PERU12_' name2 '_y' period '_mean.nc'] 
                        ncu2=netcdf(fileU2,'r');
                        U3A2=ncu2{U3d2}(:,:,1,:);
                        U3B2=squeeze(U3A2(:,:));
                        close(ncu2);
                        ncv2=netcdf(fileV2,'r');
                        V3A2=ncv2{V3d2}(:,:,1,:);
                        V3B2=squeeze(V3A2(:,:));
                        close(ncv2);
                     end
                  end
   
                  nc1=netcdf(file1,'r');
                  fieldA1=nc1{field_name1}(:,:,1,:);
                  fieldB1=squeeze(fieldA1(:,:));
                  zl=-1*nc1{'deptht'}(:);
                  xbTB1=size(fieldB1);
                  xbT=xbTB1(2);
                  close(nc1);
   
                  nc2=netcdf(file2,'r');
                  fieldA2=nc2{field_name2}(:,:,1,:);
                  fieldB2=squeeze(fieldA2(:,:));
                  close(nc2);
   
                  xb1=342;
                  xb2=844;
   
                  for num=1:3
%                 for num=1:1
                     section=plot(num).section;
                     is=eval(plot(num).is);
                     ie=eval(plot(num).ie);
%                    if     word{1} == 'T_tra' 
                     if     word{1} == 'T_tra' & not(strcmp(name1,'CARS') | strcmp(name1,'LEVITUS')) 
                        current1=eval([plot(num).current '1']);
                     end
                     if     word{1} == 'T_tra' & not(strcmp(name2,'CARS') | strcmp(name2,'LEVITUS')) 
                        current2=eval([plot(num).current '2']);
                     end
                     minmax=eval(plot(num).minmax);
                     incr=eval(plot(num).incr);
                     if strcmp(section,'West');
                        xaxis=nav_lat(is:ie);
                     else
                        xaxis=nav_lon(is:ie);
                     end
                     figure(ifile);
                     disp(['ifile = ' num2str(ifile)]);
   
                     if ZMAX == 75
                        pcolor(xaxis,1:ZMAX,fieldB2(1:ZMAX,is:ie)-fieldB1(1:ZMAX,is:ie)); shading flat;
                        set (gca,'YDir','reverse')    
                        if     word{1} == 'T_tra' 
                           if not(strcmp(name1,'CARS')) & not(strcmp(name1,'LEVITUS'))  
                              hold on
                              [ctt1,htt1]=contour(xaxis,1:ZMAX,current2(1:ZMAX,is:ie)-current1(1:ZMAX,is:ie),[-minmax:incr:0]); set(htt1,'Color','k','linewidth',1,'LineStyle','--');
                              [ctt2,htt2]=contour(xaxis,1:ZMAX,current2(1:ZMAX,is:ie)-current1(1:ZMAX,is:ie),[0:incr:minmax]); set(htt2,'Color','k','linewidth',1,'LineStyle','-');
                              [ctt3,htt3]=contour(xaxis,1:ZMAX,current2(1:ZMAX,is:ie)-current1(1:ZMAX,is:ie),[0:0]); set(htt3,'Color','k','linewidth',1,'LineStyle','-');
                              clabel(ctt1,htt1,'FontS',12,'Color','k');
                              clabel(ctt2,htt2,'FontS',12,'Color','k');
                              clabel(ctt3,htt3,'FontS',12,'Color','k');
                              line([0 0],[0 75],'color','k')
                           elseif (strcmp(name1,'CARS') | strcmp(name1,'LEVITUS')) & (strcmp(name2,'CARS') | strcmp(name2,'LEVITUS'))
                              disp('');
                           else
                              hold on
                              [ctt1,htt1]=contour(xaxis,1:ZMAX,current2(1:ZMAX,is:ie),[-minmax:incr:0]); set(htt1,'Color','k','linewidth',1,'LineStyle','--');
                              [ctt2,htt2]=contour(xaxis,1:ZMAX,current2(1:ZMAX,is:ie),[0:incr:minmax]); set(htt2,'Color','k','linewidth',1,'LineStyle','-');
                              [ctt3,htt3]=contour(xaxis,1:ZMAX,current2(1:ZMAX,is:ie),[0:0]); set(htt3,'Color','k','linewidth',1,'LineStyle','-');
                              clabel(ctt1,htt1,'FontS',12,'Color','k');
                              clabel(ctt2,htt2,'FontS',12,'Color','k');
                              clabel(ctt3,htt3,'FontS',12,'Color','k');
                              line([0 0],[0 75],'color','k')
                           end
                        end
   
                     elseif  ZMAX == 40
                        pcolor(xaxis,zl(1:ZMAX),fieldB2(1:ZMAX,is:ie)-fieldB1(1:ZMAX,is:ie)); shading flat;
%                       if strcmp(section,'South'); set (gca,'XDir','reverse'); end
                        if     word{1} == 'T_tra' 
                           if not(strcmp(name1,'CARS')) & not(strcmp(name1,'LEVITUS')) 
                              hold on
      		              [ctt1,htt1]=contour(xaxis,zl(1:ZMAX),current2(1:ZMAX,is:ie)-current1(1:ZMAX,is:ie),[-minmax:incr:0]); set(htt1,'Color','k','linewidth',1,'LineStyle','--');
                              [ctt2,htt2]=contour(xaxis,zl(1:ZMAX),current2(1:ZMAX,is:ie)-current1(1:ZMAX,is:ie),[0:incr:minmax]); set(htt2,'Color','k','linewidth',1,'LineStyle','-');
                              [ctt3,htt3]=contour(xaxis,zl(1:ZMAX),current2(1:ZMAX,is:ie)-current1(1:ZMAX,is:ie),[0:0]); set(htt3,'Color','k','linewidth',1,'LineStyle','-');
                              clabel(ctt1,htt1,'FontS',12,'Color','k');
                              clabel(ctt2,htt2,'FontS',12,'Color','k');
                              clabel(ctt3,htt3,'FontS',12,'Color','k');
                           elseif (strcmp(name1,'CARS') | strcmp(name1,'LEVITUS')) & (strcmp(name2,'CARS') | strcmp(name2,'LEVITUS'))
                              disp('');
                           else
                              hold on
                              [ctt1,htt1]=contour(xaxis,zl(1:ZMAX),current2(1:ZMAX,is:ie),[-minmax:incr:0]); set(htt1,'Color','k','linewidth',1,'LineStyle','--');
                              [ctt2,htt2]=contour(xaxis,zl(1:ZMAX),current2(1:ZMAX,is:ie),[0:incr:minmax]); set(htt2,'Color','k','linewidth',1,'LineStyle','-');
                              [ctt3,htt3]=contour(xaxis,zl(1:ZMAX),current2(1:ZMAX,is:ie),[0:0]); set(htt3,'Color','k','linewidth',1,'LineStyle','-');
                              clabel(ctt1,htt1,'FontS',12,'Color','k');
                              clabel(ctt2,htt2,'FontS',12,'Color','k');
                              clabel(ctt3,htt3,'FontS',12,'Color','k');
                              line([0 0],[0 75],'color','k')
                           end
                        end
                     end
   
                     set(gca,'FontS',14); set(gca,'PlotBoxAspectRatio',[1 0.4 1]); box on;
                     caxis ([mini maxi]);
                     hcc=colorbar; set(hcc,'FontS',14);
                     colormap(colormapA)
                     title([labelA ' ' name2_title ' - ' name1_title ' (' period ') ' section],'FontSize',12);
                     if ZMAX == 75
                        figname1=['bdy_' name2 '-' name1 '_' period '_' out_name '_' section];
                     elseif ZMAX == 40
                        figname1=['bdy_' name2 '-' name1 '_' period '_' out_name '_' section '_0-500m'];
                     end
                     print ('-dpng',figname1)
                     ifile=ifile+1;
                  end %  for num=1:3
   
               elseif  word{1} == 'T_u2d' | word{1} == 'U_u2d' | word{1} == 'V_u2d' 
                  file=[ root_dir name '/bdy' file_ext '_PERU12_' name '_y' period '.nc'] 
   
                  nc=netcdf(file,'r');
                  fieldA=nc{field_name}(:,1,:);
                  fieldB=squeeze(fieldA(:,:));
                  close(nc);
   
                  figure(ifile);
                  disp(['ifile = ' num2str(ifile)]);
                  pcolor(fieldB); shading flat;
                  caxis ([mini maxi]); hcc=colorbar; set(hcc,'FontS',14); colormap(colormapA)
                  line([342 342],[0 657],'color','k')
                  line([844 844],[0 657],'color','k')
                  text(170,-50,'south','HorizontalAlignment','center','FontSize',14)
                  text(600,-50,'west','HorizontalAlignment','center','FontSize',14) 
                  text(920,-50,'north','HorizontalAlignment','center','FontSize',14)
                  title([labelA ' ' name_title ' (' period ')'],'FontSize',12);
                  figname1=['bdy_' name '_' period '_' out_name ];
                  print ('-dpng',figname1)
                  ifile=ifile+1;
   
                  if  word{1} == 'T_u2d'
                     fieldBmean=mean(fieldB,1);
                     for kk=1:size(fieldB,2) 
                        fieldC(:,kk)=fieldB(:,kk)-fieldBmean(kk);
                     end
                     figure(ifile);
                     disp(['ifile = ' num2str(ifile)]);
                     pcolor(fieldC); shading flat;
                     caxis ([-0.2 0.2]); hcc=colorbar; set(hcc,'FontS',14); colormap(colormapA)
                     line([342 342],[0 657],'color','k')
                     line([844 844],[0 657],'color','k')
                     text(170,-50,'south','HorizontalAlignment','center','FontSize',14)
                     text(600,-50,'west','HorizontalAlignment','center','FontSize',14)
                     text(920,-50,'north','HorizontalAlignment','center','FontSize',14)
                     title([labelA '-mean ' name_title ' (' period ')'],'FontSize',12);
                     figname1=['bdy_' name '_' period '_' out_name '-mean' ];
                     print ('-dpng',figname1)
                     ifile=ifile+1;
                  end  
               end    %if     word{1} == 'T_tra' | word{1} == 'U_u3d' | word{1} == 'V_u3d' 
            end   %for ii=1:nfieldmax(2)
         end    %if (word{1} == 'U_u3d' | word{1} == 'V_u3d') & (strcmp(name,'CARS') | strcmp(name,'LEVITUS'))
      end    %for word = list
   end    %for d2=1:ndatamax(2)
end    %ZMAX