Index: trunk/src/find_field_cells.m
===================================================================
--- trunk/src/find_field_cells.m	(revision 1041)
+++ trunk/src/find_field_cells.m	(revision 1045)
@@ -18,4 +18,5 @@
 %              _vector_x,_y,_z: indices of variables giving the vector components x, y, z
 %              _warnflag: index of warnflag    
+%       .DimIndex
 %      .ProjModeRequest= 'interp_lin', 'interp_tps' indicate whether lin interpolation  or derivatives (tps) is needed to calculate the requested field
 %      .FieldName = operation to be performed to finalise the field cell after projection
@@ -64,5 +65,5 @@
 
 function [CellInfo,NbDim,errormsg]=find_field_cells(Data)
-CellInfo={};
+CellInfo={};%default output
 NbDim=0;
 errormsg='';
@@ -87,5 +88,5 @@
 FieldName=cell(size(Role)); % fieldRequest set to {} by default
 CheckSub=zeros(size(Role));% =1 for fields to substract
-Role=regexprep(Role,' ','scalar'); % Role set to 'scalar' by default
+%Role=regexprep(Role,' ','scalar'); % Role set to 'scalar' by default
 if isfield(Data,'VarAttribute')
     for ivar=1:numel(Data.VarAttribute)
@@ -105,223 +106,590 @@
 end
 
