Index: /XIOS/trunk/inputs/Regular_ugrid/context_atmosphere.xml
===================================================================
--- /XIOS/trunk/inputs/Regular_ugrid/context_atmosphere.xml	(revision 924)
+++ /XIOS/trunk/inputs/Regular_ugrid/context_atmosphere.xml	(revision 924)
@@ -0,0 +1,21 @@
+<!-- Context atm -->
+
+<context id="atmosphere">
+
+  <field_definition freq_op="1ts" enabled=".TRUE." operation="instant" >
+   <field id="temp" name="temp" long_name="Surface Temperature" unit="K" domain_ref="face" />
+  </field_definition>
+
+  <file_definition type="one_file" output_freq="1h" par_access="collective" enabled=".TRUE.">
+    <file id="output_cube_CF" name="output_cube_CF" >
+      <field field_ref="temp"/>
+    </file> 
+    <file id="output_cube_UGRID" name="output_cube_UGRID" convention="UGRID">
+      <field field_ref="temp"/>
+   </file>
+ </file_definition>
+
+  <domain_definition>
+   <domain id="face" name="mesh" nvertex="4"/>  
+ </domain_definition>
+ </context>
Index: /XIOS/trunk/inputs/Regular_ugrid/iodef.xml
===================================================================
--- /XIOS/trunk/inputs/Regular_ugrid/iodef.xml	(revision 924)
+++ /XIOS/trunk/inputs/Regular_ugrid/iodef.xml	(revision 924)
@@ -0,0 +1,20 @@
+<?xml version="1.0"?>
+<simulation>
+
+  <context id="atmosphere" src="./context_atmosphere.xml"/> 
+
+  <context id="xios">
+      <variable_definition>
+        <variable_group id="buffer">
+            <variable id="optimal_buffer_size" type="string">performance</variable>
+            <variable id="buffer_size_factor" type="double">1.0</variable>
+         </variable_group>
+
+        <variable_group id="parameters" >
+          <variable id="info_level" type="int">1</variable>
+          <variable id="print_file" type="bool">true</variable>
+        </variable_group>
+      </variable_definition>
+  </context>
+
+</simulation>
Index: /XIOS/trunk/src/client_client_dht_decl.cpp
===================================================================
--- /XIOS/trunk/src/client_client_dht_decl.cpp	(revision 923)
+++ /XIOS/trunk/src/client_client_dht_decl.cpp	(revision 924)
@@ -13,4 +13,5 @@
 
 template class CClientClientDHTTemplate<int>;
+template class CClientClientDHTTemplate<size_t>;
 template class CClientClientDHTTemplate<PairIntInt>;
 
Index: /XIOS/trunk/src/client_client_dht_template.hpp
===================================================================
--- /XIOS/trunk/src/client_client_dht_template.hpp	(revision 923)
+++ /XIOS/trunk/src/client_client_dht_template.hpp	(revision 924)
@@ -125,4 +125,5 @@
 
 typedef CClientClientDHTTemplate<int> CClientClientDHTInt;
+typedef CClientClientDHTTemplate<size_t> CClientClientDHTSizet;
 typedef CClientClientDHTTemplate<PairIntInt> CClientClientDHTPairIntInt;
 
Index: /XIOS/trunk/src/dht_auto_indexing.cpp
===================================================================
--- /XIOS/trunk/src/dht_auto_indexing.cpp	(revision 923)
+++ /XIOS/trunk/src/dht_auto_indexing.cpp	(revision 924)
@@ -12,47 +12,103 @@
 {
 
-/*!
-  Constructor with initial distribution information and the corresponding index
-  Each process holds a piece of hash value, which has no information (global index) associated.
-The constructor tries to assign this information to each process basing on their rank.
-The global index (information) are assigned across processes in the incremental manner.
-The process has lower rank will have low global index.
-  \param [in] hashValue Carray contains hash value
-  \param [in] clientIntraComm communicator of clients
-*/
+  CDHTAutoIndexing::~CDHTAutoIndexing()
+  {
+  }
 
-CDHTAutoIndexing::CDHTAutoIndexing(const CArray<size_t,1>& hashValue,
-                                   const MPI_Comm& clientIntraComm)
-  : CClientClientDHTTemplate<size_t>(clientIntraComm)
-{
-  Index2VectorInfoTypeMap indexToVecInfoMap;
-  indexToVecInfoMap.rehash(std::ceil(hashValue.size()/indexToVecInfoMap.max_load_factor()));
-  CArray<size_t,1>::const_iterator it = hashValue.begin(), ite = hashValue.end();
+  /*!
+    \param [in]
+    \param [in]
+  */
 
-  // Just use a dummy value to initialize
-  int nbLvl = this->getNbLevel();
-  for (; it != ite; ++it)
+  CDHTAutoIndexing::CDHTAutoIndexing(const CArray<size_t,1>& hashValue,
+                                     const MPI_Comm& clientIntraComm)
+    : CClientClientDHTTemplate<size_t>(clientIntraComm)
   {
-    indexToVecInfoMap[*it].resize(1);
-    indexToVecInfoMap[*it][0] = 0;
+
+    nbIndexOnProc_ = hashValue.size();
+    size_t nbIndexAccum;
+    MPI_Scan(&nbIndexOnProc_, &nbIndexAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, clientIntraComm);
+
+    // Broadcasting the total number of indexes
+    int rank, size;
+    MPI_Comm_rank(clientIntraComm, &rank);
+    MPI_Comm_size(clientIntraComm, &size);
+    if (rank == (size-1)) nbIndexesGlobal_ = nbIndexAccum;
+    MPI_Bcast(&nbIndexesGlobal_, 1, MPI_UNSIGNED_LONG, size-1, clientIntraComm);
+
+    CArray<size_t,1>::const_iterator itbIdx = hashValue.begin(), itIdx,
+                                     iteIdx = hashValue.end();
+
+    size_t idx = 0;
+    beginIndexOnProc_ = nbIndexAccum - nbIndexOnProc_;
+    globalIndex_.resize(nbIndexOnProc_);
+    for (itIdx = itbIdx; itIdx != iteIdx; ++itIdx)
+    {
+      globalIndex_[idx] = beginIndexOnProc_ + idx;
+      ++idx ;
+    }
   }
-  computeDistributedIndex(indexToVecInfoMap, clientIntraComm, nbLvl-1);
 
-  // Find out number of index on other proc
-  size_t nbIndexOnProc = index2InfoMapping_.size();
-  size_t nbIndexAccum;
+  /*!
+   * \fn void CDHTAutoIndexing::CDHTAutoIndexing(Index2VectorInfoTypeMap& hashInitMap, const MPI_Comm& clientIntraComm)
+   * Assigns a global index to unique input indexes.
+   * The returned map has unique indexes as a key and global indexes as a mapped value.
+   * \param [in] hashInitMap map<size_t, vector<size_t>> is a map of unique indexes.
+   * \param [in] clientIntraComm
+  */
+  CDHTAutoIndexing::CDHTAutoIndexing(Index2VectorInfoTypeMap& hashInitMap,
+                                     const MPI_Comm& clientIntraComm)
+    : CClientClientDHTTemplate<size_t>(clientIntraComm)
+  {
 
-  MPI_Scan(&nbIndexOnProc, &nbIndexAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, clientIntraComm);
+    nbIndexOnProc_ = hashInitMap.size();
+    size_t nbIndexAccum;
+    MPI_Scan(&nbIndexOnProc_, &nbIndexAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, clientIntraComm);
 
-  Index2VectorInfoTypeMap::iterator itbIdx = index2InfoMapping_.begin(), itIdx,
-                                    iteIdx = index2InfoMapping_.end();
-  size_t idx = 0;
-  size_t idxBegin = nbIndexAccum - nbIndexOnProc;
-  for (itIdx = itbIdx; itIdx != iteIdx; ++itIdx)
+    int rank, size;
+    MPI_Comm_rank(clientIntraComm, &rank);
+    MPI_Comm_size(clientIntraComm, &size);
+    if (rank == (size-1)) nbIndexesGlobal_ = nbIndexAccum;
+    MPI_Bcast(&nbIndexesGlobal_, 1, MPI_UNSIGNED_LONG, size-1, clientIntraComm);
+
+    Index2VectorInfoTypeMap::iterator itbIdx = hashInitMap.begin(), itIdx,
+                                      iteIdx = hashInitMap.end();
+    size_t idx = 0;
+    beginIndexOnProc_ = nbIndexAccum - nbIndexOnProc_;
+    globalIndex_.resize(nbIndexOnProc_);
+    for (itIdx = itbIdx; itIdx != iteIdx; ++itIdx)
+    {
+      (itIdx->second)[0] = beginIndexOnProc_ + idx;
+      globalIndex_[idx] = beginIndexOnProc_ + idx;
+      ++idx ;
+    }
+  }
+
+  /*!
+   * \fn size_t CDHTAutoIndexing::getNbIndexesGlobal() const
+   * Returns the total number of global indexes.
+  */
+  size_t CDHTAutoIndexing::getNbIndexesGlobal() const
   {
-    (itIdx->second)[0] = idxBegin + idx;
-    ++idx ;
+    return nbIndexesGlobal_;
   }
-}
+
+  /*!
+   * \fn size_t CDHTAutoIndexing::getIndexStart() const
+   * Returns the starting global index for a proc.
+  */
+  size_t CDHTAutoIndexing::getIndexStart() const
+  {
+    return beginIndexOnProc_;
+  }
+
+  /*!
+   * \fn size_t CDHTAutoIndexing::getIndexCount() const
+   * Returns the number of global indexes on a proc.
+  */
+  size_t CDHTAutoIndexing::getIndexCount() const
+  {
+    return nbIndexOnProc_;
+  }
 
 }
Index: /XIOS/trunk/src/dht_auto_indexing.hpp
===================================================================
--- /XIOS/trunk/src/dht_auto_indexing.hpp	(revision 923)
+++ /XIOS/trunk/src/dht_auto_indexing.hpp	(revision 924)
@@ -22,10 +22,24 @@
 class CDHTAutoIndexing: public CClientClientDHTTemplate<size_t>
 {
-public:
-  CDHTAutoIndexing(const CArray<size_t,1>& hashValue,
-                   const MPI_Comm& clientIntraComm);
+  public:
+
+    CDHTAutoIndexing(const CArray<size_t,1>& hashValue,
+                     const MPI_Comm& clientIntraComm);
+
+    CDHTAutoIndexing(Index2VectorInfoTypeMap& hashInitMap,
+                     const MPI_Comm& clientIntraComm);
+
+    size_t getNbIndexesGlobal() const;
+    size_t getIndexStart() const;
+    size_t getIndexCount() const;
 
     /** Default destructor */
     virtual ~CDHTAutoIndexing();
+
+  protected:
+    std::vector<size_t> globalIndex_;
+    size_t nbIndexesGlobal_;
+    size_t nbIndexOnProc_ ;
+    size_t beginIndexOnProc_ ;
 };
 
Index: /XIOS/trunk/src/io/nc4_data_output.cpp
===================================================================
--- /XIOS/trunk/src/io/nc4_data_output.cpp	(revision 923)
+++ /XIOS/trunk/src/io/nc4_data_output.cpp	(revision 924)
@@ -26,5 +26,5 @@
       CNc4DataOutput::CNc4DataOutput
          (const StdString & filename, bool exist, bool useClassicFormat, bool useCFConvention,
-          MPI_Comm comm_file,bool multifile, bool isCollective, const StdString& timeCounterName)
+          MPI_Comm comm_file, bool multifile, bool isCollective, const StdString& timeCounterName)
             : SuperClass()
             , SuperClassWriter(filename, exist, useClassicFormat, useCFConvention, &comm_file, multifile, timeCounterName)
@@ -449,7 +449,7 @@
       StdString domid = domain->getDomainOutputName();
       StdString domainName = domain->name;
-      domain->assignMesh(domainName);
-      domain->mesh->createMesh(domain->lonvalue_srv, domain->latvalue_srv, domain->bounds_lon_srv, domain->bounds_lat_srv);
-      //domain->mesh->createMeshEpsilon(server->intraComm, domain->lonvalue_srv, domain->latvalue_srv, domain->bounds_lon_srv, domain->bounds_lat_srv);
+      domain->assignMesh(domainName, domain->nvertex);
+      //domain->mesh->createMesh(domain->lonvalue_srv, domain->latvalue_srv, domain->bounds_lon_srv, domain->bounds_lat_srv);
+      domain->mesh->createMeshEpsilon(server->intraComm, domain->lonvalue_srv, domain->latvalue_srv, domain->bounds_lon_srv, domain->bounds_lat_srv);
 
       StdString node_x = domainName + "_node_x";
@@ -479,4 +479,6 @@
         SuperClassWriter::addVariable(domainName, NC_INT, dim0);
         SuperClassWriter::addAttribute("cf_role", StdString("mesh_topology"), &domainName);
+        SuperClassWriter::addAttribute("long_name", StdString("Topology data of 2D unstructured mesh"), &domainName);
+        SuperClassWriter::addAttribute("topology_dimension", 2, &domainName);
         SuperClassWriter::addAttribute("node_coordinates", node_x + " " + node_y, &domainName);
       }
@@ -498,9 +500,9 @@
                 SuperClassWriter::addVariable(node_x, NC_FLOAT, dim0);
                 SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &node_x);
-                SuperClassWriter::addAttribute("longname_name", StdString("Longitude of mesh nodes."), &node_x);
+                SuperClassWriter::addAttribute("long_name", StdString("Longitude of mesh nodes."), &node_x);
                 SuperClassWriter::addAttribute("units", StdString("degrees_east"), &node_x);
                 SuperClassWriter::addVariable(node_y, NC_FLOAT, dim0);
                 SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &node_y);
-                SuperClassWriter::addAttribute("longname_name", StdString("Latitude of mesh nodes."), &node_y);
+                SuperClassWriter::addAttribute("long_name", StdString("Latitude of mesh nodes."), &node_y);
                 SuperClassWriter::addAttribute("units", StdString("degrees_north"), &node_y);
               }
