clear all; close all;

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

   ZMAX=40; % couche 40 ~ 500m
%  ZMAX=0;  % pour plot complet (debug)

   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_from_obc','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','CARS','U3d','none','U2d','none','V3d','none','V2d','none','T3d','TEMP','S3d','SALT','T2d','none'); 
   data(5)=struct('name','LEVITUS','U3d','none','U2d','none','V3d','none','V2d','none','T3d','temp','S3d','sal','T2d','none'); 

% On plot d2 - d1
   d1=5;
   d2=2;
   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;

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

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

   ifile=0;
%  for word = {'T_tra' 'U_u3d' 'V_u3d'}       
   for word = {'T_tra'}       
      bdy=eval(['bdy_' word{1}]);
      file_ext=bdy.file_ext;

      nfieldmax=size(bdy.field);
      for ii=1:nfieldmax(2)
         ifile=ifile+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;

         if  strcmp(name1,'CARS') | strcmp(name1,'LEVITUS')
            file1=[ root_dir name1 '/bdy' file_ext '_PERU12_' name1 '_mean.nc'] 
         else
            file1=[ root_dir name1 '/bdy' file_ext '_PERU12_' name1 '_y' period '_mean.nc'] 
         end

         if  strcmp(name2,'CARS') | strcmp(name2,'LEVITUS')
            file2=[ root_dir name2 '/bdy' file_ext '_PERU12_' name2 '_mean.nc'] 
         else
            file2=[ root_dir name2 '/bdy' file_ext '_PERU12_' name2 '_y' period '_mean.nc'] 
         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);

         figure(ifile);
         if ZMAX == 0
            pcolor(fieldB2-fieldB1); shading flat;
            set (gca,'YDir','reverse')
            YtextA=90;
            YlineA=75;
         else
            pcolor(1:xbT,zl(1:ZMAX),fieldB2(1:ZMAX,:)-fieldB1(1:ZMAX,:)); shading flat;
            YtextA=-600;
            YlineA=-500;
         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)
         line([342 342],[0 YlineA],'color','k')
         line([844 844],[0 YlineA],'color','k')
         text(170,YtextA,'south','HorizontalAlignment','center','FontSize',14)
         text(600,YtextA,'west','HorizontalAlignment','center','FontSize',14) 
         text(920,YtextA,'north','HorizontalAlignment','center','FontSize',14)
         title([labelA ' ' name2_title ' - ' name1_title ' (' period ')'],'FontSize',18);
         if ZMAX == 0
            figname1=['bdy_' name2 '-' name1 '_' period '_' out_name];
         else
            figname1=['bdy_' name2 '-' name1 '_' period '_' out_name '_0-500m'];
         end
         print ('-dpng',figname1)

      end
   end