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','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'); 
%  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');

% Plots
   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_U_u2d=struct('file_ext','U_u2d');
   bdy_U_u2d.field(1)=struct('name','U2d','out_name','U_u2d','labelA','U (m/s)','colormapA',[cold;flipud(hot)],'mini',-0.06,'delta',0.01,'maxi',0.06);

   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_V_u2d=struct('file_ext','V_u2d');
   bdy_V_u2d.field(1)=struct('name','V2d','out_name','V_u2d','labelA','V (m/s)','colormapA',[cold;flipud(hot)],'mini',-0.06,'delta',0.01,'maxi',0.06);

   bdy_T_tra=struct('file_ext','T_tra');
   bdy_T_tra.field(1)=struct('name','T3d','out_name','T_tra','labelA','T (Celcius)','colormapA','default','mini',0,'delta',1,'maxi',30);
   bdy_T_tra.field(2)=struct('name','S3d','out_name','S_tra','labelA','Salt (PSU)','colormapA','default','mini',33,'delta',0.2,'maxi',36);
   bdy_T_u2d=struct('file_ext','T_u2d');
   bdy_T_u2d.field(1)=struct('name','T2d','out_name','T_u2d','labelA','SSH (m)','colormapA',[cold;flipud(hot)],'mini',-0.4,'delta',0.01,'maxi',0.4);

   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 = {'T_tra'};

   ndatamax=size(data);
%  for jj=1:ndatamax(2)
   for jj=1:1
      name=data(jj).name;
      U3d=data(jj).U3d;
      U2d=data(jj).U2d;
      V3d=data(jj).V3d;
      V2d=data(jj).V2d;
      T3d=data(jj).T3d;
      S3d=data(jj).S3d;
      T2d=data(jj).T2d;

      name_title=regexprep(name,'_','\\_');
      for word = list
         if (word{1} == 'U_u3d' | word{1} == 'V_u3d') & (strcmp(name,'CARS') | strcmp(name,'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_name=eval(bdy.field(ii).name);
               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(name,'CARS') | strcmp(name,'LEVITUS')
%                    if     word{1} == 'T_tra' 
                        file=[ root_dir name '/bdy' file_ext '_PERU12_' name '_mean.nc'] 
%                    else
%                       disp('Pas de U et V pour CARS et LEVITUS => quit'); quit; 
%                    end
                  else
                     file=[ root_dir name '/bdy' file_ext '_PERU12_' name '_y' period '_mean.nc'] 
                     if     word{1} == 'T_tra' 
                        fileU=[ root_dir name '/bdy' 'U_u3d' '_PERU12_' name '_y' period '_mean.nc'] 
                        fileV=[ root_dir name '/bdy' 'V_u3d' '_PERU12_' name '_y' period '_mean.nc'] 
                        ncu=netcdf(fileU,'r');
                        U3A=ncu{U3d}(:,:,1,:);
                        U3B=squeeze(U3A(:,:));
                        close(ncu);
                        ncv=netcdf(fileV,'r');
                        V3A=ncv{V3d}(:,:,1,:);
                        V3B=squeeze(V3A(:,:));
                        close(ncv);
                     end
                  end
   
                  nc=netcdf(file,'r');
                  fieldA=nc{field_name}(:,:,1,:);
                  fieldB=squeeze(fieldA(:,:));
                  zl=-1*nc{'deptht'}(:);
                  xbTB=size(fieldB);
                  xbT=xbTB(2);
                  xb1=342;
                  xb2=844;
                  close(nc);
   
% plots WEST SOUTH NORTH
                  for num=1:3  
                     section=plot(num).section;
                     is=eval(plot(num).is);
                     ie=eval(plot(num).ie);
                     if     word{1} == 'T_tra' & not(strcmp(name,'CARS') | strcmp(name,'LEVITUS')) 
                        current=eval(plot(num).current);
                     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,fieldB(1:ZMAX,is:ie)); shading flat;
                        set (gca,'YDir','reverse')    
                        line([0 0],[0 75],'color','k')
                        if     word{1} == 'T_tra' & not(strcmp(name,'CARS') | strcmp(name,'LEVITUS')) 
                           hold on
      		        [ctt1,htt1]=contour(xaxis,1:ZMAX,current(1:ZMAX,is:ie),[-minmax:incr:0]); set(htt1,'Color','k','linewidth',1,'LineStyle','--');
                           [ctt2,htt2]=contour(xaxis,1:ZMAX,current(1:ZMAX,is:ie),[0:incr:minmax]); set(htt2,'Color','k','linewidth',1,'LineStyle','-');
                           [ctt3,htt3]=contour(xaxis,1:ZMAX,current(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');
                        end
   
                     elseif  ZMAX == 40
                        pcolor(xaxis,zl(1:ZMAX),fieldB(1:ZMAX,is:ie)); shading flat;
%                       if strcmp(section,'South'); set (gca,'XDir','reverse'); end
                        line([0 0],[0 -500],'color','k')
                        if     word{1} == 'T_tra' & not(strcmp(name,'CARS') | strcmp(name,'LEVITUS')) 
                           hold on
      		        [ctt1,htt1]=contour(xaxis,zl(1:ZMAX),current(1:ZMAX,is:ie),[-minmax:incr:0]); set(htt1,'Color','k','linewidth',1,'LineStyle','--');
                           [ctt2,htt2]=contour(xaxis,zl(1:ZMAX),current(1:ZMAX,is:ie),[0:incr:minmax]); set(htt2,'Color','k','linewidth',1,'LineStyle','-');
                           [ctt3,htt3]=contour(xaxis,zl(1:ZMAX),current(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');
                        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 ' ' name_title ' (' period ') ' section],'FontSize',12);
                     if ZMAX == 75
                        figname1=['bdy_' name '_' period '_' out_name '_' section];
                     elseif ZMAX == 40
                        figname1=['bdy_' name '_' 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 jj=1:ndatamax(2)
end    %ZMAX