@@ -512,30 +514,70 @@
               if (!SuperClassWriter::varExist(node_x) || !SuperClassWriter::varExist(node_y))
               {
-                SuperClassWriter::addDimension(dimNode, domain->mesh->nbNodes);
+                SuperClassWriter::addDimension(dimNode, domain->mesh->nbNodesGlo);
                 dim0.clear();
                 dim0.push_back(dimNode);
                 SuperClassWriter::addVariable(node_x, NC_FLOAT, dim0);
                 SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &node_x);
-                SuperClassWriter::addAttribute("longname_name", StdString("Longitude of mesh nodes."), &node_x);
+                SuperClassWriter::addAttribute("long_name", StdString("Longitude of mesh nodes."), &node_x);
                 SuperClassWriter::addAttribute("units", StdString("degrees_east"), &node_x);
                 SuperClassWriter::addVariable(node_y, NC_FLOAT, dim0);
                 SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &node_y);
-                SuperClassWriter::addAttribute("longname_name", StdString("Latitude of mesh nodes."), &node_y);
+                SuperClassWriter::addAttribute("long_name", StdString("Latitude of mesh nodes."), &node_y);
+                SuperClassWriter::addAttribute("units", StdString("degrees_north"), &node_y);
+              }
+              SuperClassWriter::addAttribute("edge_node_connectivity", edge_nodes, &domainName);
+              SuperClassWriter::addAttribute("edge_coordinates", edge_x + " " + edge_y, &domainName);
+              SuperClassWriter::addDimension(dimEdge, domain->ni_glo);
+              dim0.clear();
+              dim0.push_back(dimEdge);
+              SuperClassWriter::addVariable(edge_x, NC_FLOAT, dim0);
+              SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &edge_x);
+              SuperClassWriter::addAttribute("long_name", StdString("Characteristic longitude of mesh edges."), &edge_x);
+              SuperClassWriter::addAttribute("units", StdString("degrees_east"), &edge_x);
+              SuperClassWriter::addVariable(edge_y, NC_FLOAT, dim0);
+              SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &edge_y);
+              SuperClassWriter::addAttribute("long_name", StdString("Characteristic latitude of mesh edges."), &edge_y);
+              SuperClassWriter::addAttribute("units", StdString("degrees_north"), &edge_y);
+              dim0.clear();
+              dim0.push_back(dimEdge);
+              dim0.push_back(dimTwo);
+              SuperClassWriter::addVariable(edge_nodes, NC_INT, dim0);
+              SuperClassWriter::addAttribute("cf_role", StdString("edge_node_connectivity"), &edge_nodes);
+              SuperClassWriter::addAttribute("long_name", StdString("Maps every edge/link to two nodes that it connects."), &edge_nodes);
+              SuperClassWriter::addAttribute("start_index", 0, &edge_nodes);
+            } // domain->nvertex == 2
+
+            // Adding faces, edges, and nodes, if edges and nodes have not been defined previously
+            if (domain->nvertex > 2)
+            {
+              // Nodes
+              if (!SuperClassWriter::varExist(node_x) || !SuperClassWriter::varExist(node_y))
+              {
+                SuperClassWriter::addDimension(dimNode, domain->mesh->nbNodesGlo);
+                dim0.clear();
+                dim0.push_back(dimNode);
+                SuperClassWriter::addVariable(node_x, NC_FLOAT, dim0);
+                SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &node_x);
+                SuperClassWriter::addAttribute("long_name", StdString("Longitude of mesh nodes."), &node_x);
+                SuperClassWriter::addAttribute("units", StdString("degrees_east"), &node_x);
+                SuperClassWriter::addVariable(node_y, NC_FLOAT, dim0);
+                SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &node_y);
+                SuperClassWriter::addAttribute("long_name", StdString("Latitude of mesh nodes."), &node_y);
                 SuperClassWriter::addAttribute("units", StdString("degrees_north"), &node_y);
               }
               if (!SuperClassWriter::varExist(edge_x) || !SuperClassWriter::varExist(edge_y))
               {
+                SuperClassWriter::addAttribute("edge_coordinates", edge_x + " " + edge_y, &domainName);
                 SuperClassWriter::addAttribute("edge_node_connectivity", edge_nodes, &domainName);
-                SuperClassWriter::addAttribute("edge_coordinates", edge_x + " " + edge_y, &domainName);
-                SuperClassWriter::addDimension(dimEdge, domain->mesh->nbEdges);
+                SuperClassWriter::addDimension(dimEdge, domain->mesh->nbEdgesGlo);
                 dim0.clear();
                 dim0.push_back(dimEdge);
                 SuperClassWriter::addVariable(edge_x, NC_FLOAT, dim0);
                 SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &edge_x);
-                SuperClassWriter::addAttribute("longname_name", StdString("Characteristic longitude of mesh edges."), &edge_x);
+                SuperClassWriter::addAttribute("long_name", StdString("Characteristic longitude of mesh edges."), &edge_x);
                 SuperClassWriter::addAttribute("units", StdString("degrees_east"), &edge_x);
                 SuperClassWriter::addVariable(edge_y, NC_FLOAT, dim0);
                 SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &edge_y);
-                SuperClassWriter::addAttribute("longname_name", StdString("Characteristic latitude of mesh edges."), &edge_y);
+                SuperClassWriter::addAttribute("long_name", StdString("Characteristic latitude of mesh edges."), &edge_y);
                 SuperClassWriter::addAttribute("units", StdString("degrees_north"), &edge_y);
                 dim0.clear();
@@ -547,109 +589,66 @@
                 SuperClassWriter::addAttribute("start_index", 0, &edge_nodes);
               }
-            } // domain->nvertex == 2
-
-            // Adding faces, edges, and nodes, if edges and nodes have not been defined previously
-            if (domain->nvertex > 2)
-            {
-              // Nodes
-              if (!SuperClassWriter::varExist(node_x) || !SuperClassWriter::varExist(node_y))
-              {
-                SuperClassWriter::addDimension(dimNode, domain->mesh->nbNodes);
-                dim0.clear();
-                dim0.push_back(dimNode);
-                SuperClassWriter::addVariable(node_x, NC_FLOAT, dim0);
-                SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &node_x);
-                SuperClassWriter::addAttribute("longname_name", StdString("Longitude of mesh nodes."), &node_x);
-                SuperClassWriter::addAttribute("units", StdString("degrees_east"), &node_x);
-                SuperClassWriter::addVariable(node_y, NC_FLOAT, dim0);
-                SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &node_y);
-                SuperClassWriter::addAttribute("longname_name", StdString("Latitude of mesh nodes."), &node_y);
-                SuperClassWriter::addAttribute("units", StdString("degrees_north"), &node_y);
-              }
-              if (!SuperClassWriter::varExist(edge_x) || !SuperClassWriter::varExist(edge_y))
-              {
-                SuperClassWriter::addAttribute("edge_coordinates", edge_x + " " + edge_y, &domainName);
-                SuperClassWriter::addAttribute("edge_node_connectivity", edge_nodes, &domainName);
-                SuperClassWriter::addDimension(dimEdge, domain->mesh->nbEdges);
-                dim0.clear();
-                dim0.push_back(dimEdge);
-                SuperClassWriter::addVariable(edge_x, NC_FLOAT, dim0);
-                SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &edge_x);
-                SuperClassWriter::addAttribute("longname_name", StdString("Characteristic longitude of mesh edges."), &edge_x);
-                SuperClassWriter::addAttribute("units", StdString("degrees_east"), &edge_x);
-                SuperClassWriter::addVariable(edge_y, NC_FLOAT, dim0);
-                SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &edge_y);
-                SuperClassWriter::addAttribute("longname_name", StdString("Characteristic latitude of mesh edges."), &edge_y);
-                SuperClassWriter::addAttribute("units", StdString("degrees_north"), &edge_y);
-                dim0.clear();
-                dim0.push_back(dimEdge);
-                dim0.push_back(dimTwo);
-                SuperClassWriter::addVariable(edge_nodes, NC_INT, dim0);
-                SuperClassWriter::addAttribute("cf_role", StdString("edge_node_connectivity"), &edge_nodes);
-                SuperClassWriter::addAttribute("long_name", StdString("Maps every edge/link to two nodes that it connects."), &edge_nodes);
-                SuperClassWriter::addAttribute("start_index", 0, &edge_nodes);
-              }
-                SuperClassWriter::addAttribute("face_coordinates", face_x + " " + face_y, &domainName);
-                SuperClassWriter::addAttribute("face_node_connectivity", face_nodes, &domainName);
-                SuperClassWriter::addDimension(dimFace, domain->mesh->nbFaces);
-                SuperClassWriter::addDimension(dimVertex, domain->mesh->nvertex);
-                dim0.clear();
-                dim0.push_back(dimFace);
-                SuperClassWriter::addVariable(face_x, NC_FLOAT, dim0);
-                SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &face_x);
-                SuperClassWriter::addAttribute("longname_name", StdString("Characteristic longitude of mesh faces."), &face_x);
-                SuperClassWriter::addAttribute("units", StdString("degrees_east"), &face_x);
-                SuperClassWriter::addVariable(face_y, NC_FLOAT, dim0);
-                SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &face_y);
-                SuperClassWriter::addAttribute("longname_name", StdString("Characteristic latitude of mesh faces."), &face_y);
-                SuperClassWriter::addAttribute("units", StdString("degrees_north"), &face_y);
-                dim0.clear();
-                dim0.push_back(dimFace);
-                dim0.push_back(dimVertex);
-                SuperClassWriter::addVariable(face_nodes, NC_INT, dim0);
-                SuperClassWriter::addAttribute("cf_role", StdString("face_node_connectivity"), &face_nodes);
-                SuperClassWriter::addAttribute("long_name", StdString("Maps every face to its corner nodes."), &face_nodes);
-                SuperClassWriter::addAttribute("start_index", 0, &face_nodes);
-                dim0.clear();
-                dim0.push_back(dimFace);
-                dim0.push_back(dimVertex);
-                SuperClassWriter::addVariable(face_edges, NC_INT, dim0);
-                SuperClassWriter::addAttribute("cf_role", StdString("face_edge_connectivity"), &face_edges);
-                SuperClassWriter::addAttribute("long_name", StdString("Maps every face to its edges."), &face_edges);
-                SuperClassWriter::addAttribute("start_index", 0, &face_edges);
-                dim0.clear();
-                dim0.push_back(dimEdge);
-                dim0.push_back(dimTwo);
-                SuperClassWriter::addVariable(edge_faces, NC_INT, dim0);
-                SuperClassWriter::addAttribute("cf_role", StdString("edge_face connectivity"), &edge_faces);
-                SuperClassWriter::addAttribute("long_name", StdString("neighbor faces for edges"), &edge_faces);
-                SuperClassWriter::addAttribute("start_index", 0, &edge_faces);
-                SuperClassWriter::addAttribute("_FillValue", -999, &edge_faces);
-                SuperClassWriter::addAttribute("comment", StdString("missing neighbor faces are indicated using _FillValue"), &edge_faces);
-                dim0.clear();
-                dim0.push_back(dimFace);
-                dim0.push_back(dimVertex);
-                SuperClassWriter::addVariable(face_faces, NC_INT, dim0);
-                SuperClassWriter::addAttribute("cf_role", StdString("face_face connectivity"), &face_faces);
-                SuperClassWriter::addAttribute("long_name", StdString("Indicates which other faces neighbor each face"), &face_faces);
-                SuperClassWriter::addAttribute("start_index", 0, &face_faces);
-                SuperClassWriter::addAttribute("flag_values", -1, &face_faces);
-                SuperClassWriter::addAttribute("flag_meanings", StdString("out_of_mesh"), &face_faces);
+              SuperClassWriter::addAttribute("face_coordinates", face_x + " " + face_y, &domainName);
+              SuperClassWriter::addAttribute("face_node_connectivity", face_nodes, &domainName);
+              SuperClassWriter::addDimension(dimFace, domain->ni_glo);
+              SuperClassWriter::addDimension(dimVertex, domain->nvertex);
+              dim0.clear();
+              dim0.push_back(dimFace);
+              SuperClassWriter::addVariable(face_x, NC_FLOAT, dim0);
+              SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &face_x);
+              SuperClassWriter::addAttribute("long_name", StdString("Characteristic longitude of mesh faces."), &face_x);
+              SuperClassWriter::addAttribute("units", StdString("degrees_east"), &face_x);
+              SuperClassWriter::addVariable(face_y, NC_FLOAT, dim0);
+              SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &face_y);
+              SuperClassWriter::addAttribute("long_name", StdString("Characteristic latitude of mesh faces."), &face_y);
+              SuperClassWriter::addAttribute("units", StdString("degrees_north"), &face_y);
+              dim0.clear();
+              dim0.push_back(dimFace);
+              dim0.push_back(dimVertex);
+              SuperClassWriter::addVariable(face_nodes, NC_INT, dim0);
+              SuperClassWriter::addAttribute("cf_role", StdString("face_node_connectivity"), &face_nodes);
+              SuperClassWriter::addAttribute("long_name", StdString("Maps every face to its corner nodes."), &face_nodes);
+              SuperClassWriter::addAttribute("start_index", 0, &face_nodes);
+              dim0.clear();
+              dim0.push_back(dimFace);
+              dim0.push_back(dimVertex);
+              SuperClassWriter::addVariable(face_edges, NC_INT, dim0);
+              SuperClassWriter::addAttribute("cf_role", StdString("face_edge_connectivity"), &face_edges);
+              SuperClassWriter::addAttribute("long_name", StdString("Maps every face to its edges."), &face_edges);
+              SuperClassWriter::addAttribute("start_index", 0, &face_edges);
+              dim0.clear();
+              dim0.push_back(dimEdge);
+              dim0.push_back(dimTwo);
+              SuperClassWriter::addVariable(edge_faces, NC_INT, dim0);
+              SuperClassWriter::addAttribute("cf_role", StdString("edge_face connectivity"), &edge_faces);
+              SuperClassWriter::addAttribute("long_name", StdString("neighbor faces for edges"), &edge_faces);
+              SuperClassWriter::addAttribute("start_index", 0, &edge_faces);
+              SuperClassWriter::addAttribute("_FillValue", -999, &edge_faces);
+              SuperClassWriter::addAttribute("comment", StdString("missing neighbor faces are indicated using _FillValue"), &edge_faces);
+              dim0.clear();
+              dim0.push_back(dimFace);
+              dim0.push_back(dimVertex);
+              SuperClassWriter::addVariable(face_faces, NC_INT, dim0);
+              SuperClassWriter::addAttribute("cf_role", StdString("face_face connectivity"), &face_faces);
+              SuperClassWriter::addAttribute("long_name", StdString("Indicates which other faces neighbor each face"), &face_faces);
+              SuperClassWriter::addAttribute("start_index", 0, &face_faces);
+              SuperClassWriter::addAttribute("flag_values", -1, &face_faces);
+              SuperClassWriter::addAttribute("flag_meanings", StdString("out_of_mesh"), &face_faces);
             } // domain->nvertex > 2
 
             SuperClassWriter::definition_end();
 