-%% find scattered (unstructured) coordinates
-ivar_coord_x=find(strcmp('coord_x',Role));%find variables with Role='coord_x'
-check_select=false(1,nbvar);
-check_coord=false(1,nbvar);
-CellInfo=cell(1,numel(ivar_coord_x));
-NbDim=zeros(1,numel(ivar_coord_x));
-% loop on unstructured coordinate x -> different field cells
-for icell=1:numel(ivar_coord_x)
-    DimCell=Data.VarDimName{ivar_coord_x(icell)};% cell of dimension names for ivar_coord_x(icell)
-    if ischar(DimCell),DimCell={DimCell};end % transform char to cell for a single dimension
-    % look for variables sharing dimension(s) with ivar_coord_x(icell)
-    check_cell=zeros(numel(DimCell),nbvar);
-    for idim=1:numel(DimCell);% for each variable with role coord_x, look at which other variables contain the same dimension
-        for ivar=1:nbvar
-            check_cell(idim,ivar)=max(strcmp(DimCell{idim},Data.VarDimName{ivar}));
-        end
-    end
-    check_cell=sum(check_cell,1)==numel(DimCell);%logical array=1 for variables belonging to the current cell
-    VarIndex=find(check_cell);% list of detected variable indices
-    if ~(numel(VarIndex)==1 && numel(DimCell)==1)% exclude case of isolated coord_x variable (treated later)
-        if ~(numel(VarIndex)==1 && numel(DimCell)>1)% a variable is associated to coordinate
-            CellInfo{icell}.CoordIndex=ivar_coord_x(icell);
-            % size of coordinate var
+%% detect  fields with different roles
+ind_scalar=find(strcmp('scalar',Role));
+ind_errorflag=find(strcmp('errorflag',Role));
+ind_image=find(strcmp('image',Role));
+ind_vector_x=[find(strcmp('vector_x',Role)) find(strcmp('vector_x_tps',Role))];
+if ~isempty(ind_vector_x)
+    ind_vector_y=[find(strcmp('vector_y',Role)) find(strcmp('vector_y_tps',Role))];
+    ind_vector_z=find(strcmp('vector_z',Role));
+    ind_warnflag=find(strcmp('warnflag',Role));
+    ind_ancillary=find(strcmp('ancillary',Role));
+end
+ind_discrete=find(strcmp('discrete',Role));
+ind_coord_x=find(strcmp('coord_x',Role));
+ind_coord_y=find(strcmp('coord_y',Role));
+ind_coord_z=find(strcmp('coord_z',Role));
+ind_coord_tps=find(strcmp('coord_tps',Role));
+check_string=cellfun(@ischar,Data.VarDimName)==1;
+index_string=find(check_string);
+for ivar=index_string
+    Data.VarDimName{ivar}={Data.VarDimName{ivar}};%transform char strings into cells
+end
+check_coord_names= cellfun(@numel,Data.VarDimName)==1;
+check_coord_raster=false(size(check_coord_names));% check variables describing regular mesh (raster coordinates), from two values, min and max.
+if check_var
+    for ivar=find(check_coord_names)
+        if numel(Data.(Data.ListVarName{ivar}))==2
+            check_coord_raster(ivar)=true;
+        end
+    end
+else
+    for ivar=find(check_coord_names)
+        DimIndex=find(strcmp(Data.VarDimName{ivar},Data.ListDimName));
+        if Data.DimValue(DimIndex)==2
+            check_coord_raster(ivar)=true;
+        end
+    end
+end
+
+
+%% initate cells around each scalar field
+index_remove=[];
+cell_nbre=numel(ind_scalar)+numel(ind_vector_x);
+flag_remove=false(1,cell_nbre);
+NbDim=zeros(1,cell_nbre);
+index_coord_x=zeros(size(ind_coord_x));
+for icell=1:numel(ind_scalar)
+    CellInfo{icell}.VarType='scalar';
+    CellInfo{icell}.VarIndex_scalar=ind_scalar(icell);
+    CellInfo{icell}.VarIndex=ind_scalar(icell);
+    DimCell_var=Data.VarDimName{ind_scalar(icell)};% cell of dimension names for ivar_coord_x(icell)
+    %look for errorflag
+    for ivar=ind_errorflag
+        DimCell=Data.VarDimName{ivar};
+        if isequal(DimCell,DimCell_var)
+            CellInfo{icell}.VarIndex(2)=ivar;
+            CellInfo{icell}.VarIndex_errorflag=ivar;
+            break
+        end
+    end
+end
+
+%% initate cells around each vector field
+for icell=numel(ind_scalar)+1:cell_nbre
+    CellInfo{icell}.VarType='vector';
+    CellInfo{icell}.VarIndex(1)=ind_vector_x(icell-numel(ind_scalar));
+    CellInfo{icell}.VarIndex_vector_x=ind_vector_x(icell-numel(ind_scalar));
+    DimCell_var=Data.VarDimName{ind_vector_x(icell-numel(ind_scalar))};% cell of dimension names for ivar_coord_x(icell)
+    % look for the associated y vector component
+    nbvar=1;
+    for ivar=ind_vector_y
+        DimCell=Data.VarDimName{ivar};
+        if isequal(DimCell,DimCell_var)
+            CellInfo{icell}.VarIndex(2)=ivar;
+            nbvar=2;
+            CellInfo{icell}.VarIndex_vector_y=ivar;
+            break
+        end
+    end
+    if ~isfield(CellInfo{icell},'VarIndex_vector_y')
+        flag_remove(icell)=true;% no vector_y found , mark cell to remove
+    end
+    % look for the associated z vector component
+    for ivar=ind_vector_z
+        DimCell=Data.VarDimName{ivar};
+        if isequal(DimCell,DimCell_var)
+            CellInfo{icell}.VarIndex(3)=ivar;
+            nbvar=3;
+            break
+        end
+    end
+    %look for the vector color scalar (ancillary)
+    for ivar=ind_ancillary
+        DimCell=Data.VarDimName{ivar};
+        if isequal(DimCell,DimCell_var)
+            nbvar=nbvar+1;
+            CellInfo{icell}.VarIndex(nbvar)=ivar;
+            CellInfo{icell}.VarIndex_ancillary=ivar;
+            break
+        end
+    end
+    %look for warnflag
+    for ivar=ind_warnflag
+        DimCell=Data.VarDimName{ivar};
+        if isequal(DimCell,DimCell_var)
+            nbvar=nbvar+1;
+            CellInfo{icell}.VarIndex(nbvar)=ivar;
+            CellInfo{icell}.VarIndex_warnflag=ivar;
+            break
+        end
+    end
+    %look for errorflag
+    for ivar=ind_errorflag
+        DimCell=Data.VarDimName{ivar};
+        if isequal(DimCell,DimCell_var)
+            nbvar=nbvar+1;
+            CellInfo{icell}.VarIndex(nbvar)=ivar;
+            CellInfo{icell}.VarIndex_errorflag=ivar;
+            break
+        end
+    end
+end
+
+%% find coordinates for each cell around field variables, scalars or vectors
+for icell=1:cell_nbre
+    CellInfo{icell}.CoordType='';
+    ind_var=CellInfo{icell}.VarIndex(1);
+    DimCell_var=Data.VarDimName{ind_var};% cell of dimension names for ivar_coord_x(icell)
+    if ~check_var
+        for idim=1:numel(DimCell_var)
+            CellInfo{icell}.DimIndex(idim)=find(strcmp(DimCell_var{idim},Data.ListDimName));
+        end
+    end
+    %look for z scattered coordinates
+    if isempty(ind_coord_z)
+        NbDim(icell)=2;
+        CellInfo{icell}.CoordIndex=[0 0];
+    else
+        NbDim(icell)=3;
+        CellInfo{icell}.CoordIndex=[0 0 0];
+        for ivar=ind_coord_z
+            DimCell=Data.VarDimName{ivar};
+            if isequal(DimCell,DimCell_var)
+                CellInfo{icell}.CoordType='scattered';
+                CellInfo{icell}.CoordIndex(1)=ivar;
+                CellInfo{icell}.ZName=Data.ListVarName{ivar};
+                CellInfo{icell}.ZIndex=ivar;
+                break
+            end
+        end
+    end
+    % look for y coordinate
+    for ivar=ind_coord_y
+        % detect scattered y coordinates, variable with the same dimension(s) as the field variable considered
+        DimCell=Data.VarDimName{ivar};
+        if isequal(DimCell,DimCell_var)
+            CellInfo{icell}.CoordType='scattered';
+            CellInfo{icell}.CoordIndex(NbDim(icell)-1)=ivar;
+            CellInfo{icell}.YName=Data.ListVarName{ivar};
+            CellInfo{icell}.YIndex=ivar;
+            break
+        end
+    end
+    
+    %look for x coordinates
+    if strcmp(CellInfo{icell}.CoordType,'scattered')
+        for ivar=ind_coord_x
+            DimCell=Data.VarDimName{ivar};
+            if isequal(DimCell,DimCell_var)
+                CellInfo{icell}.CoordIndex(NbDim(icell))=ivar;
+                CellInfo{icell}.XName=Data.ListVarName{ivar};
+                CellInfo{icell}.XIndex=ivar;
+                break
+            end
+        end
+    end
+    if isfield(CellInfo{icell},'ZName')
+        if isfield(CellInfo{icell},'YName')&& isfield(CellInfo{icell},'XName')
+            continue %scattered coordinates OK
+        end
+    else
+        if isfield(CellInfo{icell},'YName')
+            if isfield(CellInfo{icell},'XName')
+                NbDim(icell)=2;
+                continue %scattered coordinates OK
+            end
+        else
+            if isfield(CellInfo{icell},'XName'); % only one coordinate x, switch vector field to 1D plot
+                for ind=1:numel(CellInfo{icell}.VarIndex)
+                    Role{CellInfo{icell}.VarIndex(ind)}='coord_y';
+                end
+                continue
+            end
+        end
+    end
+    
+    %look for grid  coordinates
+    if isempty(CellInfo{icell}.CoordType)
+        NbDim(icell)=numel(DimCell_var);
+        CellInfo{icell}.DimOrder=[];
+        if NbDim(icell)==3
+            %coord z
+            for ivar=ind_coord_z
+                if check_coord_names(ivar)
+                    DimRank=find(strcmp(Data.VarDimName{ivar},DimCell_var));
+                check_coord=~isempty(DimRank);
+            elseif check_coord_raster(ivar)
+                DimRank=find(strcmp(Data.ListVarName{ivar},DimCell_var));
+                check_coord=~isempty(DimRank);
+            end
+%                 check_coord= (check_coord_names(ivar) && strcmp(Data.VarDimName{ivar},DimCell_var{1}))||...% coord varbable
+%                     (check_coord_raster(ivar) && strcmp(Data.ListVarName{ivar},DimCell_var{1})); % rasrewr coord defined by min and max
+                if check_coord
+                    CellInfo{icell}.CoordType='grid';
+                    CellInfo{icell}.CoordIndex(1)=ivar;
+                    CellInfo{icell}.ZName=Data.ListVarName{ivar};
+                    CellInfo{icell}.ZIndex=ivar;
+                    CellInfo{icell}.DimOrder=DimRank;
+                    break
+                end
+            end
+        end
+        for ivar=ind_coord_y
+              if check_coord_names(ivar)
+                    DimRank=find(strcmp(Data.VarDimName{ivar},DimCell_var));
+                check_coord=~isempty(DimRank);
+            elseif check_coord_raster(ivar)
+                DimRank=find(strcmp(Data.ListVarName{ivar},DimCell_var));
+                check_coord=~isempty(DimRank);
+            end
+%             check_coord= (check_coord_names(ivar) && strcmp(Data.VarDimName{ivar},DimCell_var{NbDim(icell)-1}))||...% coord variable
+%                 (check_coord_raster(ivar) && strcmp(Data.ListVarName{ivar},DimCell_var{NbDim(icell)-1})); % rasrewr coord defined by min and max
+            if check_coord
+                CellInfo{icell}.CoordType='grid';
+                CellInfo{icell}.CoordIndex(NbDim(icell)-1)=ivar;
+                CellInfo{icell}.YName=Data.ListVarName{ivar};
+                CellInfo{icell}.YIndex=ivar;
+                CellInfo{icell}.DimOrder=[CellInfo{icell}.DimOrder DimRank];
+                break
+            end
+        end
+        for ivar=ind_coord_x
+            if check_coord_names(ivar)
+                    DimRank=find(strcmp(Data.VarDimName{ivar},DimCell_var));
+                check_coord=~isempty(DimRank);
+            elseif check_coord_raster(ivar)
+                DimRank=find(strcmp(Data.ListVarName{ivar},DimCell_var));
+                check_coord=~isempty(DimRank);
+            end
+%             check_coord= (check_coord_names(ivar) && strcmp(Data.VarDimName{ivar},DimCell_var{NbDim(icell)}))||...% coord variable
+%                 (check_coord_raster(ivar) && strcmp(Data.ListVarName{ivar},DimCell_var{NbDim(icell)})); % raster coord defined by min and max
+            if check_coord
+                CellInfo{icell}.CoordIndex(NbDim(icell))=ivar;
+                CellInfo{icell}.XName=Data.ListVarName{ivar};
+                CellInfo{icell}.XIndex=ivar;
+                CellInfo{icell}.DimOrder=[CellInfo{icell}.DimOrder DimRank];
+                break
+            end
+        end
+    end
+    %look for tps coordinates
+    for ivar=ind_coord_tps
+        DimCell=Data.VarDimName{ivar};
+        if  numel(DimCell)==3 && strcmp(DimCell{1},DimCell_var{1})
+            CellInfo{icell}.CoordType='tps';
+            CellInfo{icell}.CoordIndex=ivar;
             if check_var