-            std::vector<StdSize> start(1) ;
-            std::vector<StdSize> count(1) ;
-            if (domain->isEmpty())
-             {
-               start[0]=0 ;
-               count[0]=0 ;
-             }
-             else
-             {
-               start[0]=domain->zoom_ibegin_srv-domain->global_zoom_ibegin;
-               count[0]=domain->zoom_ni_srv ;
-             }
+            std::vector<StdSize> startEdges(1) ;
+            std::vector<StdSize> countEdges(1) ;
+            std::vector<StdSize> startNodes(1) ;
+            std::vector<StdSize> countNodes(1) ;
+            std::vector<StdSize> startFaces(1) ;
+            std::vector<StdSize> countFaces(1) ;
+            std::vector<StdSize> startEdgeNodes(2) ;
+            std::vector<StdSize> countEdgeNodes(2) ;
+            std::vector<StdSize> startEdgeFaces(2) ;
+            std::vector<StdSize> countEdgeFaces(2) ;
+            std::vector<StdSize> startFaceConctv(2) ;
+            std::vector<StdSize> countFaceConctv(2) ;
 
             if (!isWrittenDomain(domid))
@@ -657,28 +656,100 @@
               if (domain->nvertex == 1)
               {
-                SuperClassWriter::writeData(domain->mesh->node_lat, node_y, isCollective, 0, &start, &count);
-                SuperClassWriter::writeData(domain->mesh->node_lon, node_x, isCollective, 0, &start, &count);
+                if (domain->isEmpty())
+                 {
+                   startNodes[0]=0 ;
+                   countNodes[0]=0 ;
+                 }
+                 else
+                 {
+                   startNodes[0] = domain->zoom_ibegin_srv-domain->global_zoom_ibegin;
+                   countNodes[0] = domain->zoom_ni_srv ;
+                 }
+
+                SuperClassWriter::writeData(domain->mesh->node_lat, node_y, isCollective, 0, &startNodes, &countNodes);
+                SuperClassWriter::writeData(domain->mesh->node_lon, node_x, isCollective, 0, &startNodes, &countNodes);
               }
               else if (domain->nvertex == 2)
-                            {
-                SuperClassWriter::writeData(domain->mesh->node_lat, node_y, isCollective, 0);
-                SuperClassWriter::writeData(domain->mesh->node_lon, node_x, isCollective, 0);
-                SuperClassWriter::writeData(domain->mesh->edge_lat, edge_y, isCollective, 0);
-                SuperClassWriter::writeData(domain->mesh->edge_lon, edge_x, isCollective, 0);
-                SuperClassWriter::writeData(domain->mesh->edge_nodes, edge_nodes);
-                            }
+              {
+                if (domain->isEmpty())
+                 {
+                  startEdges[0]=0 ;
+                  countEdges[0]=0 ;
+                  startNodes[0]=0 ;
+                  countNodes[0]=0 ;
+                  startEdgeNodes[0]=0;
+                  startEdgeNodes[1]=0;
+                  countEdgeNodes[0]=0;
+                  countEdgeNodes[1]=0;
+
+                 }
+                 else
+                 {
+                   startEdges[0] = domain->zoom_ibegin_srv-domain->global_zoom_ibegin;
+                   countEdges[0] = domain->zoom_ni_srv ;
+                   startNodes[0] = domain->mesh->node_start;
+                   countNodes[0] = domain->mesh->node_count;
+                   startEdgeNodes[0] = domain->zoom_ibegin_srv-domain->global_zoom_ibegin;
+                   startEdgeNodes[1] = 0;
+                   countEdgeNodes[0] = domain->zoom_ni_srv;
+                   countEdgeNodes[1]= 2;
+                 }
+                SuperClassWriter::writeData(domain->mesh->node_lat, node_y, isCollective, 0, &startNodes, &countNodes);
+                SuperClassWriter::writeData(domain->mesh->node_lon, node_x, isCollective, 0, &startNodes, &countNodes);
+                SuperClassWriter::writeData(domain->mesh->edge_lat, edge_y, isCollective, 0, &startEdges, &countEdges);
+                SuperClassWriter::writeData(domain->mesh->edge_lon, edge_x, isCollective, 0, &startEdges, &countEdges);
+                SuperClassWriter::writeData(domain->mesh->edge_nodes, edge_nodes, isCollective, 0, &startEdgeNodes, &countEdgeNodes);
+              }
               else
               {
-                SuperClassWriter::writeData(domain->mesh->node_lat, node_y, isCollective, 0);
-                SuperClassWriter::writeData(domain->mesh->node_lon, node_x, isCollective, 0);
-                SuperClassWriter::writeData(domain->mesh->edge_lat, edge_y, isCollective, 0);
-                SuperClassWriter::writeData(domain->mesh->edge_lon, edge_x, isCollective, 0);
-                SuperClassWriter::writeData(domain->mesh->edge_nodes, edge_nodes);
-                SuperClassWriter::writeData(domain->mesh->face_lat, face_y, isCollective, 0);
-                SuperClassWriter::writeData(domain->mesh->face_lon, face_x, isCollective, 0);
-                SuperClassWriter::writeData(domain->mesh->face_nodes, face_nodes);
-                SuperClassWriter::writeData(domain->mesh->face_edges, face_edges);
-                SuperClassWriter::writeData(domain->mesh->edge_faces, edge_faces);
-                SuperClassWriter::writeData(domain->mesh->face_faces, face_faces);
+                if (domain->isEmpty())
+                 {
+                   startFaces[0] = 0 ;
+                   countFaces[0] = 0 ;
+                   startNodes[0] = 0;
+                   countNodes[0] = 0;
+                   startEdges[0] = 0;
+                   countEdges[0] = 0;
+                   startEdgeFaces[0] = 0;
+                   startEdgeFaces[1] = 0;
+                   countEdgeFaces[0] = 0;
+                   countEdgeFaces[1] = 0;
+                   startFaceConctv[0] = 0;
+                   startFaceConctv[1] = 0;
+                   countFaceConctv[0] = 0;
+                   countFaceConctv[1] = 0;
+                 }
+                 else
+                 {
+                   startFaces[0] = domain->zoom_ibegin_srv-domain->global_zoom_ibegin;
+                   countFaces[0] = domain->zoom_ni_srv ;
+                   startNodes[0] = domain->mesh->node_start;
+                   countNodes[0] = domain->mesh->node_count;
+                   startEdges[0] = domain->mesh->edge_start;
+                   countEdges[0] = domain->mesh->edge_count;
+                   startEdgeNodes[0] = domain->mesh->edge_start;
+                   startEdgeNodes[1] = 0;
+                   countEdgeNodes[0] = domain->mesh->edge_count;
+                   countEdgeNodes[1]= 2;
+                   startEdgeFaces[0] = domain->mesh->edge_start;
+                   startEdgeFaces[1]= 0;
+                   countEdgeFaces[0] = domain->mesh->edge_count;
+                   countEdgeFaces[1]= 2;
+                   startFaceConctv[0] = domain->zoom_ibegin_srv-domain->global_zoom_ibegin;
+                   startFaceConctv[1] = 0;
+                   countFaceConctv[0] = domain->zoom_ni_srv;
+                   countFaceConctv[1] = domain->nvertex;
+                 }
+                SuperClassWriter::writeData(domain->mesh->node_lat, node_y, isCollective, 0, &startNodes, &countNodes);
+                SuperClassWriter::writeData(domain->mesh->node_lon, node_x, isCollective, 0, &startNodes, &countNodes);
+                SuperClassWriter::writeData(domain->mesh->edge_lat, edge_y, isCollective, 0, &startEdges, &countEdges);
+                SuperClassWriter::writeData(domain->mesh->edge_lon, edge_x, isCollective, 0, &startEdges, &countEdges);
+                SuperClassWriter::writeData(domain->mesh->edge_nodes, edge_nodes, isCollective, 0, &startEdgeNodes, &countEdgeNodes);
+                SuperClassWriter::writeData(domain->mesh->face_lat, face_y, isCollective, 0, &startFaces, &countFaces);
+                SuperClassWriter::writeData(domain->mesh->face_lon, face_x, isCollective, 0, &startFaces, &countFaces);
+                SuperClassWriter::writeData(domain->mesh->face_nodes, face_nodes, isCollective, 0, &startFaceConctv, &countFaceConctv);
+                SuperClassWriter::writeData(domain->mesh->face_edges, face_edges, isCollective, 0, &startFaceConctv, &countFaceConctv);
+                SuperClassWriter::writeData(domain->mesh->edge_faces, edge_faces, isCollective, 0, &startEdgeFaces, &countEdgeFaces);
+                SuperClassWriter::writeData(domain->mesh->face_faces, face_faces, isCollective, 0, &startFaceConctv, &countFaceConctv);
               }
               setWrittenDomain(domid);
@@ -686,45 +757,47 @@
             else
             {
-              if (domain->nvertex == 1)
-              {
-                if ( (!domain->mesh->edgesAreWritten) && (!domain->mesh->facesAreWritten) )
-                {
-                  SuperClassWriter::writeData(domain->mesh->node_lat, node_y, isCollective, 0);
-                  SuperClassWriter::writeData(domain->mesh->node_lon, node_x, isCollective, 0);
-                }
-              }
               if (domain->nvertex == 2)
               {
-                if (!domain->mesh->facesAreWritten)
-                {
-                  if (!domain->mesh->nodesAreWritten)
-                  {
-                    SuperClassWriter::writeData(domain->mesh->node_lat, node_y, isCollective, 0);
-                    SuperClassWriter::writeData(domain->mesh->node_lon, node_x, isCollective, 0);
-                  }
-                  SuperClassWriter::writeData(domain->mesh->edge_lat, edge_y, isCollective, 0);
-                  SuperClassWriter::writeData(domain->mesh->edge_lon, edge_x, isCollective, 0);
-                  SuperClassWriter::writeData(domain->mesh->edge_nodes, edge_nodes);
-                }
+                startEdges[0] = domain->zoom_ibegin_srv-domain->global_zoom_ibegin;
+                countEdges[0] = domain->zoom_ni_srv ;
+                startEdgeNodes[0] = domain->zoom_ibegin_srv-domain->global_zoom_ibegin;
+                startEdgeNodes[1] = 0;
+                countEdgeNodes[0] = domain->zoom_ni_srv;
+                countEdgeNodes[1]= 2;
+                SuperClassWriter::writeData(domain->mesh->edge_lat, edge_y, isCollective, 0, &startEdges, &countEdges);
+                SuperClassWriter::writeData(domain->mesh->edge_lon, edge_x, isCollective, 0, &startEdges, &countEdges);
+                SuperClassWriter::writeData(domain->mesh->edge_nodes, edge_nodes, isCollective, 0, &startEdgeNodes, &countEdgeNodes);
               }  
               if (domain->nvertex > 2)
               {
+
                 if (!domain->mesh->edgesAreWritten)
                 {
-                  if (!domain->mesh->nodesAreWritten)
-                  {
-                    SuperClassWriter::writeData(domain->mesh->node_lat, node_y, isCollective, 0);
-                    SuperClassWriter::writeData(domain->mesh->node_lon, node_x, isCollective, 0);
-                  }
-                  SuperClassWriter::writeData(domain->mesh->edge_lat, edge_y, isCollective, 0);
-                  SuperClassWriter::writeData(domain->mesh->edge_lon, edge_x, isCollective, 0);
-                  SuperClassWriter::writeData(domain->mesh->edge_nodes, edge_nodes);
+                  startEdges[0] = domain->mesh->edge_start;
+                  countEdges[0] = domain->mesh->edge_count;
+                  startEdgeNodes[0] = domain->mesh->edge_start;
+                  startEdgeNodes[1] = 0;
+                  countEdgeNodes[0] = domain->mesh->edge_count;
+                  countEdgeNodes[1]= 2;
+                  SuperClassWriter::writeData(domain->mesh->edge_lat, edge_y, isCollective, 0, &startEdges, &countEdges);
+                  SuperClassWriter::writeData(domain->mesh->edge_lon, edge_x, isCollective, 0, &startEdges, &countEdges);
+                  SuperClassWriter::writeData(domain->mesh->edge_nodes, edge_nodes, isCollective, 0, &startEdgeNodes, &countEdgeNodes);
                 }
-                SuperClassWriter::writeData(domain->mesh->face_lat, face_y, isCollective, 0);
-                SuperClassWriter::writeData(domain->mesh->face_lon, face_x, isCollective, 0);
-                SuperClassWriter::writeData(domain->mesh->face_nodes, face_nodes);
-                SuperClassWriter::writeData(domain->mesh->face_edges, face_edges);
-                SuperClassWriter::writeData(domain->mesh->edge_faces, edge_faces);
-                SuperClassWriter::writeData(domain->mesh->face_faces, face_faces);
+                startFaces[0] = domain->zoom_ibegin_srv-domain->global_zoom_ibegin;
+                countFaces[0] = domain->zoom_ni_srv;
+                startEdgeFaces[0] = domain->mesh->edge_start;
+                startEdgeFaces[1]= 0;
+                countEdgeFaces[0] = domain->mesh->edge_count;
+                countEdgeFaces[1]= 2;
+                startFaceConctv[0] = domain->zoom_ibegin_srv-domain->global_zoom_ibegin;
+                startFaceConctv[1] = 0;
+                countFaceConctv[0] = domain->zoom_ni_srv;
+                countFaceConctv[1]= domain->nvertex;
+                SuperClassWriter::writeData(domain->mesh->face_lat, face_y, isCollective, 0, &startFaces, &countFaces);
+                SuperClassWriter::writeData(domain->mesh->face_lon, face_x, isCollective, 0, &startFaces, &countFaces);
+                SuperClassWriter::writeData(domain->mesh->face_nodes, face_nodes, isCollective, 0, &startFaceConctv, &countFaceConctv);
+                SuperClassWriter::writeData(domain->mesh->face_edges, face_edges, isCollective, 0, &startFaceConctv, &countFaceConctv);
+                SuperClassWriter::writeData(domain->mesh->edge_faces, edge_faces, isCollective, 0, &startEdgeFaces, &countEdgeFaces);
+                SuperClassWriter::writeData(domain->mesh->face_faces, face_faces, isCollective, 0, &startFaceConctv, &countFaceConctv);
               }
             }// isWrittenDomain
@@ -754,5 +827,5 @@
           msg.append(context->getId()); msg.append("\n");
           msg.append(e.what());
-          ERROR("CNc4DataOutput::writeUnstructuredDomain(CDomain* domain)", << msg);
+          ERROR("CNc4DataOutput::writeUnstructuredDomainUgrid(CDomain* domain)", << msg);
         }
 
Index: /XIOS/trunk/src/io/onetcdf4_decl.cpp
===================================================================
--- /XIOS/trunk/src/io/onetcdf4_decl.cpp	(revision 923)
+++ /XIOS/trunk/src/io/onetcdf4_decl.cpp	(revision 924)
@@ -10,4 +10,5 @@
  
   macro(int, 1)
+  macro(int, 2)
   macro(double, 1)
   macro(double, 2)
Index: /XIOS/trunk/src/node/domain.cpp
===================================================================
--- /XIOS/trunk/src/node/domain.cpp	(revision 923)
+++ /XIOS/trunk/src/node/domain.cpp	(revision 924)
@@ -36,6 +36,5 @@
       , isRedistributed_(false)
    {
-	   //mesh = new CMesh();
-	}
+   }
 
    CDomain::CDomain(const StdString & id)
@@ -47,17 +46,15 @@
       , isRedistributed_(false)
    {
-	   //mesh = new CMesh();
-	}
+	 }
 
    CDomain::~CDomain(void)
    {
-	   //delete mesh;
    }
 
    ///---------------------------------------------------------------
 