-                CellInfo{icell}.CoordSize=numel(Data.(Data.ListVarName{ivar_coord_x(icell)}));
+                NbDim(icell)=size(Data.(Data.ListVarName{ivar}),2);
             else
-                for idim=1:numel(DimCell)
-                    check_index= strcmp(DimCell{idim},Data.ListDimName);
-                    CellInfo{icell}.CoordSize(idim)=Data.DimValue(check_index);
+                DimIndex=find(strcmp(Data.VarDimName{ivar},Data.ListDimName));
+                NbDim(icell)= Data.DimValue(DimIndex);
+            end
+            for ivardim=1:numel(Data.VarDimName)
+                if strcmp(Data.VarDimName{ivardim},DimCell{3})
+                    CellInfo{icell}.NbCentres_tps= ivardim;% nbre of sites for each tps subdomain
+                elseif strcmp(Data.VarDimName{ivardim}{1},DimCell{2}) && numel(Data.VarDimName{ivardim})>=3 && strcmp(Data.VarDimName{ivardim}{3},DimCell{3})
+                    CellInfo{icell}.SubRange_tps=ivardim;% subrange definiton for tps
                 end
-                CellInfo{icell}.CoordSize=prod(CellInfo{icell}.CoordSize);
-            end
-            ind_scalar=find(strcmp('scalar',Role(VarIndex)));
-            ind_vector_x=find(strcmp('vector_x',Role(VarIndex)));
-            ind_vector_y=find(strcmp('vector_y',Role(VarIndex)));
-            ind_y=find(strcmp('coord_y',Role(VarIndex)));
-            if numel([ind_scalar ind_vector_x ind_vector_y])==0           
-%             if numel(VarIndex)==2||isempty(ind_y)% no variable, except possibly y
-                NbDim(icell)=1;
-            else
-                CellInfo{icell}.CoordType='scattered';
-                ind_z=find(strcmp('coord_z',Role(VarIndex)));
-                if numel(VarIndex)==3||isempty(ind_z)% no z variable, except possibly as a fct z(x,y)
-                    CellInfo{icell}.CoordIndex=[VarIndex(ind_y) CellInfo{icell}.CoordIndex];
-                    NbDim(icell)=2;
+            end
+        end
+        break
+    end
+end
+
+%% get number of coordinate points for each cell
+if check_var
+    for icell=1:numel(CellInfo)
+        switch CellInfo{icell}.CoordType
+            case 'scattered'
+                CellInfo{icell}.CoordSize=numel(Data.(CellInfo{icell}.XName));
+            case 'grid'
+                if NbDim(icell)==3
+                    CellInfo{icell}.CoordSize=[numel(Data.(CellInfo{icell}.XName)) numel(Data.(CellInfo{icell}.YName)) numel(Data.(CellInfo{icell}.YName))];
                 else