-   void CDomain::assignMesh(const StdString meshName)
-   {
-     mesh = CMesh::getMesh(meshName);
+   void CDomain::assignMesh(const StdString meshName, const int nvertex)
+   {
+     mesh = CMesh::getMesh(meshName, nvertex);
    }
 
@@ -1897,4 +1894,5 @@
 
     buffer >> lon;
+
     if (hasBounds) buffer >> boundslon;
 
Index: /XIOS/trunk/src/node/domain.hpp
===================================================================
--- /XIOS/trunk/src/node/domain.hpp	(revision 923)
+++ /XIOS/trunk/src/node/domain.hpp	(revision 924)
@@ -67,5 +67,5 @@
          
          CMesh* mesh;
-         void assignMesh(const StdString);
+         void assignMesh(const StdString, const int);
         
          virtual void parse(xml::CXMLNode & node);
Index: /XIOS/trunk/src/node/mesh.cpp
===================================================================
--- /XIOS/trunk/src/node/mesh.cpp	(revision 923)
+++ /XIOS/trunk/src/node/mesh.cpp	(revision 924)
@@ -11,6 +11,8 @@
 /// ////////////////////// DÃ©finitions ////////////////////// ///
 
-  CMesh::CMesh(void) :  nbFaces{0}, nbNodes{0}, nbEdges{0}
-            ,   nvertex{0}
+  CMesh::CMesh(void) :  nbNodesGlo{0}, nbEdgesGlo{0}
+            ,  node_start{0}, node_count{0}
+            ,  edge_start{0}, edge_count{0}
+            ,  nbFaces{0}, nbNodes{0}, nbEdges{0}
             ,  nodesAreWritten{false}, edgesAreWritten{false}, facesAreWritten{false}
             ,  node_lon(), node_lat()
@@ -18,4 +20,5 @@
             ,  face_lon(), face_lat()
             ,  face_nodes()
+            ,  pNodeGlobalIndex{NULL}, pEdgeGlobalIndex{NULL}
   {
   }
@@ -24,8 +27,10 @@
   CMesh::~CMesh(void)
   {
+    if (pNodeGlobalIndex != NULL) delete pNodeGlobalIndex;
+    if (pEdgeGlobalIndex != NULL) delete pEdgeGlobalIndex;
   }
 
   std::map <StdString, CMesh> CMesh::meshList = std::map <StdString, CMesh>();
-  //CMesh* CMesh::getMesh;
+  std::map <StdString, vector<int> > CMesh::domainList = std::map <StdString, vector<int> >();
 
 ///---------------------------------------------------------------
@@ -34,7 +39,10 @@
  * Returns a pointer to a mesh. If a mesh has not been created, creates it and adds its name to the list of meshes meshList.
  * \param [in] meshName  The name of a mesh ("name" attribute of a domain).
+ * \param [in] nvertex Number of verteces (1 for nodes, 2 for edges, 3 and up for faces).
  */
-  CMesh* CMesh::getMesh (StdString meshName)
+  CMesh* CMesh::getMesh (StdString meshName, int nvertex)
   {
+    CMesh::domainList[meshName].push_back(nvertex);
+
     if ( CMesh::meshList.begin() != CMesh::meshList.end() )
     {
@@ -60,8 +68,75 @@
 
 ///----------------------------------------------------------------
-  int hashPair(int first, int second)
+  size_t hashPair(size_t first, size_t second)
   {
-    int seed = first + 0x9e3779b9 ;
-    seed ^= second + 0x9e3779b9 + (seed << 6) + (seed >> 2);
+    HashXIOS<size_t> sizetHash;
+    size_t seed = sizetHash(first) + 0x9e3779b9 ;
+    seed ^= sizetHash(second) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
+    return seed ;
+  }
+
+///----------------------------------------------------------------
+  size_t hashPairOrdered(size_t first, size_t second)
+  {
+    size_t seed;
+    HashXIOS<size_t> sizetHash;
+    if (first < second)
+    {
+      seed = sizetHash(first) + 0x9e3779b9 ;
+      seed ^= sizetHash(second) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
+    }
+    else
+    {
+      seed = sizetHash(second) + 0x9e3779b9 ;
+      seed ^= sizetHash(first) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
+    }
+    return seed ;
+  }
+
+///----------------------------------------------------------------
+/*!
+ * \fn size_t generateEdgeIndex(size_t first, size_t second, int rank)
+ * Generates an edge index.
+ * If the same edge is generated by two processes, each process will have its own edge index.
+ * \param [in] first Edge node.
+ * \param [in] second Edge node.
+ * \param [in] rank MPI process rank.
+ */
+  size_t generateEdgeIndex(size_t first, size_t second, int rank)
+  {
+    size_t seed = rank ;
+    if (first < second)
+    {
+      seed = hashPair(seed, first);
+      seed = hashPair(seed, second);
+    }
+    else
+    {
+      seed = hashPair(seed, second);
+      seed = hashPair(seed, first);
+    }
+
+    return seed ;
+  }
+
+///----------------------------------------------------------------
+/*!
+ * \fn size_t generateNodeIndex(vector<size_t>& valList, int rank)
+ * Generates a node index.
+ * If the same node is generated by two processes, each process will have its own node index.
+ * \param [in] valList Vector storing four node hashes.
+ * \param [in] rank MPI process rank.
+ */
+  size_t generateNodeIndex(vector<size_t>& valList, int rank)
+  {
+    // Sort is needed to avoid problems for nodes with lon = 0 generated by faces in east and west semisphere
+    vector<size_t> vec = valList;
+    sort (vec.begin(), vec.end());
+    size_t seed = rank ;
+    int it = 0;
+    for(; it != vec.size(); ++it)
+    {
+       seed = hashPair(seed, vec[it]);
+    }
     return seed ;
   }
@@ -144,19 +219,23 @@
 ///----------------------------------------------------------------
 /*!
- * \fn CArray<size_t,1>& CMesh::createHashes (double lon, double lat)
+ * \fn CArray<size_t,1>& CMesh::createHashes (const double longitude, const double latitude)
  * Creates two hash values for each dimension, longitude and latitude.
- * \param [in] lon Node longitude in degrees.
- * \param [in] lat Node latitude in degrees ranged from 0 to 360.
+ * \param [in] longitude Node longitude in degrees.
+ * \param [in] latitude Node latitude in degrees ranged from 0 to 360.
  */
 
-  vector<size_t> CMesh::createHashes (double lon, double lat)
+  vector<size_t> CMesh::createHashes (const double longitude, const double latitude)
   {
     double minBoundLon = 0. ;
     double maxBoundLon = 360. ;
-    double minBoundLat = -90 ;
-    double maxBoundLat = 90 ;
+    double minBoundLat = -90. ;
+    double maxBoundLat = 90. ;
     double prec=1e-11 ;
     double precLon=prec ;
     double precLat=prec ;
+    double lon = longitude;
+    double lat = latitude;
+
+    if (lon > (360.- prec)) lon = 0.;
 
     size_t maxsize_t=numeric_limits<size_t>::max() ;
@@ -213,14 +292,4 @@
       return std::pair<int,int>(b,a);
   }
-
-///----------------------------------------------------------------
-//#include <random>
-//  size_t randNum()
-//  {
-//    typedef mt19937 rng;
-//    size_t max = numeric_limits<size_t>::max();
-//    uniform_int_distribution<> distr(0, max);
-//    return distr(rng);
-//  }
 
 ///----------------------------------------------------------------
@@ -237,5 +306,5 @@
             const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
   {
-    nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows();
+    int nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows();
 
     if (nvertex == 1)
@@ -357,5 +426,5 @@
             edge_lon(nbEdges) = ( abs( node_lon(face_nodes(nv1,nf)) - node_lon(face_nodes(nv2,nf))) < 180.) ?
                         (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5) :
-                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5 -180.);;
+                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5 -180.);
             edge_lat(nbEdges) = ( node_lat(face_nodes(nv1,nf)) + node_lat(face_nodes(nv2,nf)) ) * 0.5;
             ++nbEdges;
@@ -418,4 +487,5 @@
  * Creates or updates a mesh for the three types of mesh elements: nodes, edges, and faces.
  * Precision check is implemented with two hash values for each dimension, longitude and latitude.
+ * \param [in] comm
  * \param [in] lonvalue  Array of longitudes.
  * \param [in] latvalue  Array of latitudes.
@@ -423,73 +493,17 @@
  * \param [in] bounds_lat Array of boundary latitudes. Its size depends on the element type.
  */
-  void CMesh::createMeshEpsilon(const MPI_Comm& comm, const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
-            const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
+  void CMesh::createMeshEpsilon(const MPI_Comm& comm,
+                                const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
+                                const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
   {
 
-    nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows();
+    int nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows();
+    int mpiRank, mpiSize;
+    MPI_Comm_rank(comm, &mpiRank);
+    MPI_Comm_size(comm, &mpiSize);
 
     if (nvertex == 1)
     {
       nbNodes = lonvalue.numElements();
-
-//      CClientClientDHTSizet::Index2VectorInfoTypeMap hash2Idx;      // map <hash, idx>
-      vector<size_t> hashValues(4);                                 // temporary vector for storing hashes for each node
-      vector<size_t> globalIndexes(nbNodes);                        // Probably a map ?
-      CArray<size_t,1> hashList(nbNodes*4);                         // Do I need it?
-      CArray<size_t,1> idxList(nbNodes);
-
-      // Assign a unique index for each node represented by four hashes
-      // The first hash will serve as the unique index
-      size_t randomIdx;
-      int rank, size;
-      MPI_Comm_rank(comm, &rank);
-      MPI_Comm_size(comm, &size);
-      srand((unsigned)time(NULL)+rank*size);
-
-      for (size_t nn = 0; nn < nbNodes; ++nn)
-      {
-        hashValues = CMesh::createHashes(lonvalue(nn), latvalue(nn));
-        idxList(nn) = hashValues[nn];
-//        randomIdx = rand() % numeric_limits<size_t>::max();
-//        idxList(nn) = randomIdx;
-//        for (size_t nh = 0; nh < 4; ++nh)
-//        {
-//          hash2Idx[hashValues[nh]].push_back(randomIdx);
-//          hashList(nn*4 + nh) = hashValues[nh];
-//        }
-      }
-
-
-//      CClientClientDHTSizet dhtSizet(hash2Idx, comm);
-//      dhtSizet.computeIndexInfoMapping(hashList);
- //     CClientClientDHTSizet::Index2VectorInfoTypeMap& hashIdxList = dhtSizet.getInfoIndexMap();
-
-//      // (2.1) Create a list of indexes for each hush
-//      CClientClientDHTSizet dhtSizet(hash2Idx, comm);
-//      dhtSizet.computeIndexInfoMapping(hashList);
-//      CClientClientDHTSizet::Index2VectorInfoTypeMap& hashIdxList = dhtSizet.getInfoIndexMap();
-//
-//      // (2.2) If more than one index assigned to a hush, set the larger value to be equal to the smaller
-//      int iidx = 0;
-//      for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = hashIdxList.begin(); it != hashIdxList.end(); ++it, ++iidx)
-//      {
-//        vector<size_t> tmp = it->second;
-//        size_t idxMinValue = it->second[0];
-//        for (int i = 1; i < it->second.size(); ++i)
-//        {
-//          if (it->second[i] < idxMinValue )
-//              idxMinValue = it->second[i];
-//        }
-//        if ( (iidx % 4) == 0 )
-//          idxList(iidx/4) = idxMinValue;
-//      }
-
-      // Unique global indexing
-//      CDHTAutoIndexing dhtSizetIdx = CDHTAutoIndexing(idxList, comm);
-//      CClientClientDHTSizet* pDhtSizet = &dhtSizetIdx;
-//      pDhtSizet->computeIndexInfoMapping(idxList);
-
-      //globalIndexes = dhtSizetIdx.getGlobalIndex();
-
       node_lon.resize(nbNodes);
       node_lat.resize(nbNodes);
@@ -497,6 +511,17 @@
       node_lat = latvalue;
 
+      // Global node indexes
+      vector<size_t> hashValues(4);
+      CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2IdxGlo;
+      for (size_t nn = 0; nn < nbNodes; ++nn)
+      {
+        hashValues = CMesh::createHashes(lonvalue(nn), latvalue(nn));
+        for (size_t nh = 0; nh < 4; ++nh)
+        {
+          nodeHash2IdxGlo[hashValues[nh]].push_back(mpiRank*nbNodes + nn);
+        }
+      }
+      pNodeGlobalIndex = new CClientClientDHTSizet (nodeHash2IdxGlo, comm);
       nodesAreWritten = true;
-
     }
 
@@ -504,121 +529,914 @@
     {
       nbEdges = bounds_lon.shape()[1];
-
-      vector<size_t> hashValues(4);
-      for (int ne = 0; ne < nbEdges; ++ne)
+      edge_lon.resize(nbEdges);
+      edge_lat.resize(nbEdges);
+      edge_lon = lonvalue;
+      edge_lat = latvalue;
+      edge_nodes.resize(nvertex, nbEdges);
+      CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2IdxGlo;
+      CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx;
+
+      // Case (1): node indexes have been generated by domain "nodes"
+      if (nodesAreWritten)
       {
-        for (int nv = 0; nv < nvertex; ++nv)
-        {
-          hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
-//          for (size_t nh = 0; nh < 4; ++nh)
-//          {
-//            hash2Idx[hashValues[nh]].push_back(randomIdx);
-//            hashList(nn*4 + nh) = hashValues[nh];
-//          }
-
-
+        vector<size_t> hashValues(4);
+        CArray<size_t,1> nodeHashList(nbEdges*nvertex*4);
+        for (int ne = 0; ne < nbEdges; ++ne)      // size_t doesn't work with CArray<double, 2>
+        {
+          for (int nv = 0; nv < nvertex; ++nv)
+          {
+            hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
+            for (int nh = 0; nh < 4; ++nh)
+            {
+              nodeHashList((ne*nvertex + nv)*4 + nh) = hashValues[nh];
+            }
+          }
+        }
+
+        // Recuperating the node global indexing and writing edge_nodes
+        // Creating map edgeHash2IdxGlo <hash, idxGlo>
+        pNodeGlobalIndex->computeIndexInfoMapping(nodeHashList);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap();
+        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it;
+        vector <size_t> edgeNodes;
+        for (int ne = 0; ne < nbEdges; ++ne)
+        {
+          for (int nv = 0; nv < nvertex; ++nv)
+          {
+            int nh = 0;
+            it = nodeHash2IdxGlo.find(nodeHashList((ne*nvertex + nv)*4 + nh));
+            // The loop below is needed in case if a hash generated by domain "edges" differs
+            // from that generated by domain "nodes" because of possible precision issues
+            while (it == nodeHash2IdxGlo.end())
+            {
+              ++nh;
+              it = nodeHash2IdxGlo.find(nodeHashList((ne*nvertex + nv)*4 + nh));
+            }
+            edge_nodes(nv,ne) = it->second[0];
+            edgeNodes.push_back(it->second[0]);
+          }
+          edgeHash2IdxGlo[ hashPairOrdered(edgeNodes[0], edgeNodes[1]) ].push_back(nbEdges*mpiRank + ne);
+        }
+      } // nodesAreWritten
+
+      // Case (2): node indexes have not been generated previously
+      else
+      {
+        // (2.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> >
+        vector<size_t> hashValues(4);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx;
+        CArray<size_t,1> nodeHashList(nbEdges*nvertex*4);
+        for (int ne = 0; ne < nbEdges; ++ne)
+        {
+          for (int nv = 0; nv < nvertex; ++nv)
+          {
+            hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
+            for (int nh = 0; nh < 4; ++nh)
+            {
+              if (nodeHash2Idx[hashValues[nh]].size() == 0)
+              {
+                nodeHash2Idx[hashValues[nh]].push_back(generateNodeIndex(hashValues, mpiRank));
+                nodeHash2Idx[hashValues[nh]].push_back(mpiRank);
+              }
+              nodeHashList((ne*nvertex + nv)*4 + nh) = hashValues[nh];
+            }
+          }
+        }
+
+        // (2.2) Generating global node indexes
+        // The ownership criterion: priority of the process holding the smaller index
+        // Maps generated in this step are:
+        // nodeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]>
+        // nodeIdx2IdxMin = <idx, idxMin>
+        // nodeIdx2IdxGlo = <idxMin, idxMin>
+
+        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin;
+        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxGlo;
+        CArray<size_t,1> nodeIdxMinList(nbEdges*nvertex);
+
+        CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
+        dhtNodeHash.computeIndexInfoMapping(nodeHashList);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();
+        size_t iIdxMin = 0;
+
+        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
+        {
+          size_t idx = (it->second)[0];
+          size_t idxMin = (it->second)[0];
+          for (int i = 2; i < (it->second).size();)
+          {
+            if (mpiRank == (it->second)[i+1])
+            {
+              idx = (it->second)[i];
+            }
+            if ((it->second)[i] < idxMin)
+            {
+                idxMin = (it->second)[i];
+                (it->second)[i] = (it->second)[i-2];
+            }
+            i += 2;
+          }
+          (it->second)[0] = idxMin;
+          if (nodeIdx2IdxMin.count(idx) == 0)
+          {
+            nodeIdx2IdxMin[idx].push_back(idxMin);
+            if (idx == idxMin)
+              nodeIdx2IdxGlo[idxMin].push_back(idxMin);
+            nodeIdxMinList(iIdxMin) = idxMin;
+            ++iIdxMin;
+          }
+        }
+        nodeIdxMinList.resizeAndPreserve(iIdxMin);
+        CDHTAutoIndexing dhtNodeIdxGlo = CDHTAutoIndexing(nodeIdx2IdxGlo, comm);
+        CClientClientDHTSizet dhtNodeIdx(nodeIdx2IdxGlo, comm);
+        dhtNodeIdx.computeIndexInfoMapping(nodeIdxMinList);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdxMin2IdxGlo = dhtNodeIdx.getInfoIndexMap();
+        // nodeIdx2IdxGlo holds global indexes only for nodes owned by a process
+        // nodeIdxMin2IdxGlo holds global indexes for all nodes generated by a process
+
+        // (2.3) Saving variables: node_lon, node_lat, edge_nodes
+        // Creating map nodeHash2IdxGlo <hash, idxGlo>
+        // Creating map edgeHash2IdxGlo <hash, idxGlo>
+        nbNodesGlo = dhtNodeIdxGlo.getNbIndexesGlobal();
+        node_count = dhtNodeIdxGlo.getIndexCount();
+        node_start = dhtNodeIdxGlo.getIndexStart();
+        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2IdxGlo;
+        node_lon.resize(node_count);
+        node_lat.resize(node_count);
+        vector <size_t> edgeNodes;
+        size_t idxGlo = 0;
+
+        for (int ne = 0; ne < nbEdges; ++ne)
+        {
+          for (int nv = 0; nv < nvertex; ++nv)
+          {
+            hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
+            size_t myIdx = generateNodeIndex(hashValues, mpiRank);
+            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = nodeIdx2IdxMin.find(myIdx);
+            size_t ownerIdx = (itIdx->second)[0];
+
+            if (myIdx == ownerIdx)
+            {
+              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = nodeIdx2IdxGlo.find(myIdx);
+              idxGlo = (itIdxGlo->second)[0];
+//              node_lon(idxGlo - node_start) = (bounds_lon(nv, ne) == 360.) ? (0.) : (bounds_lon(nv, ne));
+              node_lon(idxGlo - node_start) = bounds_lon(nv, ne);
+              node_lat(idxGlo - node_start) = bounds_lat(nv, ne);
+            }
+            else
+            {
+              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = nodeIdxMin2IdxGlo.find(ownerIdx);
+              idxGlo = (itIdxGlo->second)[0];
+            }
+
+            edge_nodes(nv,ne) = idxGlo;
+            for (int nh = 0; nh < 4; ++nh)
+              nodeHash2IdxGlo[hashValues[nh]].push_back(idxGlo);
+            edgeNodes.push_back(idxGlo);
+          }
+          edgeHash2IdxGlo[ hashPairOrdered(edgeNodes[0], edgeNodes[1]) ].push_back(nbEdges*mpiRank + ne);
+          edgeNodes.clear();
+        }
+        pNodeGlobalIndex = new CClientClientDHTSizet (nodeHash2IdxGlo, comm);
+      } // !nodesAreWritten
+
+      pEdgeGlobalIndex = new CClientClientDHTSizet (edgeHash2IdxGlo, comm);
+      edgesAreWritten = true;
+    } //nvertex = 2
+
+    else
+    {
+      nbFaces = bounds_lon.shape()[1];
+      face_lon.resize(nbFaces);
+      face_lat.resize(nbFaces);
+      face_lon = lonvalue;
+      face_lat = latvalue;
+      face_nodes.resize(nvertex, nbFaces);
+      face_edges.resize(nvertex, nbFaces);
+
+      // Case (1): edges have been previously generated
+      if (edgesAreWritten)
+      {
+        // (1.1) Recuperating node global indexing and saving face_nodes
+        vector<size_t> hashValues(4);
+        CArray<size_t,1> nodeHashList(nbFaces*nvertex*4);
+        for (int nf = 0; nf < nbFaces; ++nf)
+        {
+          for (int nv = 0; nv < nvertex; ++nv)
+            {
+              hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
+              for (int nh = 0; nh < 4; ++nh)
+                nodeHashList((nf*nvertex + nv)*4 + nh) = hashValues[nh];
+            }
+        }
+        pNodeGlobalIndex->computeIndexInfoMapping(nodeHashList);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap();
+        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2;
+        CArray<size_t,1> edgeHashList(nbFaces*nvertex);
+        for (int nf = 0; nf < nbFaces; ++nf)
+        {
+          for (int nv1 = 0; nv1 < nvertex; ++nv1)
+          {
+            int nh1 = 0;
+            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
+            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
+            while (it1 == nodeHash2IdxGlo.end())
+            {
+              ++nh1;
+              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
+            }
+            int nh2 = 0;
+            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
+            while (it2 == nodeHash2IdxGlo.end())
+            {
+              ++nh2;
+              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
+            }
+            face_nodes(nv1,nf) = it1->second[0];
+            edgeHashList(nf*nvertex + nv1) = hashPairOrdered(it1->second[0], it2->second[0]);
+
+          }
+        }
+
+        // (1.2) Recuperating edge global indexing and saving face_edges
+        pEdgeGlobalIndex->computeIndexInfoMapping(edgeHashList);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2IdxGlo = pEdgeGlobalIndex->getInfoIndexMap();
+        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itEdgeHash;
+        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Rank;
+
+        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face;
+        CArray<size_t,1> edgeIdxGloList(nbFaces*nvertex);
+        size_t iIdx = 0;
+
+        for (int nf = 0; nf < nbFaces; ++nf)
+        {
+          for (int nv1 = 0; nv1 < nvertex; ++nv1)
+          {
+            int nh1 = 0;
+            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
+            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
+            while (it1 == nodeHash2IdxGlo.end())
+            {
+              ++nh1;
+              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
+            }
+            int nh2 = 0;
+            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
+            while (it2 == nodeHash2IdxGlo.end())
+            {
+              ++nh2;
+              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
+            }
+            size_t faceIdxGlo = mpiRank * nbFaces + nf;
+            size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]);
+            itEdgeHash = edgeHash2IdxGlo.find(edgeHash);
+            size_t edgeIdxGlo = itEdgeHash->second[0];
+            face_edges(nv1,nf) = edgeIdxGlo;
+            if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0)
+            {
+              edgeIdxGloList(iIdx) = edgeIdxGlo;
+              ++iIdx;
+            }
+            edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo);
+            edgeHash2Rank[edgeHash].push_back(itEdgeHash->first);
+            edgeHash2Rank[edgeHash].push_back(mpiRank);
+          }
+        }
+        edgeIdxGloList.resizeAndPreserve(iIdx);
+
+        // (1.3) Saving remaining variables edge_faces and face_faces
+
+        // Establishing edge ownership
+        // The ownership criterion: priority of the process with smaller rank
+        CClientClientDHTSizet dhtEdgeHash (edgeHash2Rank, comm);
+        dhtEdgeHash.computeIndexInfoMapping(edgeHashList);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap();
+        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHashMin2IdxGlo;            // map is needed purely for counting
+
+        // edgeHash2Info = <edgeHash, <idxGlo, rank1, idxGlo, rank2>>
+        // edgeHashMin2IdxGlo edgeHash2Info = <edgeHashMin, idxGlo>
+        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it)
+        {
+          vector <size_t> edgeInfo = it->second;
+          if (edgeInfo.size() == 4)                       // two processes generate the same edge
+            if (edgeInfo[1] > edgeInfo[3])
+              edgeInfo[1] = edgeInfo[3];
+          if (edgeInfo[1] == mpiRank)
+            edgeHashMin2IdxGlo[it->first].push_back(edgeInfo[0]);
+        }
+
+        CDHTAutoIndexing dhtEdgeIdxGlo = CDHTAutoIndexing(edgeHashMin2IdxGlo, comm);
+        edge_count = dhtEdgeIdxGlo.getIndexCount();
+        edge_start = dhtEdgeIdxGlo.getIndexStart();
+
+        edge_faces.resize(2, edge_count);
+        face_faces.resize(nvertex, nbFaces);
+
+        CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm);
+        dhtEdge2Face.computeIndexInfoMapping(edgeIdxGloList);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap();
+        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itFace1, itFace2;
+
+        for (int nf = 0; nf < nbFaces; ++nf)
+        {
+          for (int nv1 = 0; nv1 < nvertex; ++nv1)
+          {
+            int nh1 = 0;
+            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
+            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
+            while (it1 == nodeHash2IdxGlo.end())
+            {
+              ++nh1;
+              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
+            }
+            int nh2 = 0;
+            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
+            while (it2 == nodeHash2IdxGlo.end())
+            {
+              ++nh2;
+              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
+            }
+            size_t faceIdxGlo = mpiRank * nbFaces + nf;
+            size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]);
+            itEdgeHash = edgeHash2Info.find(edgeHash);
+            size_t edgeOwnerRank = itEdgeHash->second[1];
+            int edgeIdxGlo = itEdgeHash->second[0];
+
+            if (mpiRank == edgeOwnerRank)
+            {
+              itFace1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
+              int face1 = itFace1->second[0];
+              if (itFace1->second.size() == 1)
+              {
+                edge_faces(0, edgeIdxGlo - edge_start) = face1;
+                edge_faces(1, edgeIdxGlo - edge_start) = -999;
+                face_faces(nv1, nf) = -1;
+              }
+              else
+              {
+                int face2 = itFace1->second[1];
+                edge_faces(0, edgeIdxGlo - edge_start) = face1;
+                edge_faces(1, edgeIdxGlo - edge_start) = face2;
+                face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
+              }
+            }
+            else
+            {
+              itFace1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
+              int face1 = itFace1->second[0];
+              int face2 = itFace1->second[1];
+              face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
+            }
+          }
+        }
+      } // edgesAreWritten
+
+      // Case (2): nodes have been previously generated
+      else if (nodesAreWritten)
+      {
+        // (2.1) Generating nodeHashList
+        vector<size_t> hashValues(4);
+        CArray<size_t,1> nodeHashList(nbFaces*nvertex*4);
+        for (int nf = 0; nf < nbFaces; ++nf)
+        {
+          for (int nv = 0; nv < nvertex; ++nv)
+            {
+              hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
+              for (int nh = 0; nh < 4; ++nh)
+                nodeHashList((nf*nvertex + nv)*4 + nh) = hashValues[nh];
+            }
+        }
+
+        // (2.2) Recuperating node global indexing and saving face_nodes
+        // Generating edgeHash2Info = <hash, <idx, rank>> and edgeHashList
+        pNodeGlobalIndex->computeIndexInfoMapping(nodeHashList);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap();
+        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx;
+        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2;
+        CArray<size_t,1> edgeHashList(nbFaces*nvertex);
+        for (int nf = 0; nf < nbFaces; ++nf)
+        {
+          for (int nv1 = 0; nv1 < nvertex; ++nv1)
+          {
+            int nh1 = 0;
+            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
+            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
+            while (it1 == nodeHash2IdxGlo.end())
+            {
+              ++nh1;
+              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
+            }
+            int nh2 = 0;
+            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
+            while (it2 == nodeHash2IdxGlo.end())
+            {
+              ++nh2;
+              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
+            }
+            face_nodes(nv1,nf) = it1->second[0];
+            size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]);
+            edgeHash2Idx[edgeHash].push_back(generateEdgeIndex(it1->second[0], it2->second[0], mpiRank));
+            edgeHash2Idx[edgeHash].push_back(mpiRank);
+            edgeHashList(nf*nvertex + nv1) = edgeHash;
+          }
+        }
+
+        // (2.2) Generating global edge indexes
+        // The ownership criterion: priority of the process holding the smaller index
+        // Maps generated in this step are:
+        // edgeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]>
+        // edgeIdx2IdxMin = = <idx, idxMin>
+        // edgeIdx2IdxGlo = <idxMin, idxGlo>
+
+        CClientClientDHTSizet dhtEdgeHash(edgeHash2Idx, comm);
+        dhtEdgeHash.computeIndexInfoMapping(edgeHashList);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap();
+
+        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxMin;
+        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxGlo;
+        CArray<size_t,1> edgeIdxMinList(nbFaces*nvertex);
+        size_t iIdxMin = 0;
+
+        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it)
+        {
+          size_t idxMin = (it->second)[0];
+          size_t idx = (it->second)[0];
+          for (int i = 2; i < (it->second).size();)
+          {
+            if (mpiRank == (it->second)[i+1])
+            {
+              idx = (it->second)[i];
+            }
+            if ((it->second)[i] < idxMin)
+            {
+                idxMin = (it->second)[i];
+                (it->second)[i] = (it->second)[i-2];
+            }
+            i += 2;
+          }
+          (it->second)[0] = idxMin;
+          if (edgeIdx2IdxMin.count(idx) == 0)
+          {
+            edgeIdx2IdxMin[idx].push_back(idxMin);
+            if (idx == idxMin)
+              edgeIdx2IdxGlo[idxMin].push_back(idxMin);
+            edgeIdxMinList(iIdxMin) = idxMin;
+            ++iIdxMin;
+          }
+        }
+        edgeIdxMinList.resizeAndPreserve(iIdxMin);
+        CDHTAutoIndexing dhtEdgeIdxGlo = CDHTAutoIndexing(edgeIdx2IdxGlo, comm);
+        CClientClientDHTSizet dhtEdgeIdx(edgeIdx2IdxGlo, comm);
+        dhtEdgeIdx.computeIndexInfoMapping(edgeIdxMinList);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxMin2IdxGlo = dhtEdgeIdx.getInfoIndexMap();
+        // edgeIdx2IdxGlo holds global indexes only for edges owned by a process
+        // edgeIdxMin2IdxGlo holds global indexes for all edges generated by a process
+
+        // (2.3) Saving variables: edge_lon, edge_lat, face_edges
+        nbEdgesGlo = dhtEdgeIdxGlo.getNbIndexesGlobal();
+        edge_count = dhtEdgeIdxGlo.getIndexCount();
+        edge_start = dhtEdgeIdxGlo.getIndexStart();
+        edge_lon.resize(edge_count);
+        edge_lat.resize(edge_count);
+        edge_nodes.resize(2, edge_count);
+        face_edges.resize(nvertex, nbFaces);
+
+        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face;
+        CArray<size_t,1> edgeIdxGloList(nbFaces*nvertex);
+        size_t iIdx = 0;
+
+        for (int nf = 0; nf < nbFaces; ++nf)
+        {
+          for (int nv1 = 0; nv1 < nvertex; ++nv1)
+          {
+            // Getting global indexes of edge's nodes
+            int nh1 = 0;
+            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
+            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
+            while (it1 == nodeHash2IdxGlo.end())
+            {
+              ++nh1;
+              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
+            }
+            int nh2 = 0;
+            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
+            while (it2 == nodeHash2IdxGlo.end())
+            {
+              ++nh2;
+              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
+            }
+            // Getting edge global index
+            size_t nodeIdxGlo1 = it1->second[0];
+            size_t nodeIdxGlo2 = it2->second[0];
+            size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank);
+            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx);
+            size_t ownerIdx = (itIdx->second)[0];
+            int edgeIdxGlo = 0;
+            size_t faceIdxGlo = mpiRank * nbFaces + nf;
+
+            if (myIdx == ownerIdx)
+            {
+              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx);
+              edgeIdxGlo = (itIdxGlo->second)[0];
+              edge_lon(edgeIdxGlo - edge_start) = ( abs( bounds_lon(nv1, nf) - bounds_lon(nv2, nf)) < 180.) ?
+                          (( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5) :
+                          (( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5 -180.);
+              edge_lat(edgeIdxGlo - edge_start) = ( bounds_lat(nv1, nf) + bounds_lat(nv2, nf) ) * 0.5;
+              edge_nodes(0, edgeIdxGlo - edge_start) = nodeIdxGlo1;
+              edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2;
+            }
+            else
+            {
+              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx);
+              edgeIdxGlo = (itIdxGlo->second)[0];
+            }
+            face_edges(nv1,nf) = edgeIdxGlo;
+            if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0)
+            {
+              edgeIdxGloList(iIdx) = edgeIdxGlo;
+              ++iIdx;
+            }
+            edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo);
+          }
+        }
+        edgeIdxGloList.resizeAndPreserve(iIdx);
+
+        // (2.4) Saving remaining variables edge_faces and face_faces
+        edge_faces.resize(2, edge_count);
+        face_faces.resize(nvertex, nbFaces);
+
+        CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm);
+        dhtEdge2Face.computeIndexInfoMapping(edgeIdxGloList);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap();
+        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdxGlo1, itNodeIdxGlo2;
+
+        for (int nf = 0; nf < nbFaces; ++nf)
+        {
+          for (int nv1 = 0; nv1 < nvertex; ++nv1)
+          {
+            // Getting global indexes of edge's nodes
+            int nh1 = 0;
+            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
+            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
+            while (it1 == nodeHash2IdxGlo.end())
+            {
+              ++nh1;
+              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
+            }
+            int nh2 = 0;
+            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
+            while (it2 == nodeHash2IdxGlo.end())
+            {
+              ++nh2;
+              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
+            }
+            size_t nodeIdxGlo1 = it1->second[0];
+            size_t nodeIdxGlo2 = it2->second[0];
+
+            size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank);
+            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx);
+            size_t ownerIdx = (itIdx->second)[0];
+            size_t faceIdxGlo = mpiRank * nbFaces + nf;
+
+            if (myIdx == ownerIdx)
+            {
+              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx);
+              int edgeIdxGlo = (itIdxGlo->second)[0];
+              it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
+              int face1 = it1->second[0];
+              if (it1->second.size() == 1)
+              {
+                edge_faces(0, edgeIdxGlo - edge_start) = face1;
+                edge_faces(1, edgeIdxGlo - edge_start) = -999;
+                face_faces(nv1, nf) = -1;
+              }
+              else
+              {
+                int face2 = it1->second[1];
+                edge_faces(0, edgeIdxGlo - edge_start) = face1;
+                edge_faces(1, edgeIdxGlo - edge_start) = face2;
+                face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
+              }
+            }
+            else
+            {
+              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx);
+              size_t edgeIdxGlo = (itIdxGlo->second)[0];
+              it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
+              int face1 = it1->second[0];
+              int face2 = it1->second[1];
+              face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
+            }
+          }
+        }
+      } // nodesAreWritten
+
+      // Case (3): Neither nodes nor edges have been previously generated
+      else
+      {
+        // (3.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> >
+        vector<size_t> hashValues(4);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx;
+        CArray<size_t,1> nodeHashList(nbFaces*nvertex*4);
+        size_t iHash = 0;
+        for (int nf = 0; nf < nbFaces; ++nf)
+        {
+          for (int nv = 0; nv < nvertex; ++nv)
+          {
+            hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
+            size_t nodeIndex = generateNodeIndex(hashValues, mpiRank);
+            for (int nh = 0; nh < 4; ++nh)
+            {
+              if (nodeHash2Idx.count(hashValues[nh])==0)
+              {
+                nodeHash2Idx[hashValues[nh]].push_back(nodeIndex);
+                nodeHash2Idx[hashValues[nh]].push_back(mpiRank);
+                nodeHashList(iHash) = hashValues[nh];
+                ++iHash;
+              }
+            }
+          }
+        }
+        nodeHashList.resizeAndPreserve(iHash);
+
+        // (3.2) Generating global node indexes
+        // The ownership criterion: priority of the process holding the smaller index
+        // Maps generated in this step are:
+        // nodeHash2Info = <hash, [[idxMin, rankMin], [idx1, rank1], [idx3, rank3]..]>
+        // nodeIdx2IdxMin = = <idx, idxMin>
+        // nodeIdx2IdxGlo = <idxMin, idxGlo>
+
+        CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
+        dhtNodeHash.computeIndexInfoMapping(nodeHashList);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();
+
+        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin;
+        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxGlo;
+        CArray<size_t,1> nodeIdxMinList(nbFaces*nvertex*4);
+        size_t iIdxMin = 0;
+
+        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
+        {
+          size_t idxMin = (it->second)[0];
+          size_t idx = (it->second)[0];
+          for (int i = 2; i < (it->second).size();)
+          {
+            if (mpiRank == (it->second)[i+1])
+            {
+              idx = (it->second)[i];
+            }
+            if ((it->second)[i] < idxMin)
+            {
+                idxMin = (it->second)[i];
+                (it->second)[i] = (it->second)[i-2];
+            }
+            i += 2;
+          }
+          (it->second)[0] = idxMin;
+          if (nodeIdx2IdxMin.count(idx) == 0)
+          {
+            nodeIdx2IdxMin[idx].push_back(idxMin);
+            if (idx == idxMin)
+              nodeIdx2IdxGlo[idxMin].push_back(idxMin);
+            nodeIdxMinList(iIdxMin) = idxMin;
+            ++iIdxMin;
+          }
+        }
+
+        nodeIdxMinList.resizeAndPreserve(iIdxMin);
+        CDHTAutoIndexing dhtNodeIdxGlo = CDHTAutoIndexing(nodeIdx2IdxGlo, comm);
+        CClientClientDHTSizet dhtNodeIdx(nodeIdx2IdxGlo, comm);
+        dhtNodeIdx.computeIndexInfoMapping(nodeIdxMinList);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdxMin2IdxGlo = dhtNodeIdx.getInfoIndexMap();
+
+        // (3.3) Saving node data: node_lon, node_lat, and face_nodes
+        // Generating edgeHash2Info = <hash, <idx, rank>> and edgeHashList
+        nbNodesGlo = dhtNodeIdxGlo.getNbIndexesGlobal();
+        node_count = dhtNodeIdxGlo.getIndexCount();
+        node_start = dhtNodeIdxGlo.getIndexStart();
+        node_lon.resize(node_count);
+        node_lat.resize(node_count);
+        size_t nodeIdxGlo1 = 0;
+        size_t nodeIdxGlo2 = 0;
+        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx;
+        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdxGlo1, itNodeIdxGlo2;
+        CArray<size_t,1> edgeHashList(nbFaces*nvertex);
+
+        for (int nf = 0; nf < nbFaces; ++nf)
+        {
+          for (int nv1 = 0; nv1 < nvertex; ++nv1)
+          {
+            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
+            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
+            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
+            size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank);
+            size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank);
+            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx1 = nodeIdx2IdxMin.find(myNodeIdx1);
+            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx2 = nodeIdx2IdxMin.find(myNodeIdx2);
+            size_t ownerNodeIdx = (itNodeIdx1->second)[0];
+
+            if (myNodeIdx1 == ownerNodeIdx)
+            {
+              itNodeIdxGlo1 = nodeIdx2IdxGlo.find(myNodeIdx1);
+              nodeIdxGlo1 = (itNodeIdxGlo1->second)[0];
+              node_lon(nodeIdxGlo1 - node_start) = bounds_lon(nv1, nf);
+              node_lat(nodeIdxGlo1 - node_start) = bounds_lat(nv1, nf);
+            }
+            else
+            {
+              itNodeIdxGlo1 = nodeIdxMin2IdxGlo.find(ownerNodeIdx);
+              nodeIdxGlo1 = (itNodeIdxGlo1->second)[0];
+            }
+            itNodeIdxGlo2 = nodeIdxMin2IdxGlo.find((itNodeIdx2->second)[0]);
+            nodeIdxGlo2 = (itNodeIdxGlo2->second)[0];
+            size_t edgeHash = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
+            edgeHash2Idx[edgeHash].push_back(generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank));
+            edgeHash2Idx[edgeHash].push_back(mpiRank);
+            edgeHashList(nf*nvertex + nv1) = edgeHash;
+            face_nodes(nv1,nf) = nodeIdxGlo1;
+          }
+        }
+
+        // (3.4) Generating global edge indexes
+        // Maps generated in this step are:
+        // edgeIdx2IdxMin = = <idx, idxMin>
+        // edgeIdx2IdxGlo = <idxMin, idxGlo>
+
+        CClientClientDHTSizet dhtEdgeHash(edgeHash2Idx, comm);
+        dhtEdgeHash.computeIndexInfoMapping(edgeHashList);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap();
+        // edgeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]>
+
+        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxMin;
+        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxGlo;
+        CArray<size_t,1> edgeIdxMinList(nbFaces*nvertex);
+        iIdxMin = 0;
+
+        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it)
+        {
+          size_t idxMin = (it->second)[0];
+          size_t idx = (it->second)[0];
+
+          for (int i = 2; i < (it->second).size();)
+          {
+            if (mpiRank == (it->second)[i+1])
+            {
+              idx = (it->second)[i];
+            }
+            if ((it->second)[i] < idxMin)
+            {
+                idxMin = (it->second)[i];
+                (it->second)[i] = (it->second)[i-2];
+            }
+            i += 2;
+          }
+          (it->second)[0] = idxMin;
+          if (edgeIdx2IdxMin.count(idx) == 0)
+          {
+            edgeIdx2IdxMin[idx].push_back(idxMin);
+            if (idx == idxMin)
+              edgeIdx2IdxGlo[idxMin].push_back(idxMin);
+            edgeIdxMinList(iIdxMin) = idxMin;
+            ++iIdxMin;
+          }
+        }
+        edgeIdxMinList.resizeAndPreserve(iIdxMin);
+        CDHTAutoIndexing dhtEdgeIdxGlo = CDHTAutoIndexing(edgeIdx2IdxGlo, comm);
+        CClientClientDHTSizet dhtEdgeIdx(edgeIdx2IdxGlo, comm);
+        dhtEdgeIdx.computeIndexInfoMapping(edgeIdxMinList);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxMin2IdxGlo = dhtEdgeIdx.getInfoIndexMap();
+
+        // (3.5) Saving variables: edge_lon, edge_lat, face_edges
+        // Creating map edgeIdxGlo2Face <idxGlo, face>
+
+        nbEdgesGlo = dhtEdgeIdxGlo.getNbIndexesGlobal();
+        edge_count = dhtEdgeIdxGlo.getIndexCount();
+        edge_start = dhtEdgeIdxGlo.getIndexStart();
+        edge_lon.resize(edge_count);
+        edge_lat.resize(edge_count);
+        edge_nodes.resize(2, edge_count);
+        face_edges.resize(nvertex, nbFaces);
+
+        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2;
+        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face;
+        CArray<size_t,1> edgeIdxGloList(nbFaces*nvertex);
+        size_t iIdx = 0;
+
+        for (int nf = 0; nf < nbFaces; ++nf)
+        {
+          for (int nv1 = 0; nv1 < nvertex; ++nv1)
+          {
+            // Getting global indexes of edge's nodes
+            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
+            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
+            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
+
+            size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank);
+            size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank);
+            it1 = nodeIdx2IdxMin.find(myNodeIdx1);
+            it2 = nodeIdx2IdxMin.find(myNodeIdx2);
+            itNodeIdxGlo1 = nodeIdxMin2IdxGlo.find((it1->second)[0]);
+            itNodeIdxGlo2 = nodeIdxMin2IdxGlo.find((it2->second)[0]);
+            size_t nodeIdxGlo1 = (itNodeIdxGlo1->second)[0];
+            size_t nodeIdxGlo2 = (itNodeIdxGlo2->second)[0];
+
+            size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank);
+            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx);
+            size_t ownerIdx = (itIdx->second)[0];
+
+            int edgeIdxGlo = 0;
+            size_t faceIdxGlo = mpiRank * nbFaces + nf;
+
+            if (myIdx == ownerIdx)
+            {
+              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx);
+              edgeIdxGlo = (itIdxGlo->second)[0];
+              edge_lon(edgeIdxGlo - edge_start) = ( abs( bounds_lon(nv1, nf) - bounds_lon(nv2, nf)) < 180.) ?
+                          (( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5) :
+                          (( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5 -180.);
+              edge_lat(edgeIdxGlo-edge_start) = ( bounds_lat(nv1, nf) + bounds_lat(nv2, nf) ) * 0.5;
+              edge_nodes(0, edgeIdxGlo - edge_start) = nodeIdxGlo1;
+              edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2;
+            }
+            else
+            {
+              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx);
+              edgeIdxGlo = (itIdxGlo->second)[0];
+            }
+            face_edges(nv1,nf) = edgeIdxGlo;
+            if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0)
+            {
+              edgeIdxGloList(iIdx) = edgeIdxGlo;
+              ++iIdx;
+            }
+            edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo);
+          }
+        }
+        edgeIdxGloList.resizeAndPreserve(iIdx);
+
+        // (3.6) Saving remaining variables edge_faces and face_faces
+        edge_faces.resize(2, edge_count);
+        face_faces.resize(nvertex, nbFaces);
+
+        CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm);
+        dhtEdge2Face.computeIndexInfoMapping(edgeIdxGloList);
+        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap();
+
+        for (int nf = 0; nf < nbFaces; ++nf)
+        {
+          for (int nv1 = 0; nv1 < nvertex; ++nv1)
+          {
+            // Getting global indexes of edge's nodes
+            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
+            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
+            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
+
+            size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank);
+            size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank);
+            it1 = nodeIdx2IdxMin.find(myNodeIdx1);
+            it2 = nodeIdx2IdxMin.find(myNodeIdx2);
+            itNodeIdxGlo1 = nodeIdxMin2IdxGlo.find((it1->second)[0]);
+            itNodeIdxGlo2 = nodeIdxMin2IdxGlo.find((it2->second)[0]);
+            size_t nodeIdxGlo1 = (itNodeIdxGlo1->second)[0];
+            size_t nodeIdxGlo2 = (itNodeIdxGlo2->second)[0];
+
+            size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank);
+            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx);
+            size_t ownerIdx = (itIdx->second)[0];
+            size_t faceIdxGlo = mpiRank * nbFaces + nf;
+
+            if (myIdx == ownerIdx)
+            {
+              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx);
+              int edgeIdxGlo = (itIdxGlo->second)[0];
+              it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
+              int face1 = it1->second[0];
+              if (it1->second.size() == 1)
+              {
+                edge_faces(0, edgeIdxGlo - edge_start) = face1;
+                edge_faces(1, edgeIdxGlo - edge_start) = -999;
+                face_faces(nv1, nf) = -1;
+              }
+              else
+              {
+                size_t face2 = it1->second[1];
+                edge_faces(0, edgeIdxGlo - edge_start) = face1;
+                edge_faces(1, edgeIdxGlo - edge_start) = face2;
+                face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
+              }
+            }
+            else
+            {
+              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx);
+              size_t edgeIdxGlo = (itIdxGlo->second)[0];
+              it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
+              int face1 = it1->second[0];
+              int face2 = it1->second[1];
+              face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
+            }
+          }
         }
       }
-      // Create nodes and edge_node connectivity
-
-      node_lon.resizeAndPreserve(nbEdges*nvertex); // Max possible number of nodes
-      node_lat.resizeAndPreserve(nbEdges*nvertex);
-      edge_nodes.resizeAndPreserve(nvertex, nbEdges);
-
-//      for (int ne = 0; ne < nbEdges; ++ne)
-//      {
-//        for (int nv = 0; nv < nvertex; ++nv)
-//        {
-//          if ( nodeIndex(bounds_lon(nv, ne), bounds_lat(nv ,ne)) == -1)
-//          {
-//            edge_nodes(nv,ne) = nbNodes ;
-//            node_lon(nbNodes) = bounds_lon(nv, ne);
-//            node_lat(nbNodes) = bounds_lat(nv, ne);
-//            ++nbNodes ;
-//          }
-//          else
-//            edge_nodes(nv,ne) = nodeIndex(bounds_lon(nv, ne), bounds_lat(nv ,ne));
-//        }
-//      }
-//      node_lon.resizeAndPreserve(nbNodes);
-//      node_lat.resizeAndPreserve(nbNodes);
-
-      // Create edges
-      edge_lon.resizeAndPreserve(nbEdges);
-      edge_lat.resizeAndPreserve(nbEdges);
-
-      for (int ne = 0; ne < nbEdges; ++ne)
-      {
-        if (map_edges.find(make_ordered_pair (edge_nodes(0,ne), edge_nodes(1,ne))) == map_edges.end())
-        {
-          map_edges[make_ordered_pair ( edge_nodes(0,ne), edge_nodes(1,ne) )] = ne ;
-          edge_lon(ne) = lonvalue(ne);
-          edge_lat(ne) = latvalue(ne);
-        }
-      }
-      edgesAreWritten = true;
-    } // nvertex = 2
-    else
-    {
-      nbFaces = bounds_lon.shape()[1];
-
-      // Create nodes and face_node connectivity
-      node_lon.resizeAndPreserve(nbFaces*nvertex);  // Max possible number of nodes
-      node_lat.resizeAndPreserve(nbFaces*nvertex);
-      face_nodes.resize(nvertex, nbFaces);
-
-      /*for (int nf = 0; nf < nbFaces; ++nf)
-      {
-        for (int nv = 0; nv < nvertex; ++nv)
-        {
-          if ( nodeIndex(bounds_lon(nv, nf), bounds_lat(nv ,nf)) == -1)
-          {
-            face_nodes(nv,nf) = nbNodes ;
-            node_lon(nbNodes) = bounds_lon(nv, nf);
-            node_lat(nbNodes) = bounds_lat(nv ,nf);
-            ++nbNodes ;
-          }
-          else
-          {
-            face_nodes(nv,nf) = nodeIndex(bounds_lon(nv, nf), bounds_lat(nv ,nf));
-          }
-        }
-      }
-      node_lon.resizeAndPreserve(nbNodes);
-      node_lat.resizeAndPreserve(nbNodes);*/
-
-
-      // Create edges and edge_nodes connectivity
-//      edge_lon.resizeAndPreserve(nbFaces*nvertex); // Max possible number of edges
-//      edge_lat.resizeAndPreserve(nbFaces*nvertex);
-//      edge_nodes.resizeAndPreserve(2, nbFaces*nvertex);
-//      for (int nf = 0; nf < nbFaces; ++nf)
-//      {
-//        for (int nv1 = 0; nv1 < nvertex; ++nv1)
-//        {
-//          int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
-//          if (map_edges.find(make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))) == map_edges.end())
-//          {
-//            map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))] = nbEdges ;
-//            edge_nodes(Range::all(),nbEdges) = face_nodes(nv1,nf), face_nodes(nv2,nf);
-//            edge_lon(nbEdges) = ( abs( node_lon(face_nodes(nv1,nf)) - node_lon(face_nodes(nv2,nf))) < 180.) ?
-//                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5) :
-//                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5 -180.);;
-//            edge_lat(nbEdges) = ( node_lat(face_nodes(nv1,nf)) + node_lat(face_nodes(nv2,nf)) ) * 0.5;
-//            ++nbEdges;
-//          }
-//        }
-//      }
-//      edge_nodes.resizeAndPreserve(2, nbEdges);
-//      edge_lon.resizeAndPreserve(nbEdges);
-//      edge_lat.resizeAndPreserve(nbEdges);
-
-      // Create faces
-//      face_lon.resize(nbFaces);
-//      face_lat.resize(nbFaces);
-//      face_lon = lonvalue;
-//      face_lat = latvalue;
-//      facesAreWritten = true;
-
+      facesAreWritten = true;
     } // nvertex >= 3
 