-                    CellInfo{icell}.CoordIndex=[VarIndex(ind_z) CellInfo{icell}.CoordIndex];
-                    NbDim(icell)=3;
+                    CellInfo{icell}.CoordSize=[numel(Data.(CellInfo{icell}.XName)) numel(Data.(CellInfo{icell}.YName))];
                 end
-            end
-        end
-        CellInfo{icell}.VarIndex=VarIndex;
-        check_select=check_select|check_cell;
-    end
-end
-
-%% look for tps coordinates
-ivar_remain=find(~check_select);% indices of remaining variables (not already selected)
-check_coord_tps= strcmp('coord_tps',Role(~check_select));
-ivar_tps=ivar_remain(check_coord_tps);% variable indices corresponding to tps coordinates
-
-% loop on the tps coordinate sets
-for icell_tps=1:numel(ivar_tps)
-    check_cell=zeros(1,nbvar);% =1 for the variables selected in the current cell
-    check_cell(ivar_tps(icell_tps))=1;% mark the coordinate variable as selected
-    DimCell=Data.VarDimName{ivar_tps(icell_tps)};% dimension names for the current tps coordinate variable
-    icell=numel(CellInfo)+icell_tps; % new field cell index
-    CellInfo{icell}.CoordIndex=ivar_tps(icell_tps);% index of the  tps coordinate variable
-    if numel(DimCell)==3
-        VarDimName=Data.VarDimName(~check_select);
-        for ivardim=1:numel(VarDimName)
-            if strcmp(VarDimName{ivardim},DimCell{3})
-                CellInfo{icell}.NbCentres_tps= ivar_remain(ivardim);% nbre of sites for each tps subdomain
-                check_cell(ivar_remain(ivardim))=1;% mark the variable as selected
-            elseif strcmp(VarDimName{ivardim}{1},DimCell{2}) && numel(VarDimName{ivardim})>=3 && strcmp(VarDimName{ivardim}{3},DimCell{3})
-                CellInfo{icell}.SubRange_tps=ivar_remain(ivardim);% subrange definiton for tps
-                check_cell(ivar_remain(ivardim))=1;% mark the variable as selected
-            elseif strcmp(VarDimName{ivardim}{1},DimCell{1}) && strcmp(VarDimName{ivardim}{2},DimCell{3})% variable
-                check_cell(ivar_remain(ivardim))=1;% mark the variable as selected
-            end
-        end
-    end
-    CellInfo{icell}.CoordType='tps';
-    CellInfo{icell}.VarIndex=find(check_cell);
-    if check_var
-        NbDim(icell)=size(Data.(Data.ListVarName{CellInfo{icell}.CoordIndex}),2);
-        CellInfo{icell}.CoordSize=size(Data.(Data.ListVarName{CellInfo{icell}.CoordIndex}),1)*size(Data.(Data.ListVarName{CellInfo{icell}.CoordIndex}),3);
-    else
-        check_index_1= strcmp(DimCell{1},Data.ListDimName);
-        check_index_2= strcmp(DimCell{2},Data.ListDimName);
-        NbDim(icell)=Data.DimValue(check_index_2);
-        if numel(DimCell)>=3
-        check_index_3= strcmp(DimCell{3},Data.ListDimName);  
-        CellInfo{icell}.CoordSize=Data.DimValue(check_index_1)*Data.DimValue(check_index_3);
-        end
-    end
-    check_select=check_select|check_cell;
-end
-
-%% look for coordinate variables and corresponding gridded data:
-% coordinate variables are variables associated with a single dimension, defining the coordinate values
-% two cases: 1)the coordiante variable represents the set of coordiante values
-%            2)the coordinate variable contains only two elements, representing the coordinate bounds for the dimension with the same name as the cordinate
-ivar_remain=find(~check_select);% indices of remaining variables, not already taken into account
-ListVarName=Data.ListVarName(~check_select);%list of names of remaining variables
-VarDimName=Data.VarDimName(~check_select);%dimensions of remaining variables
-check_coord_select= cellfun(@numel,VarDimName)==1|cellfun(@ischar,VarDimName)==1;% find remaining variables with a single dimension
-check_coord_select=check_coord_select & ~strcmp('ancillary',Role(~check_select));% do not select ancillary variables as coordinates
-ListCoordIndex=ivar_remain(check_coord_select);% indices of remaining variables with a single dimension
-ListCoordName=ListVarName(check_coord_select);% corresponding names of remaining variables with a single dimension
-ListDimName=VarDimName(check_coord_select);% dimension names of remaining variables with a single dimension
-
-%remove redondant variables -> keep only one variable per dimension
-check_keep=true(size(ListDimName));
-for idim=1:numel(ListDimName)
-    prev_ind=find(strcmp(ListDimName{idim},ListDimName(1:idim-1)));% check whether the dimension is already taken into account
-    if ~isempty(prev_ind)% in case of multiple coord variable
-        if strcmp(ListCoordName{idim},ListDimName{idim}) %variable with the same name as the coordinate taken in priority
-            check_keep(prev_ind)=0;% choose a variable with the same name as coordinate in priority
-        else
-           check_keep(idim)=0; %keep the first coordiante variable found
-        end
-    end
-end
-ListCoordIndex=ListCoordIndex(check_keep);% list of coordinate variable indices
-ListCoordName=ListCoordName(check_keep);% list of coordinate variable names
-ListDimName=ListDimName(check_keep);% list of coordinate dimension names
-
-% determine dimension sizes
-CoordSize=zeros(size(ListCoordIndex));
-for ilist=1:numel(ListCoordIndex)
-    if iscell(ListDimName{ilist})
-        ListDimName(ilist)=ListDimName{ilist};%transform cell to string
-    end
-    if check_var% if the list of dimensions has been directly defined, no variable data available
-        CoordSize(ilist)=numel(Data.(ListCoordName{ilist}));% number of elements in the variable corresponding to the dimension #ilist
-    else
-        check_index= strcmp(ListDimName{ilist},Data.ListDimName);% find the  index in the list of dimensions
-        CoordSize(ilist)=Data.DimValue(check_index);% find the  corresponding dimension value
-    end
-    if CoordSize(ilist)==2% case of uniform grid coordinate defined by lower and upper bounds only
-        ListDimName{ilist}=ListCoordName{ilist};% replace the dimension name by the coordinate variable name 
-    end
-end
-
-%% group the remaining variables in cells sharing the same coordinate variables
-NewCellInfo={};
-NewCellDimIndex={};
-NewNbDim=[];
-for ivardim=1:numel(VarDimName) % loop at the list of dimensions for the remaining variables
-    DimCell=VarDimName{ivardim};% dimension names of the current variable 
-    if ischar(DimCell), DimCell={DimCell}; end %transform char to cell if needed
-    DimIndices=[];
-    for idim=1:numel(DimCell)
-        ind_dim=find(strcmp(DimCell{idim},ListDimName));%find the dim index in the list of dimensions ListDimName
-        if ~isempty(ind_dim)
-            DimIndices=[DimIndices ind_dim]; %update the list of dim indices included in DimCell
-            if check_var & CoordSize(ind_dim)==2 % determine the size of the coordinate in case of coordinate definition limited to lower and upper bounds
-                if isvector(Data.(ListVarName{ivardim})) 
-                    if numel(Data.(ListVarName{ivardim}))>2
-                        CoordSize(ind_dim)=numel(Data.(ListVarName{ivardim}));
-                    end
-                else
-                    CoordSize(ind_dim)=size(Data.(ListVarName{ivardim}),idim);
-                end
-            end
-        end
-    end
-    % look for cells of variables with the same coordinate variables
-    check_previous=0;
-    for iprev=1:numel(NewCellInfo)
-        if isequal(DimIndices,NewCellDimIndex{iprev})
-            check_previous=1;
-            NewCellInfo{iprev}.VarIndex=[NewCellInfo{iprev}.VarIndex ivar_remain(ivardim)];%append the current variable index to the found field cell
-            break
-        end
-    end
-    % create a new cell if no previous one contains the coordinate variables
-    if ~check_previous
-        nbcell=numel(NewCellInfo)+1;
-        NewCellDimIndex{nbcell}=DimIndices;
-        NewCellInfo{nbcell}.VarIndex=ivar_remain(ivardim);% create a new field cell with the current variable index
-        NewNbDim(nbcell)=numel(DimIndices);
-        NewCellInfo{nbcell}.CoordType='grid';
-        NewCellInfo{nbcell}.CoordSize=CoordSize(DimIndices);
-        NewCellInfo{nbcell}.CoordIndex=ListCoordIndex(DimIndices);
-    end
-end
-NbDim=[NbDim NewNbDim];
-CellInfo=[CellInfo NewCellInfo];
-
-%% suppress empty cells or cells with a single coordinate variable 
-check_remove=false(size(CellInfo));
-for icell=1:numel(check_remove)
-    if isempty(CellInfo{icell})||(numel(CellInfo{icell}.VarIndex)==1 && numel(check_coord)>=icell && check_coord(icell))
-        check_remove(icell)=1;
-    end
-end
-CellInfo(check_remove)=[];
-NbDim(check_remove)=[];
+            case 'tps'
+                NbDim(icell)=size(Data.(Data.ListVarName{CellInfo{icell}.CoordIndex}),2);
+                CellInfo{icell}.CoordSize=size(Data.(Data.ListVarName{CellInfo{icell}.CoordIndex}),1);
+        end
+    end
+else
+    for icell=1:numel(CellInfo)
+        CellInfo{icell}.CoordSize=size(Data.DimValue(CellInfo{icell}.DimIndex));
+    end
+end
+%
+% %% loop on the tps coordinate sets
+%
+%     for icell_tps=1:numel(ind_coord_tps)
+%         check_cell=zeros(1,nbvar);% =1 for the variables selected in the current cell
+%         check_cell(ivar_tps(icell_tps))=1;% mark the coordinate variable as selected
+%         DimCell=Data.VarDimName{ivar_tps(icell_tps)};% dimension names for the current tps coordinate variable
+%         icell=numel(CellInfo)+icell_tps; % new field cell index
+%         CellInfo{icell}.CoordIndex=ivar_tps(icell_tps);% index of the  tps coordinate variable
+%         if numel(DimCell)==3
+%             VarDimName=Data.VarDimName(~check_select);
+%             for ivardim=1:numel(VarDimName)
+%                 if strcmp(VarDimName{ivardim},DimCell{3})
+%                     CellInfo{icell}.NbCentres_tps= ivar_remain(ivardim);% nbre of sites for each tps subdomain
+%                     check_cell(ivar_remain(ivardim))=1;% mark the variable as selected
+%                 elseif strcmp(VarDimName{ivardim}{1},DimCell{2}) && numel(VarDimName{ivardim})>=3 && strcmp(VarDimName{ivardim}{3},DimCell{3})
+%                     CellInfo{icell}.SubRange_tps=ivar_remain(ivardim);% subrange definiton for tps
+%                     check_cell(ivar_remain(ivardim))=1;% mark the variable as selected
+%                 elseif strcmp(VarDimName{ivardim}{1},DimCell{1}) && strcmp(VarDimName{ivardim}{2},DimCell{3})% variable
+%                     check_cell(ivar_remain(ivardim))=1;% mark the variable as selected
+%                 end
+%             end
+%         end
+%         CellInfo{icell}.CoordType='tps';
+%         CellInfo{icell}.VarIndex=find(check_cell);
+%         if check_var
+%             NbDim(icell)=size(Data.(Data.ListVarName{CellInfo{icell}.CoordIndex}),2);
+%             CellInfo{icell}.CoordSize=size(Data.(Data.ListVarName{CellInfo{icell}.CoordIndex}),1)*size(Data.(Data.ListVarName{CellInfo{icell}.CoordIndex}),3);
+%         else
+%             check_index_1= strcmp(DimCell{1},Data.ListDimName);
+%             check_index_2= strcmp(DimCell{2},Data.ListDimName);
+%             NbDim(icell)=Data.DimValue(check_index_2);
+%             if numel(DimCell)>=3
+%                 check_index_3= strcmp(DimCell{3},Data.ListDimName);
+%                 CellInfo{icell}.CoordSize=Data.DimValue(check_index_1)*Data.DimValue(check_index_3);
+%             end
+%         end
+%         check_select=check_select|check_cell;
+%     end
+
+
+%% cell for ordinary plots
+iremove=false(1,numel(ind_coord_y));
+for ilist=1:numel(ind_coord_y)% remove the y coordinates which have been used yet in scalar or vector fields
+    for icell=1:numel(CellInfo)
+        if isfield(CellInfo{icell},'YIndex')&& isequal(CellInfo{icell}.YIndex,ind_coord_y(ilist))
+            iremove(ilist)=true;
+            continue
+        end
+    end
+end
+ind_coord_y(iremove)=[];
+if ~isempty(ind_coord_x)
+    y_nbre=zeros(1,numel(ind_coord_x));
+    for icell=1:numel(ind_coord_x)
+        Cell1DPlot{icell}.VarType='1DPlot';
+        Cell1DPlot{icell}.XIndex=ind_coord_x(icell);
+        Cell1DPlot{icell}.XName=Data.ListVarName{ind_coord_x(icell)};
+        DimCell_x=Data.VarDimName{ind_coord_x(icell)};
+        for ivar=[ind_coord_y ind_discrete]
+            DimCell=Data.VarDimName{ivar};
+            if  numel(DimCell)==1 && strcmp(DimCell_x{1},DimCell{1})
+                y_nbre(icell)=y_nbre(icell)+1;
+                Cell1DPlot{icell}.YIndex(y_nbre(icell))=ivar;
+                break
+            end
+        end
+    end
+    Cell1DPlot(find(y_nbre==0))=[];
+    CellInfo=[CellInfo Cell1DPlot];
+    NbDim=[NbDim ones(1,numel(Cell1DPlot))];
+end
 
 %% document roles of non-coordinate variables
 for icell=1:numel(CellInfo)