@@ -626,130 +1444,2 @@
 
 } // namespace xios
-
-//  void CMesh::createMeshEpsilon(const MPI_Comm& comm, const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
-//            const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
-//  {
-//
-//    nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows();
-//
-//    if (nvertex == 1)
-//    {
-//      nbNodes = lonvalue.numElements();
-//      node_lon.resizeAndPreserve(nbNodes);
-//      node_lat.resizeAndPreserve(nbNodes);
-//      for (int nn = 0; nn < nbNodes; ++nn)
-//      {
-//        if ( nodeIndex(lonvalue(nn), latvalue(nn)) == -1 )
-//        {
-//          node_lon(nn) = lonvalue(nn);
-//          node_lat(nn) = latvalue(nn);
-//        }
-//      }
-//
-//      nodesAreWritten = true;
-//    }
-//    else if (nvertex == 2)
-//    {
-//      nbEdges = bounds_lon.shape()[1];
-//
-//      // Create nodes and edge_node connectivity
-//      node_lon.resizeAndPreserve(nbEdges*nvertex); // Max possible number of nodes
-//      node_lat.resizeAndPreserve(nbEdges*nvertex);
-//      edge_nodes.resizeAndPreserve(nvertex, nbEdges);
-//
-//      for (int ne = 0; ne < nbEdges; ++ne)
-//      {
-//        for (int nv = 0; nv < nvertex; ++nv)
-//        {
-//          if ( nodeIndex(bounds_lon(nv, ne), bounds_lat(nv ,ne)) == -1)
-//          {
-//            edge_nodes(nv,ne) = nbNodes ;
-//            node_lon(nbNodes) = bounds_lon(nv, ne);
-//            node_lat(nbNodes) = bounds_lat(nv, ne);
-//            ++nbNodes ;
-//          }
-//          else
-//            edge_nodes(nv,ne) = nodeIndex(bounds_lon(nv, ne), bounds_lat(nv ,ne));
-//        }
-//      }
-//      node_lon.resizeAndPreserve(nbNodes);
-//      node_lat.resizeAndPreserve(nbNodes);
-//
-//      // Create edges
-//      edge_lon.resizeAndPreserve(nbEdges);
-//      edge_lat.resizeAndPreserve(nbEdges);
-//
-//      for (int ne = 0; ne < nbEdges; ++ne)
-//      {
-//        if (map_edges.find(make_ordered_pair (edge_nodes(0,ne), edge_nodes(1,ne))) == map_edges.end())
-//        {
-//          map_edges[make_ordered_pair ( edge_nodes(0,ne), edge_nodes(1,ne) )] = ne ;
-//          edge_lon(ne) = lonvalue(ne);
-//          edge_lat(ne) = latvalue(ne);
-//        }
-//      }
-//      edgesAreWritten = true;
-//    } // nvertex = 2
-//    else
-//    {
-//      nbFaces = bounds_lon.shape()[1];
-//
-//      // Create nodes and face_node connectivity
-//      node_lon.resizeAndPreserve(nbFaces*nvertex);  // Max possible number of nodes
-//      node_lat.resizeAndPreserve(nbFaces*nvertex);
-//      face_nodes.resize(nvertex, nbFaces);
-//
-//      for (int nf = 0; nf < nbFaces; ++nf)
-//      {
-//        for (int nv = 0; nv < nvertex; ++nv)
-//        {
-//          if ( nodeIndex(bounds_lon(nv, nf), bounds_lat(nv ,nf)) == -1)
-//          {
-//            face_nodes(nv,nf) = nbNodes ;
-//            node_lon(nbNodes) = bounds_lon(nv, nf);
-//            node_lat(nbNodes) = bounds_lat(nv ,nf);
-//            ++nbNodes ;
-//          }
-//          else
-//          {
-//            face_nodes(nv,nf) = nodeIndex(bounds_lon(nv, nf), bounds_lat(nv ,nf));
-//          }
-//        }
-//      }
-//      node_lon.resizeAndPreserve(nbNodes);
-//      node_lat.resizeAndPreserve(nbNodes);
-//
-//      // Create edges and edge_nodes connectivity
-//      edge_lon.resizeAndPreserve(nbFaces*nvertex); // Max possible number of edges
-//      edge_lat.resizeAndPreserve(nbFaces*nvertex);
-//      edge_nodes.resizeAndPreserve(2, nbFaces*nvertex);
-//      for (int nf = 0; nf < nbFaces; ++nf)
-//      {
-//        for (int nv1 = 0; nv1 < nvertex; ++nv1)
-//        {
-//          int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
-//          if (map_edges.find(make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))) == map_edges.end())
-//          {
-//            map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))] = nbEdges ;
-//            edge_nodes(Range::all(),nbEdges) = face_nodes(nv1,nf), face_nodes(nv2,nf);
-//            edge_lon(nbEdges) = ( abs( node_lon(face_nodes(nv1,nf)) - node_lon(face_nodes(nv2,nf))) < 180.) ?
-//                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5) :
-//                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5 -180.);;
-//            edge_lat(nbEdges) = ( node_lat(face_nodes(nv1,nf)) + node_lat(face_nodes(nv2,nf)) ) * 0.5;
-//            ++nbEdges;
-//          }
-//        }
-//      }
-//      edge_nodes.resizeAndPreserve(2, nbEdges);
-//      edge_lon.resizeAndPreserve(nbEdges);
-//      edge_lat.resizeAndPreserve(nbEdges);
-//
-//      // Create faces
-//      face_lon.resize(nbFaces);
-//      face_lat.resize(nbFaces);
-//      face_lon = lonvalue;
-//      face_lat = latvalue;
-//      facesAreWritten = true;
-//    } // nvertex >= 3
-//
-//  } // createMeshEpsilon
Index: /XIOS/trunk/src/node/mesh.hpp
===================================================================
--- /XIOS/trunk/src/node/mesh.hpp	(revision 923)
+++ /XIOS/trunk/src/node/mesh.hpp	(revision 924)
@@ -9,5 +9,4 @@
  
 #include "array_new.hpp"