-    VarIndex=CellInfo{icell}.VarIndex;
-    for ivar=VarIndex
-        if isfield(CellInfo{icell},['VarIndex_' Role{ivar}])
-            CellInfo{icell}.(['VarIndex_' Role{ivar}])=[CellInfo{icell}.(['VarIndex_' Role{ivar}]) ivar];
-        else
-            CellInfo{icell}.(['VarIndex_' Role{ivar}])= ivar;
-        end
-        if ~isempty(ProjModeRequest{ivar})
-            CellInfo{icell}.ProjModeRequest=ProjModeRequest{ivar};
-        end
-        if ~isempty(FieldName{ivar})
-            CellInfo{icell}.FieldName=FieldName{ivar};
-        end
-        if CheckSub(ivar)==1
-            CellInfo{icell}.CheckSub=1;
-        end
-    end
-end
+    if isfield(CellInfo{icell},'VarIndex')
+        VarIndex=CellInfo{icell}.VarIndex;
+        for ivar=VarIndex
+            %         if isfield(CellInfo{icell},['VarIndex_' Role{ivar}])
+            %             CellInfo{icell}.(['VarIndex_' Role{ivar}])=[CellInfo{icell}.(['VarIndex_' Role{ivar}]) ivar];
+            %         else
+            %             CellInfo{icell}.(['VarIndex_' Role{ivar}])= ivar;
+            %         end
+            if ~isempty(ProjModeRequest{ivar})
+                CellInfo{icell}.ProjModeRequest=ProjModeRequest{ivar};
+            end
+            if ~isempty(FieldName{ivar})
+                CellInfo{icell}.FieldName=FieldName{ivar};
+            end
+            if CheckSub(ivar)==1
+                CellInfo{icell}.CheckSub=1;
+            end
+        end
+    end
+end
+% for icell=ind_coord_tps
+%     VarIndex=CellInfo{icell}.VarIndex;
+%     for ivar=VarIndex
+%         if isfield(CellInfo{icell},['VarIndex_' Role{ivar}])
+%             CellInfo{icell}.(['VarIndex_' Role{ivar}])=[CellInfo{icell}.(['VarIndex_' Role{ivar}]) ivar];
+%         else
+%             CellInfo{icell}.(['VarIndex_' Role{ivar}])= ivar;
+%         end
+%         if ~isempty(ProjModeRequest{ivar})
+%             CellInfo{icell}.ProjModeRequest=ProjModeRequest{ivar};
+%         end
+%         if ~isempty(FieldName{ivar})
+%             CellInfo{icell}.FieldName=FieldName{ivar};
+%         end
+%         if CheckSub(ivar)==1
+%             CellInfo{icell}.CheckSub=1;
+%         end
+%     end
+% end
+
+
+
+% 
+% %% analyse vector fields
+% if ~isempty(ind_vector_x) && ~isempty(ind_vector_y)
+%     if numel(ind_vector_x)>1
+%         errormsg='multiply defined vector x component'
+%         return
+%     end
+%     DimCell_vec=Data.VarDimName{ind_vector_x};% cell of dimension names for ivar_coord_x(icell)
+%     if ischar(DimCell),DimCell={DimCell};end % transform char to cell for a single dimension
+%     DimCell_y=Data.VarDimName{ind_vector_y};% cell of dimension names for ivar_coord_x(icell)
+%     if ischar(DimCell_y),DimCell_y={DimCell_y};end % transform char to cell for a single dimension
+%     if ~isequal(DimCell,DimCell_y)
+%         errormsg='inconsistent x and y vector components';
+%         return
+%     end
+%     %look for coordinates
+%     for ivar=ind_coord_y
+%         DimCell=Data.VarDimName{ivar};
+%         if ischar(DimCell),DimCell={DimCell};end % transform char to cell for a single dimension
+%         if isequal(DimCell,DimCell_vec)
+%             CoordType='scattered';
+%             coordy=ivar;
+%         else
+%             if isempty(ind_coord_z) && strcmp(DimCell{1},DimCell_vec{1})
+%                 CoordType='grid';
+%                 coordy=ivar;
+%             elseif ~isempty(ind_coord_z) && strcmp(DimCell{1},DimCell_vec{2})
+%                 CoordType='grid';
+%                 coordy=ivar;
+%                 coordz=ind_coord_z;
+%             end
+%         end
+%         
+%         %% find scattered (unstructured) coordinates
+%         ivar_coord_x=find(strcmp('coord_x',Role));%find variables with Role='coord_x'
+%         check_select=false(1,nbvar);
+%         check_coord=false(1,nbvar);
+%         CellInfo=cell(1,numel(ivar_coord_x));
+%         NbDim=zeros(1,numel(ivar_coord_x));
+%         % loop on unstructured coordinate x -> different field cells
+%         for icell=1:numel(ivar_coord_x)
+%             DimCell=Data.VarDimName{ivar_coord_x(icell)};% cell of dimension names for ivar_coord_x(icell)
+%             if ischar(DimCell),DimCell={DimCell};end % transform char to cell for a single dimension
+%             % look for variables sharing dimension(s) with ivar_coord_x(icell)
+%             check_cell=zeros(numel(DimCell),nbvar);
+%             for idim=1:numel(DimCell);% for each variable with role coord_x, look at which other variables contain the same dimension
+%                 for ivar=1:nbvar
+%                     check_cell(idim,ivar)=max(strcmp(DimCell{idim},Data.VarDimName{ivar}));
+%                 end
+%             end
+%             check_cell=sum(check_cell,1)==numel(DimCell);%logical array=1 for variables belonging to the current cell
+%             VarIndex=find(check_cell);% list of detected variable indices
+%             if ~(numel(VarIndex)==1 && numel(DimCell)==1)% exclude case of isolated coord_x variable (treated later)
+%                 if ~(numel(VarIndex)==1 && numel(DimCell)>1)% a variable is associated to coordinate
+%                     CellInfo{icell}.CoordIndex=ivar_coord_x(icell);
+%                     % size of coordinate var
+%                     if check_var
+%                         CellInfo{icell}.CoordSize=numel(Data.(Data.ListVarName{ivar_coord_x(icell)}));
+%                     else
+%                         for idim=1:numel(DimCell)
+%                             check_index= strcmp(DimCell{idim},Data.ListDimName);
+%                             CellInfo{icell}.CoordSize(idim)=Data.DimValue(check_index);
+%                         end
+%                         CellInfo{icell}.CoordSize=prod(CellInfo{icell}.CoordSize);
+%                     end
+%                     %             ind_scalar=find(strcmp('scalar',Role(VarIndex)));
+%                     %             ind_vector_x=find(strcmp('vector_x',Role(VarIndex)));
+%                     %             ind_vector_y=find(strcmp('vector_y',Role(VarIndex)));
+%                     ind_y=find(strcmp('coord_y',Role(VarIndex)));
+%                     if numel([ind_scalar ind_vector_x ind_vector_y])==0
+%                         %             if numel(VarIndex)==2||isempty(ind_y)% no variable, except possibly y
+%                         NbDim(icell)=1;
+%                     else
+%                         CellInfo{icell}.CoordType='scattered';
+%                         ind_z=find(strcmp('coord_z',Role(VarIndex)));
+%                         if numel(VarIndex)==3||isempty(ind_z)% no z variable, except possibly as a fct z(x,y)
+%                             CellInfo{icell}.CoordIndex=[VarIndex(ind_y) CellInfo{icell}.CoordIndex];
+%                             NbDim(icell)=2;
+%                         else
+%                             CellInfo{icell}.CoordIndex=[VarIndex(ind_z) CellInfo{icell}.CoordIndex];
+%                             NbDim(icell)=3;
+%                         end
+%                     end
+%                 end
+%                 CellInfo{icell}.VarIndex=VarIndex;
+%                 check_select=check_select|check_cell;
+%             end
+%         end
+%         
+%         %% look for tps coordinates
+%         ivar_remain=find(~check_select);% indices of remaining variables (not already selected)
+%         check_coord_tps= strcmp('coord_tps',Role(~check_select));
+%         ivar_tps=ivar_remain(check_coord_tps);% variable indices corresponding to tps coordinates
+%         
+%         % loop on the tps coordinate sets
+%         for icell_tps=1:numel(ivar_tps)
+%             check_cell=zeros(1,nbvar);% =1 for the variables selected in the current cell
+%             check_cell(ivar_tps(icell_tps))=1;% mark the coordinate variable as selected
+%             DimCell=Data.VarDimName{ivar_tps(icell_tps)};% dimension names for the current tps coordinate variable
+%             icell=numel(CellInfo)+icell_tps; % new field cell index
+%             CellInfo{icell}.CoordIndex=ivar_tps(icell_tps);% index of the  tps coordinate variable
+%             if numel(DimCell)==3
+%                 VarDimName=Data.VarDimName(~check_select);
+%                 for ivardim=1:numel(VarDimName)
+%                     if strcmp(VarDimName{ivardim},DimCell{3})
+%                         CellInfo{icell}.NbCentres_tps= ivar_remain(ivardim);% nbre of sites for each tps subdomain
+%                         check_cell(ivar_remain(ivardim))=1;% mark the variable as selected
+%                     elseif strcmp(VarDimName{ivardim}{1},DimCell{2}) && numel(VarDimName{ivardim})>=3 && strcmp(VarDimName{ivardim}{3},DimCell{3})
+%                         CellInfo{icell}.SubRange_tps=ivar_remain(ivardim);% subrange definiton for tps
+%                         check_cell(ivar_remain(ivardim))=1;% mark the variable as selected
+%                     elseif strcmp(VarDimName{ivardim}{1},DimCell{1}) && strcmp(VarDimName{ivardim}{2},DimCell{3})% variable
+%                         check_cell(ivar_remain(ivardim))=1;% mark the variable as selected
+%                     end
+%                 end
+%             end
+%             CellInfo{icell}.CoordType='tps';
+%             CellInfo{icell}.VarIndex=find(check_cell);
+%             if check_var
+%                 NbDim(icell)=size(Data.(Data.ListVarName{CellInfo{icell}.CoordIndex}),2);
+%                 CellInfo{icell}.CoordSize=size(Data.(Data.ListVarName{CellInfo{icell}.CoordIndex}),1)*size(Data.(Data.ListVarName{CellInfo{icell}.CoordIndex}),3);
+%             else
+%                 check_index_1= strcmp(DimCell{1},Data.ListDimName);
+%                 check_index_2= strcmp(DimCell{2},Data.ListDimName);
+%                 NbDim(icell)=Data.DimValue(check_index_2);
+%                 if numel(DimCell)>=3
+%                     check_index_3= strcmp(DimCell{3},Data.ListDimName);
+%                     CellInfo{icell}.CoordSize=Data.DimValue(check_index_1)*Data.DimValue(check_index_3);
+%                 end
+%             end
+%             check_select=check_select|check_cell;
+%         end
+%         
+%      
+%         
+%         % determine dimension sizes
+%         CoordSize=zeros(size(ListCoordIndex));
+%         for ilist=1:numel(ListCoordIndex)
+%             if iscell(ListDimName{ilist})
+%                 ListDimName(ilist)=ListDimName{ilist};%transform cell to string
+%             end
+%             if check_var% if the list of dimensions has been directly defined, no variable data available
+%                 CoordSize(ilist)=numel(Data.(ListCoordName{ilist}));% number of elements in the variable corresponding to the dimension #ilist
+%             else
+%                 check_index= strcmp(ListDimName{ilist},Data.ListDimName);% find the  index in the list of dimensions
+%                 CoordSize(ilist)=Data.DimValue(check_index);% find the  corresponding dimension value
+%             end
+%             if CoordSize(ilist)==2% case of uniform grid coordinate defined by lower and upper bounds only
+%                 ListDimName{ilist}=ListCoordName{ilist};% replace the dimension name by the coordinate variable name
+%             end
+%         end
+%     end
+% end
+% 
+% 
+% 
+% %% suppress empty cells or cells with a single coordinate variable 
+% check_remove=false(size(CellInfo));
+% for icell=1:numel(check_remove)
+%     if isempty(CellInfo{icell})||(numel(CellInfo{icell}.VarIndex)==1 && numel(check_coord)>=icell && check_coord(icell))
+%         check_remove(icell)=1;
+%     end
+% end
+% CellInfo(check_remove)=[];
+% NbDim(check_remove)=[];
+% 
+% 