-#include "client_client_dht_template_impl.hpp"
 #include "dht_auto_indexing.hpp"
 
@@ -32,8 +31,11 @@
       ~CMesh(void);
     
-      int nbNodes;
-      int nbEdges;
-      int nbFaces;
-      int nvertex;  
+      int nbNodesGlo;
+      int nbEdgesGlo;
+
+      int node_start;
+      int node_count;
+      int edge_start;
+      int edge_count;
 
       bool nodesAreWritten;
@@ -56,15 +58,24 @@
 
       void createMesh(const CArray<double, 1>&, const CArray<double, 1>&, 
-            const CArray<double, 2>&, const CArray<double, 2>& );
+                      const CArray<double, 2>&, const CArray<double, 2>& );
                         
-      void createMeshEpsilon(const MPI_Comm&, const CArray<double, 1>&, const CArray<double, 1>&,
-            const CArray<double, 2>&, const CArray<double, 2>& );
+      void createMeshEpsilon(const MPI_Comm&,
+                             const CArray<double, 1>&, const CArray<double, 1>&,
+                             const CArray<double, 2>&, const CArray<double, 2>& );
             
-      static CMesh* getMesh(StdString);
+      static CMesh* getMesh(StdString, int);
 
     private:
 
+      int nbNodes;
+      int nbEdges;
+      int nbFaces;
+
       static std::map <StdString, CMesh> meshList;
-      vector<size_t> createHashes (double, double);
+      static std::map <StdString, vector<int> > domainList;
+      CClientClientDHTSizet* pNodeGlobalIndex;                    // pointer to a map <nodeHash, nodeIdxGlo>
+      CClientClientDHTSizet* pEdgeGlobalIndex;                    // pointer to a map <edgeHash, edgeIdxGlo>
+
+      vector<size_t> createHashes (const double, const double);
 
       size_t nodeIndex (double, double);                           // redundant in parallel version with epsilon precision
Index: /XIOS/trunk/src/test/test_regular.f90
===================================================================
--- /XIOS/trunk/src/test/test_regular.f90	(revision 924)
+++ /XIOS/trunk/src/test/test_regular.f90	(revision 924)
@@ -0,0 +1,195 @@
+PROGRAM test_regular
+
+  USE xios
+  IMPLICIT NONE
+  INCLUDE "mpif.h"
+
+  INTEGER, PARAMETER :: ntime=1
+
+  CHARACTER(len=*),PARAMETER :: id="client"
+  TYPE(xios_duration)  :: dtime
+  TYPE(xios_context) :: ctx_hdl
+  TYPE(xios_field) :: field_hdl
+  TYPE(xios_fieldgroup) :: fieldgroup_hdl
+  TYPE(xios_file) :: file_hdl
+  LOGICAL :: ok
+  CHARACTER(len=256) :: crname
+  INTEGER :: ts
+  INTEGER :: comm  
+  INTEGER :: ierr
+  INTEGER :: size, rank
+
+  INTEGER :: nlon = 100 
+  INTEGER :: nlat = 100
+  INTEGER :: ncell 
+  INTEGER :: ilat, ilon, ind
+  DOUBLE PRECISION :: lon1, lon2, lat1, lat2
+  INTEGER :: ni, ibegin
+
+  DOUBLE PRECISION,ALLOCATABLE :: lon_glo(:)
+  DOUBLE PRECISION,ALLOCATABLE :: lat_glo(:)
+  DOUBLE PRECISION,ALLOCATABLE :: bounds_lon_glo(:,:)
+  DOUBLE PRECISION,ALLOCATABLE :: bounds_lat_glo(:,:)
+  DOUBLE PRECISION,ALLOCATABLE :: field_temp_glo(:,:)
+  DOUBLE PRECISION,ALLOCATABLE :: lon(:)
+  DOUBLE PRECISION,ALLOCATABLE :: lat(:)
+  DOUBLE PRECISION,ALLOCATABLE :: bounds_lon(:,:)
+  DOUBLE PRECISION,ALLOCATABLE :: bounds_lat(:,:)
+  DOUBLE PRECISION,ALLOCATABLE :: field_temp(:,:)
+
+!!! MPI Initialization
+  CALL MPI_INIT(ierr)
+
+!!! XIOS Initialization (get the local communicator)
+  CALL xios_initialize(id,return_comm=comm)
+
+!###########################################################################
+! Contexte ATM
+!###########################################################################
+
+!!! Initialisation des coordonnÃ©es globales et locales pour la grille rÃ©guliÃšre
+
+! Regions around the poles are not included into the grid
+! The whole grid is rectangular (nvertex=4)
+
+  ncell = nlon * (nlat-1)
+  ALLOCATE(lon_glo(ncell))
+  ALLOCATE(lat_glo(ncell))
+  ALLOCATE(bounds_lon_glo(4,ncell))
+  ALLOCATE(bounds_lat_glo(4,ncell))
+  ALLOCATE(field_temp_glo(ncell,ntime))
+
+
+  ind = 0
+  DO ilat = 1, nlat-1
+   DO ilon = 1, nlon
+
+      ind=ind+1
+
+      lon1 = 360./DBLE(nlon) * DBLE(ilon-1)
+      lon2 = lon1 + 360./DBLE(nlon)
+
+      lat1 = (90. + 90./DBLE(nlat)) - 180./DBLE(nlat)*DBLE(ilat)
+      lat2 = lat1 - 180./DBLE(nlat)
+
+      lon_glo(ind) = (lon1+lon2)*0.5
+      lat_glo(ind) = (lat1+lat2)*0.5 
+
+      bounds_lon_glo(1,ind) = lon1
+      bounds_lon_glo(2,ind) = lon2
+      bounds_lon_glo(3,ind) = lon2
+      bounds_lon_glo(4,ind) = lon1
+
+      bounds_lat_glo(1,ind) = lat1
+      bounds_lat_glo(2,ind) = lat1
+      bounds_lat_glo(3,ind) = lat2      
+      bounds_lat_glo(4,ind) = lat2     
+
+      field_temp_glo(ind,1) = DBLE(ind)
+
+    ENDDO
+  ENDDO
+
+! Initialization of local variables
+
+  CALL MPI_COMM_RANK(comm,rank,ierr)
+  CALL MPI_COMM_SIZE(comm,size,ierr)
+
+  IF (MOD(ncell, size) == 0) THEN
+    ni = ncell/size
+    ibegin = rank*ni
+  ELSE
+    IF (rank < MOD(ncell, size)) THEN
+      ni = ncell/size + 1
+      ibegin = rank*(ncell/size + 1)
+    ELSE
+      ni = ncell/size
+      IF (rank == MOD(ncell, size)) THEN
+        ibegin = rank*(ncell/size + 1)
+      ELSE
+        ibegin = MOD(ncell,size)*(ncell/size + 1) + (rank-MOD(ncell,size))*ncell/size
+      END IF
+    END IF
+  END IF
+
+  ALLOCATE(lon(ni))
+  ALLOCATE(lat(ni))
+  ALLOCATE(bounds_lon(4,ni))
+  ALLOCATE(bounds_lat(4,ni))
+  ALLOCATE(field_temp(ni,ntime)) 
+  lon = lon_glo(1+ibegin:1+ibegin+ni)
+  lat = lat_glo(1+ibegin:1+ibegin+ni)
+  bounds_lon(:,:) = bounds_lon_glo(:,1+ibegin:1+ibegin+ni)
+  bounds_lat(:,:) = bounds_lat_glo(:,1+ibegin:1+ibegin+ni)
+  field_temp(:,:) = rank
+
+
+!!! Context ATMOSPHERE
+
+  CALL xios_context_initialize("atmosphere",comm)
+  CALL xios_get_handle("atmosphere",ctx_hdl)
+  CALL xios_set_current_context(ctx_hdl)
+
+  CALL xios_define_calendar(type="Gregorian", &
+                            start_date=xios_date(2000, 01, 01, 00, 00, 00), &
+                            time_origin=xios_date(1999, 01, 01, 15, 00, 00))
+
+  CALL xios_set_domain_attr("face", ni_glo=ncell, ibegin=ibegin, ni=ni, type='unstructured')
+  CALL xios_set_domain_attr("face", lonvalue_1d=lon, latvalue_1d=lat)
+  CALL xios_set_domain_attr("face", bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat)
+
+   
+!!! Definition du timestep
+
+  dtime%second=3600
+  CALL xios_set_timestep(timestep=dtime)
+
+!!! Fin de la definition du contexte
+
+  CALL xios_close_context_definition()
+
+
+!####################################################################################
+!!! Boucle temporelle
+!####################################################################################
+    
+    DO ts=1,1
+
+!!! Mise a jour du pas de temps
+      CALL xios_update_calendar(ts)
+
+!!! On donne la valeur du champ atm
+     CALL xios_send_field("temp",field_temp(:,1))
+
+    ENDDO
+
+!####################################################################################
+!!! Finalisation
+!####################################################################################
+
+!!! Fin des contextes
+
+    CALL xios_context_finalize()
+
+    DEALLOCATE(lon_glo, lat_glo)
+    DEALLOCATE(bounds_lon_glo, bounds_lat_glo)
+    DEALLOCATE(field_temp_glo)
+    DEALLOCATE(lon, lat)
+    DEALLOCATE(bounds_lon, bounds_lat)
+    DEALLOCATE(field_temp)
+
+!!! Fin de XIOS
+
+    CALL MPI_COMM_FREE(comm, ierr)
+
+    CALL xios_finalize()
+
+    CALL MPI_FINALIZE(ierr)
+
+  END PROGRAM test_regular
+
+
+
+
+
+
