Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_duplicate_scalar.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_duplicate_scalar.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_duplicate_scalar.cpp	(revision 1985)
@@ -0,0 +1,78 @@
+/*!
+   \file axis_algorithm_duplicate_scalar.cpp
+
+   \brief Algorithm to duplicate scalar into  axis 
+ */
+#include "axis_algorithm_duplicate_scalar.hpp"
+#include "duplicate_scalar_to_axis.hpp"
+#include "axis.hpp"
+#include "scalar.hpp"
+#include "grid.hpp"
+#include "grid_transformation_factory_impl.hpp"
+
+namespace xios {
+CGenericAlgorithmTransformation* CAxisAlgorithmDuplicateScalar::create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                                     CTransformation<CAxis>* transformation,
+                                                                     int elementPositionInGrid,
+                                                                     std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                                     std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                                     std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2DomainPosition)
+TRY
+{
+  std::vector<CAxis*> axisListDestP = gridDst->getAxis();
+  std::vector<CScalar*> scalarListSrcP = gridSrc->getScalars();
+
+  CDuplicateScalarToAxis* duplicateScalar = dynamic_cast<CDuplicateScalarToAxis*> (transformation);
+  int axisDstIndex = elementPositionInGridDst2AxisPosition[elementPositionInGrid];
+  int scalarSrcIndex = elementPositionInGridSrc2ScalarPosition[elementPositionInGrid];
+
+  return (new CAxisAlgorithmDuplicateScalar(axisListDestP[axisDstIndex], scalarListSrcP[scalarSrcIndex], duplicateScalar));
+}
+CATCH
+
+bool CAxisAlgorithmDuplicateScalar::dummyRegistered_ = CAxisAlgorithmDuplicateScalar::registerTrans();
+bool CAxisAlgorithmDuplicateScalar::registerTrans()
+TRY
+{
+  return CGridTransformationFactory<CAxis>::registerTransformation(TRANS_DUPLICATE_SCALAR_TO_AXIS, create);
+}
+CATCH
+
+
+CAxisAlgorithmDuplicateScalar::CAxisAlgorithmDuplicateScalar(CAxis* axisDestination, CScalar* scalarSource, CDuplicateScalarToAxis* algo)
+ : CAxisAlgorithmTransformation(axisDestination, scalarSource)
+{
+
+}
+
+
+CAxisAlgorithmDuplicateScalar::~CAxisAlgorithmDuplicateScalar()
+{
+}
+
+void CAxisAlgorithmDuplicateScalar::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+TRY
+{
+  this->transformationMapping_.resize(1);
+  this->transformationWeight_.resize(1);
+
+  TransformationIndexMap& transMap = this->transformationMapping_[0];
+  TransformationWeightMap& transWeight = this->transformationWeight_[0];
+
+  CArray<int,1>& axisDstIndex = axisDest_->index;
+
+  int nbAxisIdx = axisDstIndex.numElements();
+  for (int idxAxis = 0; idxAxis < nbAxisIdx; ++idxAxis)
+  {
+    int globalAxisIdx = axisDstIndex(idxAxis);
+    transMap[globalAxisIdx].resize(1);
+    transWeight[globalAxisIdx].resize(1);
+    transMap[globalAxisIdx][0] = 0 ;
+    transWeight[globalAxisIdx][0] = 1.0 ;
+  }
+}
+CATCH
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_duplicate_scalar.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_duplicate_scalar.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_duplicate_scalar.hpp	(revision 1985)
@@ -0,0 +1,52 @@
+/*!
+   \file axis_algorithm_reduce_domain.hpp
+   \author Ha NGUYEN
+   \since 23 June 2016
+   \date 23 June 2016
+
+   \brief Algorithm for reduce an axis to a axis
+ */
+#ifndef __XIOS_AXIS_ALGORITHM_DUPLICATE_SCALAR_HPP__
+#define __XIOS_AXIS_ALGORITHM_DUPLICATE_SCALAR_HPP__
+
+#include "axis_algorithm_transformation.hpp"
+#include "transformation.hpp"
+
+namespace xios {
+
+class CAxis;
+class CScalar;
+class CDuplicateScalarToAxis;
+
+
+/*!
+  \class CAxisAlgorithmDuplicateScalar
+  Duplicate scalar into axis destination
+*/
+class CAxisAlgorithmDuplicateScalar : public CAxisAlgorithmTransformation
+{
+public:
+  CAxisAlgorithmDuplicateScalar(CAxis* axisDestination, CScalar* scalarSource, CDuplicateScalarToAxis* algo);
+
+  virtual ~CAxisAlgorithmDuplicateScalar();
+
+  static bool registerTrans();
+
+protected:
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs);
+
+private:
+  static CGenericAlgorithmTransformation* create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                CTransformation<CAxis>* transformation,
+                                                int elementPositionInGrid,
+                                                std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridDst2DomainPosition);
+  static bool dummyRegistered_;
+};
+
+}
+#endif // __XIOS_AXIS_ALGORITHM_DUPLICATE_SCALAR_HPP__
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_extract.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_extract.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_extract.cpp	(revision 1985)
@@ -0,0 +1,122 @@
+/*!
+   \file axis_algorithm_extract.cpp
+   \brief Algorithm for extracting an axis.
+ */
+#include "axis_algorithm_extract.hpp"
+#include "axis.hpp"
+#include "grid.hpp"
+#include "grid_transformation_factory_impl.hpp"
+#include "extract_axis.hpp"
+
+namespace xios {
+CGenericAlgorithmTransformation* CAxisAlgorithmExtract::create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                           CTransformation<CAxis>* transformation,
+                                                           int elementPositionInGrid,
+                                                           std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                           std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                           std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                           std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                           std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                           std::map<int, int>& elementPositionInGridDst2DomainPosition)
+TRY
+{
+  std::vector<CAxis*> axisListDestP = gridDst->getAxis();
+  std::vector<CAxis*> axisListSrcP  = gridSrc->getAxis();
+
+  CExtractAxis* extractAxis = dynamic_cast<CExtractAxis*> (transformation);
+  int axisDstIndex = elementPositionInGridDst2AxisPosition[elementPositionInGrid];
+  int axisSrcIndex = elementPositionInGridSrc2AxisPosition[elementPositionInGrid];
+
+  return (new CAxisAlgorithmExtract(axisListDestP[axisDstIndex], axisListSrcP[axisSrcIndex], extractAxis));
+}
+CATCH
+
+bool CAxisAlgorithmExtract::dummyRegistered_ = CAxisAlgorithmExtract::registerTrans();
+bool CAxisAlgorithmExtract::registerTrans()
+TRY
+{
+  return CGridTransformationFactory<CAxis>::registerTransformation(TRANS_EXTRACT_AXIS, create);
+}
+CATCH
+
+CAxisAlgorithmExtract::CAxisAlgorithmExtract(CAxis* axisDestination, CAxis* axisSource, CExtractAxis* extractAxis)
+: CAxisAlgorithmTransformation(axisDestination, axisSource)
+TRY
+{
+  extractAxis->checkValid(axisSource);
+  extractBegin_ = extractAxis->begin.getValue();
+  extractN_  = extractAxis->n.getValue();
+  extractEnd_   = extractBegin_ + extractN_ - 1;
+
+  if (extractN_ > axisSource->n_glo.getValue())
+  {
+    ERROR("CAxisAlgorithmExtract::CAxisAlgorithmExtract(CAxis* axisDestination, CAxis* axisSource, CExtractAxis* extractAxis)",
+           << "Extract size is greater than global size of source axis"
+           << "Global size of source axis " <<axisSource->getId() << " is " << axisSource->n_glo.getValue()  << std::endl
+           << "Extract size is " << extractN_ );
+  }
+
+  int idxSrc, nDest = 0, beginDestLoc, beginDestGlo = 0 ;
+  int indGloDest, indGloSrc, iSrc;
+  for (int i = 0; i < axisSrc_->n.getValue(); i++)
+  {
+    idxSrc = axisSrc_->index(i);
+    if ((idxSrc >= extractBegin_) && (idxSrc <= extractEnd_))
+    {
+      if (nDest == 0) beginDestLoc = i;
+      ++nDest;
+    }
+  }
+  beginDestGlo = beginDestLoc + axisSrc_->begin - extractBegin_;
+  axisDest_->n_glo.setValue(extractN_);
+  axisDest_->n.setValue(nDest);
+  axisDest_->begin.setValue(beginDestGlo);
+  axisDest_->index.resize(nDest);
+
+  axisDest_->data_n.setValue(nDest);
+  axisDest_->data_begin.setValue(0);
+  axisDest_->data_index.resize(nDest);
+
+  axisDest_->mask.resize(nDest);
+  if (axisSrc_->hasValue) axisDest_->value.resize(nDest);
+  if (axisSrc_->hasLabel) axisDest_->label.resize(nDest);
+  if (axisSrc_->hasBounds) axisDest_->bounds.resize(2,nDest);
+
+  this->transformationMapping_.resize(1);
+  this->transformationWeight_.resize(1);
+  TransformationIndexMap& transMap = this->transformationMapping_[0];
+  TransformationWeightMap& transWeight = this->transformationWeight_[0];
+
+  for (int iDest = 0; iDest < nDest; iDest++)
+  {
+    iSrc = iDest + beginDestLoc;
+    axisDest_->index(iDest) = iDest + beginDestGlo;
+    axisDest_->data_index(iDest) = axisSrc_->data_index(iSrc) - beginDestLoc;
+    axisDest_->mask(iDest) = axisSrc_->mask(iSrc);
+
+    if (axisSrc_->hasValue)
+      axisDest_->value(iDest) = axisSrc_->value(iSrc);
+    if (axisSrc_->hasLabel)
+      axisDest_->label(iDest) = axisSrc_->label(iSrc);
+    if (axisSrc_->hasBounds)
+    {
+      axisDest_->bounds(0,iDest) = axisSrc_->bounds(0,iSrc);
+      axisDest_->bounds(1,iDest) = axisSrc_->bounds(1,iSrc);
+    }
+    indGloDest = axisDest_->index(iDest);
+    indGloSrc = axisSrc_->index(iSrc);
+    transMap[indGloDest].push_back(indGloSrc);
+    transWeight[indGloDest].push_back(1.0);
+
+  }
+}
+CATCH
+
+/*!
+  Compute the index mapping between domain on grid source and one on grid destination
+*/
+void CAxisAlgorithmExtract::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+{
+}
+
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_extract.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_extract.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_extract.hpp	(revision 1985)
@@ -0,0 +1,61 @@
+/*!
+   \file axis_algorithm_extract.hpp
+   \brief Algorithm for extracting an axis.
+ */
+#ifndef __XIOS_AXIS_ALGORITHM_EXTRACT_HPP__
+#define __XIOS_AXIS_ALGORITHM_EXTRACT_HPP__
+
+#include "axis_algorithm_transformation.hpp"
+#include "transformation.hpp"
+
+namespace xios {
+class CAxis;
+class CExtractAxis;
+
+/*!
+  \class CAxisAlgorithmExtract
+  Implementing extract on axis
+  A extracted region can be considered as region that isn't masked.
+  Only this extracted region is extracted to write on Netcdf.
+*/
+class CAxisAlgorithmExtract : public CAxisAlgorithmTransformation
+{
+public:
+  CAxisAlgorithmExtract(CAxis* axisDestination, CAxis* axisSource, CExtractAxis* extractAxis);
+
+  virtual ~CAxisAlgorithmExtract() {}
+
+  static bool registerTrans();
+protected:
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs);
+
+private:
+
+private:
+  //! Global extract begin on axis
+  StdSize extractBegin_;
+
+  //! Global extract end on axis
+  StdSize extractEnd_;
+
+  //! Global extract size on axis
+  StdSize extractN_;
+
+  std::vector<int> extractIndex_;
+
+private:
+
+  static CGenericAlgorithmTransformation* create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                CTransformation<CAxis>* transformation,
+                                                int elementPositionInGrid,
+                                                std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridDst2DomainPosition);
+  static bool dummyRegistered_;
+};
+
+}
+#endif // __XIOS_AXIS_ALGORITHM_EXTRACT_HPP__
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_extract_domain.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_extract_domain.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_extract_domain.cpp	(revision 1985)
@@ -0,0 +1,131 @@
+/*!
+   \file axis_algorithm_reduce_domain.cpp
+   \author Ha NGUYEN
+   \since 23 June 2016
+   \date 23 June 2016
+
+   \brief Algorithm for extract a domain to an axis
+ */
+#include "axis_algorithm_extract_domain.hpp"
+#include "extract_domain_to_axis.hpp"
+#include "axis.hpp"
+#include "domain.hpp"
+#include "grid.hpp"
+#include "grid_transformation_factory_impl.hpp"
+#include "reduction.hpp"
+
+namespace xios {
+CGenericAlgorithmTransformation* CAxisAlgorithmExtractDomain::create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                                     CTransformation<CAxis>* transformation,
+                                                                     int elementPositionInGrid,
+                                                                     std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                                     std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                                     std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2DomainPosition)
+TRY
+{
+  std::vector<CAxis*> axisListDestP = gridDst->getAxis();
+  std::vector<CDomain*> domainListSrcP = gridSrc->getDomains();
+
+  CExtractDomainToAxis* extractDomain = dynamic_cast<CExtractDomainToAxis*> (transformation);
+  int axisDstIndex = elementPositionInGridDst2AxisPosition[elementPositionInGrid];
+  int domainSrcIndex = elementPositionInGridSrc2DomainPosition[elementPositionInGrid];
+
+  return (new CAxisAlgorithmExtractDomain(axisListDestP[axisDstIndex], domainListSrcP[domainSrcIndex], extractDomain));
+}
+CATCH
+
+bool CAxisAlgorithmExtractDomain::dummyRegistered_ = CAxisAlgorithmExtractDomain::registerTrans();
+bool CAxisAlgorithmExtractDomain::registerTrans()
+TRY
+{
+  return CGridTransformationFactory<CAxis>::registerTransformation(TRANS_EXTRACT_DOMAIN_TO_AXIS, create);
+}
+CATCH
+
+
+CAxisAlgorithmExtractDomain::CAxisAlgorithmExtractDomain(CAxis* axisDestination, CDomain* domainSource, CExtractDomainToAxis* algo)
+ : CAxisAlgorithmTransformation(axisDestination, domainSource), pos_(-1), reduction_(0)
+TRY
+{
+  algo->checkValid(axisDestination, domainSource);
+  StdString op = "extract";
+
+  switch (algo->direction)
+  {
+    case CExtractDomainToAxis::direction_attr::jDir:
+      dir_ = jDir;
+      break;
+    case CExtractDomainToAxis::direction_attr::iDir:
+      dir_ = iDir;
+      break;
+    default:
+      break;
+  }
+
+  pos_ = algo->position;
+  reduction_ = CReductionAlgorithm::createOperation(CReductionAlgorithm::ReductionOperations[op]);
+}
+CATCH
+
+void CAxisAlgorithmExtractDomain::apply(const std::vector<std::pair<int,double> >& localIndex,
+                                        const double* dataInput,
+                                        CArray<double,1>& dataOut,
+                                        std::vector<bool>& flagInitial,                     
+                                        bool ignoreMissingValue, bool firstPass)
+TRY
+{
+  reduction_->apply(localIndex, dataInput, dataOut, flagInitial, ignoreMissingValue, firstPass);
+}
+CATCH
+
+CAxisAlgorithmExtractDomain::~CAxisAlgorithmExtractDomain()
+TRY
+{
+  if (0 != reduction_) delete reduction_;
+}
+CATCH
+
+void CAxisAlgorithmExtractDomain::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+TRY
+{
+  this->transformationMapping_.resize(1);
+  this->transformationWeight_.resize(1);
+
+  TransformationIndexMap& transMap = this->transformationMapping_[0];
+  TransformationWeightMap& transWeight = this->transformationWeight_[0];
+
+  CArray<int,1>& axisDstIndex = axisDest_->index;
+  int ni_glo = domainSrc_->ni_glo, nj_glo = domainSrc_->nj_glo;
+  if (jDir == dir_)
+  {
+    int nbAxisIdx = axisDstIndex.numElements();
+    for (int idxAxis = 0; idxAxis < nbAxisIdx; ++idxAxis)
+    {
+      int globalAxisIdx = axisDstIndex(idxAxis);
+      transMap[globalAxisIdx].resize(1);
+      transWeight[globalAxisIdx].resize(1);
+      transMap[globalAxisIdx][0] = globalAxisIdx * ni_glo + pos_;
+      transWeight[globalAxisIdx][0] = 1.0;
+
+    }
+  }
+  else if (iDir == dir_)
+  {
+    int nbAxisIdx = axisDstIndex.numElements();
+    for (int idxAxis = 0; idxAxis < nbAxisIdx; ++idxAxis)
+    {
+      int globalAxisIdx = axisDstIndex(idxAxis);
+      transMap[globalAxisIdx].resize(1);
+      transWeight[globalAxisIdx].resize(1);
+      transMap[globalAxisIdx][0] = globalAxisIdx + ni_glo * pos_;
+      transWeight[globalAxisIdx][0] = 1.0;
+    }
+  }
+  else
+  {}
+}
+CATCH
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_extract_domain.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_extract_domain.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_extract_domain.hpp	(revision 1985)
@@ -0,0 +1,70 @@
+/*!
+   \file axis_algorithm_reduce_domain.hpp
+   \author Ha NGUYEN
+   \since 23 June 2016
+   \date 23 June 2016
+
+   \brief Algorithm for reduce an axis to a axis
+ */
+#ifndef __XIOS_AXIS_ALGORITHM_EXTRACT_DOMAIN_HPP__
+#define __XIOS_AXIS_ALGORITHM_EXTRACT_DOMAIN_HPP__
+
+#include "axis_algorithm_transformation.hpp"
+#include "transformation.hpp"
+
+namespace xios {
+
+class CAxis;
+class CDomain;
+class CExtractDomainToAxis;
+class CReductionAlgorithm;
+
+
+/*!
+  \class CAxisAlgorithmExtractDomain
+  Extract a domain to an axis
+*/
+class CAxisAlgorithmExtractDomain : public CAxisAlgorithmTransformation
+{
+public:
+  CAxisAlgorithmExtractDomain(CAxis* axisDestination, CDomain* domainSource, CExtractDomainToAxis* algo);
+
+  virtual void apply(const std::vector<std::pair<int,double> >& localIndex,
+                     const double* dataInput,
+                     CArray<double,1>& dataOut,
+                     std::vector<bool>& flagInitial,                     
+                     bool ignoreMissingValue, bool firstPass);
+
+  virtual ~CAxisAlgorithmExtractDomain();
+
+  static bool registerTrans();
+protected:
+  enum ExtractDirection {
+    undefined = 0,
+    iDir = 1,
+    jDir = 2
+  };
+
+  ExtractDirection dir_;
+  int pos_; //! Position to extract
+protected:
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs);
+
+protected:
+  CReductionAlgorithm* reduction_;
+
+private:
+  static CGenericAlgorithmTransformation* create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                CTransformation<CAxis>* transformation,
+                                                int elementPositionInGrid,
+                                                std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridDst2DomainPosition);
+  static bool dummyRegistered_;
+};
+
+}
+#endif // __XIOS_AXIS_ALGORITHM_EXTRACT_DOMAIN_HPP__
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_interpolate.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_interpolate.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_interpolate.cpp	(revision 1985)
@@ -0,0 +1,394 @@
+/*!
+   \file axis_algorithm_interpolate.cpp
+   \author Ha NGUYEN
+   \since 23 June 2015
+   \date 02 Jul 2015
+
+   \brief Algorithm for interpolation on an axis.
+ */
+#include "axis_algorithm_interpolate.hpp"
+#include "axis.hpp"
+#include "interpolate_axis.hpp"
+#include <algorithm>
+#include "context.hpp"
+#include "context_client.hpp"
+#include "utils.hpp"
+#include "grid.hpp"
+#include "grid_transformation_factory_impl.hpp"
+#include "distribution_client.hpp"
+#include "timer.hpp"
+
+namespace xios {
+CGenericAlgorithmTransformation* CAxisAlgorithmInterpolate::create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                                   CTransformation<CAxis>* transformation,
+                                                                   int elementPositionInGrid,
+                                                                   std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                                   std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                                   std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                                   std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                                   std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                                   std::map<int, int>& elementPositionInGridDst2DomainPosition)
+TRY
+{
+  std::vector<CAxis*> axisListDestP = gridDst->getAxis();
+  std::vector<CAxis*> axisListSrcP  = gridSrc->getAxis();
+
+  CInterpolateAxis* interpolateAxis = dynamic_cast<CInterpolateAxis*> (transformation);
+  int axisDstIndex = elementPositionInGridDst2AxisPosition[elementPositionInGrid];
+  int axisSrcIndex = elementPositionInGridSrc2AxisPosition[elementPositionInGrid];
+
+  return (new CAxisAlgorithmInterpolate(axisListDestP[axisDstIndex], axisListSrcP[axisSrcIndex], interpolateAxis));
+}
+CATCH
+
+bool CAxisAlgorithmInterpolate::dummyRegistered_ = CAxisAlgorithmInterpolate::registerTrans();
+bool CAxisAlgorithmInterpolate::registerTrans()
+TRY
+{
+  return CGridTransformationFactory<CAxis>::registerTransformation(TRANS_INTERPOLATE_AXIS, create);
+}
+CATCH
+
+CAxisAlgorithmInterpolate::CAxisAlgorithmInterpolate(CAxis* axisDestination, CAxis* axisSource, CInterpolateAxis* interpAxis)
+: CAxisAlgorithmTransformation(axisDestination, axisSource), coordinate_(), transPosition_()
+TRY
+{
+  interpAxis->checkValid(axisSource);
+  order_ = interpAxis->order.getValue();
+  if (!interpAxis->coordinate.isEmpty())
+  {
+    coordinate_ = interpAxis->coordinate.getValue();
+    this->idAuxInputs_.resize(1);
+    this->idAuxInputs_[0] = coordinate_;
+  }
+}
+CATCH
+
+/*!
+  Compute the index mapping between axis on grid source and one on grid destination
+*/
+void CAxisAlgorithmInterpolate::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+TRY
+{
+  CTimer::get("CAxisAlgorithmInterpolate::computeIndexSourceMapping_").resume() ;
+  CContext* context = CContext::getCurrent();
+  int nbClient = context->intraCommSize_;
+  CArray<bool,1>& axisMask = axisSrc_->mask;
+  int srcSize  = axisSrc_->n_glo.getValue();
+  std::vector<CArray<double,1> > vecAxisValue;
+
+  // Fill in axis value from coordinate
+  fillInAxisValue(vecAxisValue, dataAuxInputs);
+  std::vector<double> valueSrc(srcSize);
+  std::vector<double> recvBuff(srcSize);
+  std::vector<int> indexVec(srcSize);
+
+  for (int idx = 0; idx < vecAxisValue.size(); ++idx)
+  {
+    CArray<double,1>& axisValue = vecAxisValue[idx];
+    retrieveAllAxisValue(axisValue, axisMask, recvBuff, indexVec);
+    XIOSAlgorithms::sortWithIndex<double, CVectorStorage>(recvBuff, indexVec);
+    for (int i = 0; i < srcSize; ++i) valueSrc[i] = recvBuff[indexVec[i]];
+    computeInterpolantPoint(valueSrc, indexVec, idx);
+  }
+  CTimer::get("CAxisAlgorithmInterpolate::computeIndexSourceMapping_").suspend() ;
+}
+CATCH
+
+/*!
+  Compute the interpolant points
+  Assume that we have all value of axis source, with these values, need to calculate weight (coeff) of Lagrange polynomial
+  \param [in] axisValue all value of axis source
+  \param [in] tranPos position of axis on a domain
+*/
+void CAxisAlgorithmInterpolate::computeInterpolantPoint(const std::vector<double>& axisValue,
+                                                        const std::vector<int>& indexVec,
+                                                        int transPos)
+TRY
+{
+  std::vector<double>::const_iterator itb = axisValue.begin(), ite = axisValue.end();
+  std::vector<double>::const_iterator itLowerBound, itUpperBound, it, iteRange, itfirst, itsecond;
+  const double sfmax = NumTraits<double>::sfmax();
+  const double precision = NumTraits<double>::dummy_precision();
+
+  int ibegin = axisDest_->begin.getValue();
+  CArray<double,1>& axisDestValue = axisDest_->value;
+  int numValue = axisDestValue.numElements();
+  std::map<int, std::vector<std::pair<int,double> > > interpolatingIndexValues;
+
+  for (int idx = 0; idx < numValue; ++idx)
+  {
+    bool outOfRange = false;
+    double destValue = axisDestValue(idx);
+    if (destValue < *itb) outOfRange = true;
+
+    itLowerBound = std::lower_bound(itb, ite, destValue);
+    itUpperBound = std::upper_bound(itb, ite, destValue);
+    if ((ite != itUpperBound) && (sfmax == *itUpperBound)) itUpperBound = ite;
+
+    if ((ite == itLowerBound) || (ite == itUpperBound)) outOfRange = true;
+
+    // We don't do extrapolation FOR NOW, maybe in the future
+    if (!outOfRange)
+    {
+      if ((itLowerBound == itUpperBound) && (itb != itLowerBound)) --itLowerBound;
+      double distanceToLower = destValue - *itLowerBound;
+      double distanceToUpper = *itUpperBound - destValue;
+      int order = (order_ + 1) - 2;
+      bool down = (distanceToLower < distanceToUpper) ? true : false;
+      for (int k = 0; k < order; ++k)
+      {
+        if ((itb != itLowerBound) && down)
+        {
+          --itLowerBound;
+          distanceToLower = destValue - *itLowerBound;
+          down = (distanceToLower < distanceToUpper) ? true : false;
+          continue;
+        }
+        if ((ite != itUpperBound) && (sfmax != *itUpperBound))
+        {
+          ++itUpperBound;
+          distanceToUpper = *itUpperBound - destValue;
+          down = (distanceToLower < distanceToUpper) ? true : false;
+
+        }
+      }
+
+      iteRange = (ite == itUpperBound) ? itUpperBound : itUpperBound + 1;
+      itsecond = it = itLowerBound; ++itsecond;
+      while (it < iteRange)
+      {
+        while ( (itsecond < ite) && ((*itsecond -*it) < precision) )
+        { ++itsecond; ++it; }
+        int index = std::distance(itb, it);
+        interpolatingIndexValues[idx+ibegin].push_back(make_pair(indexVec[index],*it));
+        ++it; ++itsecond;
+      }
+
+    }
+  }
+  computeWeightedValueAndMapping(interpolatingIndexValues, transPos);
+}
+CATCH
+
+/*!
+  Compute weight (coeff) of Lagrange's polynomial
+  \param [in] interpolatingIndexValues the necessary axis value to calculate the coeffs
+*/
+void CAxisAlgorithmInterpolate::computeWeightedValueAndMapping(const std::map<int, std::vector<std::pair<int,double> > >& interpolatingIndexValues, int transPos)
+TRY
+{
+  TransformationIndexMap& transMap = this->transformationMapping_[transPos];
+  TransformationWeightMap& transWeight = this->transformationWeight_[transPos];
+  std::map<int, std::vector<std::pair<int,double> > >::const_iterator itb = interpolatingIndexValues.begin(), it,
+                                                                      ite = interpolatingIndexValues.end();
+  int ibegin = axisDest_->begin.getValue();
+  for (it = itb; it != ite; ++it)
+  {
+    int globalIndexDest = it->first;
+    double localValue = axisDest_->value(globalIndexDest - ibegin);
+    const std::vector<std::pair<int,double> >& interpVal = it->second;
+    int interpSize = interpVal.size();
+    transMap[globalIndexDest].resize(interpSize);
+    transWeight[globalIndexDest].resize(interpSize);
+    for (int idx = 0; idx < interpSize; ++idx)
+    {
+      int index = interpVal[idx].first;
+      double weight = 1.0;
+
+      for (int k = 0; k < interpSize; ++k)
+      {
+        if (k == idx) continue;
+        weight *= (localValue - interpVal[k].second);
+        weight /= (interpVal[idx].second - interpVal[k].second);
+      }
+      transMap[globalIndexDest][idx] = index;
+      transWeight[globalIndexDest][idx] = weight;
+      if (!transPosition_.empty())
+      {
+        (this->transformationPosition_[transPos])[globalIndexDest] = transPosition_[transPos];
+      }
+    }
+  }
+  if (!transPosition_.empty() && this->transformationPosition_[transPos].empty())
+    (this->transformationPosition_[transPos])[0] = transPosition_[transPos];
+
+}
+CATCH
+
+/*!
+  Each client retrieves all values of an axis
+  \param [in/out] recvBuff buffer for receiving values (already allocated)
+  \param [in/out] indexVec mapping between values and global index of axis
+*/
+void CAxisAlgorithmInterpolate::retrieveAllAxisValue(const CArray<double,1>& axisValue, const CArray<bool,1>& axisMask,
+                                                     std::vector<double>& recvBuff, std::vector<int>& indexVec)
+TRY
+{
+  CContext* context = CContext::getCurrent();
+  int nbClient = context->intraCommSize_;
+
+  int srcSize  = axisSrc_->n_glo.getValue();
+  int numValue = axisValue.numElements();
+
+  if (srcSize == numValue)  // Only one client or axis not distributed
+  {
+    for (int idx = 0; idx < srcSize; ++idx)
+    {
+      if (axisMask(idx))
+      {
+        recvBuff[idx] = axisValue(idx);
+        indexVec[idx] = idx;
+      }
+      else
+      {
+        recvBuff[idx] = NumTraits<double>::sfmax();
+        indexVec[idx] = -1;
+      }
+    }
+
+  }
+  else // Axis distributed
+  {
+    double* sendValueBuff = new double [numValue];
+    int* sendIndexBuff = new int [numValue];
+    int* recvIndexBuff = new int [srcSize];
+
+    int ibegin = axisSrc_->begin.getValue();
+    for (int idx = 0; idx < numValue; ++idx)
+    {
+      if (axisMask(idx))
+      {
+        sendValueBuff[idx] = axisValue(idx);
+        sendIndexBuff[idx] = idx + ibegin;
+      }
+      else
+      {
+        sendValueBuff[idx] = NumTraits<double>::sfmax();
+        sendIndexBuff[idx] = -1;
+      }
+    }
+
+    int* recvCount=new int[nbClient];
+    MPI_Allgather(&numValue,1,MPI_INT,recvCount,1,MPI_INT,context->intraComm_);
+
+    int* displ=new int[nbClient];
+    displ[0]=0 ;
+    for(int n=1;n<nbClient;n++) displ[n]=displ[n-1]+recvCount[n-1];
+
+    // Each client have enough global info of axis
+    MPI_Allgatherv(sendIndexBuff,numValue,MPI_INT,recvIndexBuff,recvCount,displ,MPI_INT,context->intraComm_);
+    MPI_Allgatherv(sendValueBuff,numValue,MPI_DOUBLE,&(recvBuff[0]),recvCount,displ,MPI_DOUBLE,context->intraComm_);
+
+    for (int idx = 0; idx < srcSize; ++idx)
+    {
+      indexVec[idx] = recvIndexBuff[idx];
+    }
+
+    delete [] displ;
+    delete [] recvCount;
+    delete [] recvIndexBuff;
+    delete [] sendIndexBuff;
+    delete [] sendValueBuff;
+  }
+}
+CATCH
+
+/*!
+  Fill in axis value dynamically from a field whose grid is composed of a domain and an axis
+  \param [in/out] vecAxisValue vector axis value filled in from input field
+*/
+void CAxisAlgorithmInterpolate::fillInAxisValue(std::vector<CArray<double,1> >& vecAxisValue,
+                                                const std::vector<CArray<double,1>* >& dataAuxInputs)
+TRY
+{
+  if (coordinate_.empty())
+  {
+    vecAxisValue.resize(1);
+    vecAxisValue[0].resize(axisSrc_->value.numElements());
+    vecAxisValue[0] = axisSrc_->value;
+    this->transformationMapping_.resize(1);
+    this->transformationWeight_.resize(1);
+  }
+  else
+  {
+    CField* field = CField::get(coordinate_);
+    CGrid* grid = field->getGrid();
+
+    std::vector<CDomain*> domListP = grid->getDomains();
+    std::vector<CAxis*> axisListP = grid->getAxis();
+    if (domListP.empty() || axisListP.empty() || (1 < domListP.size()) || (1 < axisListP.size()))
+      ERROR("CAxisAlgorithmInterpolate::fillInAxisValue(std::vector<CArray<double,1> >& vecAxisValue)",
+             << "XIOS only supports dynamic interpolation with coordinate (field) associated with grid composed of a domain and an axis"
+             << "Coordinate (field) id = " <<field->getId() << std::endl
+             << "Associated grid id = " << grid->getId());
+
+    CDomain* dom = domListP[0];
+    size_t vecAxisValueSize = dom->i_index.numElements();
+    size_t vecAxisValueSizeWithMask = 0;
+    for (size_t idx = 0; idx < vecAxisValueSize; ++idx)
+    {
+      if (dom->domainMask(idx)) ++vecAxisValueSizeWithMask;
+    }
+
+    int niGlobDom = dom->ni_glo.getValue();
+    vecAxisValue.resize(vecAxisValueSizeWithMask);
+    if (transPosition_.empty())
+    {
+      size_t indexMask = 0;
+      transPosition_.resize(vecAxisValueSizeWithMask);
+      for (size_t idx = 0; idx < vecAxisValueSize; ++idx)
+      {
+        if (dom->domainMask(idx))
+        {
+          transPosition_[indexMask].resize(1);
+          transPosition_[indexMask][0] = (dom->i_index)(idx) + niGlobDom * (dom->j_index)(idx);
+          ++indexMask;
+        }
+
+      }
+    }
+    this->transformationMapping_.resize(vecAxisValueSizeWithMask);
+    this->transformationWeight_.resize(vecAxisValueSizeWithMask);
+    this->transformationPosition_.resize(vecAxisValueSizeWithMask);
+
+    const CDistributionClient::GlobalLocalDataMap& globalLocalIndexSendToServer = grid->getClientDistribution()->getGlobalLocalDataSendToServer();
+    CDistributionClient::GlobalLocalDataMap::const_iterator itIndex, iteIndex = globalLocalIndexSendToServer.end();
+    size_t axisSrcSize = axisSrc_->index.numElements();
+    std::vector<int> globalDimension = grid->getGlobalDimension();
+
+    size_t indexMask = 0;
+    for (size_t idx = 0; idx < vecAxisValueSize; ++idx)
+    {
+      if (dom->domainMask(idx))
+      {
+        size_t axisValueSize = 0;
+        for (size_t jdx = 0; jdx < axisSrcSize; ++jdx)
+        {
+          size_t globalIndex = ((dom->i_index)(idx) + (dom->j_index)(idx)*globalDimension[0]) + (axisSrc_->index)(jdx)*globalDimension[0]*globalDimension[1];
+          if (iteIndex != globalLocalIndexSendToServer.find(globalIndex))
+          {
+            ++axisValueSize;
+          }
+        }
+
+        vecAxisValue[indexMask].resize(axisValueSize);
+        axisValueSize = 0;
+        for (size_t jdx = 0; jdx < axisSrcSize; ++jdx)
+        {
+          size_t globalIndex = ((dom->i_index)(idx) + (dom->j_index)(idx)*globalDimension[0]) + (axisSrc_->index)(jdx)*globalDimension[0]*globalDimension[1];
+          itIndex = globalLocalIndexSendToServer.find(globalIndex);
+          if (iteIndex != itIndex)
+          {
+            vecAxisValue[indexMask](axisValueSize) = (*dataAuxInputs[0])(itIndex->second);
+            ++axisValueSize;
+          }
+        }
+        ++indexMask;
+      }
+    }
+  }
+}
+CATCH
+
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_interpolate.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_interpolate.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_interpolate.hpp	(revision 1985)
@@ -0,0 +1,66 @@
+/*!
+   \file axis_algorithm_interpolate.hpp
+   \author Ha NGUYEN
+   \since 23 June 2015
+   \date 23 June 2015
+
+   \brief Algorithm for interpolation on an axis.
+ */
+#ifndef __XIOS_AXIS_ALGORITHM_INTERPOLATE_HPP__
+#define __XIOS_AXIS_ALGORITHM_INTERPOLATE_HPP__
+
+#include "axis_algorithm_transformation.hpp"
+#include "transformation.hpp"
+
+namespace xios {
+
+class CAxis;
+class CGrid;
+class CInterpolateAxis;
+
+/*!
+  \class CAxisAlgorithmInterpolate
+  Implementing interpolation on axis
+  The values on axis source are assumed monotonic
+*/
+class CAxisAlgorithmInterpolate : public CAxisAlgorithmTransformation
+{
+public:
+  CAxisAlgorithmInterpolate(CAxis* axisDestination, CAxis* axisSource, CInterpolateAxis* interpAxis);
+
+  virtual ~CAxisAlgorithmInterpolate() {}
+
+  static bool registerTrans();
+protected:
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs);
+
+private:
+  void retrieveAllAxisValue(const CArray<double,1>& axisValue, const CArray<bool,1>& axisMask,
+                            std::vector<double>& recvBuff, std::vector<int>& indexVec);
+  void computeInterpolantPoint(const std::vector<double>& recvBuff, const std::vector<int>&, int transPos = 0);
+  void computeWeightedValueAndMapping(const std::map<int, std::vector<std::pair<int,double> > >& interpolatingIndexValues, int transPos = 0);
+  void fillInAxisValue(std::vector<CArray<double,1> >& vecAxisValue,
+                       const std::vector<CArray<double,1>* >& dataAuxInputs);
+
+private:
+  // Interpolation order
+  int order_;
+  StdString coordinate_;
+  std::vector<std::vector<int> > transPosition_;
+
+private:
+  static CGenericAlgorithmTransformation* create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                CTransformation<CAxis>* transformation,
+                                                int elementPositionInGrid,
+                                                std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridDst2DomainPosition);
+  static bool dummyRegistered_;
+};
+
+}
+
+#endif // __XIOS_AXIS_ALGORITHM_INTERPOLATE_HPP__
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_inverse.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_inverse.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_inverse.cpp	(revision 1985)
@@ -0,0 +1,294 @@
+/*!
+   \file axis_algorithm_inverse.hpp
+   \author Ha NGUYEN
+   \since 14 May 2015
+   \date 29 June 2015
+
+   \brief Algorithm for inversing an axis..
+ */
+#include "axis_algorithm_inverse.hpp"
+#include "context.hpp"
+#include "context_client.hpp"
+#include "axis.hpp"
+#include "grid.hpp"
+#include "grid_transformation_factory_impl.hpp"
+#include "inverse_axis.hpp"
+#include "client_client_dht_template.hpp"
+
+namespace xios {
+
+CGenericAlgorithmTransformation* CAxisAlgorithmInverse::create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                               CTransformation<CAxis>* transformation,
+                                                               int elementPositionInGrid,
+                                                               std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                               std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                               std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                               std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                               std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                               std::map<int, int>& elementPositionInGridDst2DomainPosition)
+TRY
+{
+  std::vector<CAxis*> axisListDestP = gridDst->getAxis();
+  std::vector<CAxis*> axisListSrcP  = gridSrc->getAxis();
+
+  CInverseAxis* inverseAxis = dynamic_cast<CInverseAxis*> (transformation);
+  int axisDstIndex = elementPositionInGridDst2AxisPosition[elementPositionInGrid];
+  int axisSrcIndex = elementPositionInGridSrc2AxisPosition[elementPositionInGrid];
+
+  return (new CAxisAlgorithmInverse(axisListDestP[axisDstIndex], axisListSrcP[axisSrcIndex], inverseAxis));
+}
+CATCH
+
+bool CAxisAlgorithmInverse::dummyRegistered_ = CAxisAlgorithmInverse::registerTrans();
+bool CAxisAlgorithmInverse::registerTrans()
+TRY
+{
+  return CGridTransformationFactory<CAxis>::registerTransformation(TRANS_INVERSE_AXIS, create);
+}
+CATCH
+
+CAxisAlgorithmInverse::CAxisAlgorithmInverse(CAxis* axisDestination, CAxis* axisSource, CInverseAxis* inverseAxis)
+ : CAxisAlgorithmTransformation(axisDestination, axisSource)
+TRY
+{
+  if (axisDestination->n_glo.getValue() != axisSource->n_glo.getValue())
+  {
+    ERROR("CAxisAlgorithmInverse::CAxisAlgorithmInverse(CAxis* axisDestination, CAxis* axisSource)",
+           << "Two axis have different global size"
+           << "Size of axis source " <<axisSource->getId() << " is " << axisSource->n_glo.getValue()  << std::endl
+           << "Size of axis destination " <<axisDestination->getId() << " is " << axisDestination->n_glo.getValue());
+  }
+}
+CATCH
+
+void CAxisAlgorithmInverse::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+TRY
+{
+  this->transformationMapping_.resize(1);
+  this->transformationWeight_.resize(1);
+
+  TransformationIndexMap& transMap = this->transformationMapping_[0];
+  TransformationWeightMap& transWeight = this->transformationWeight_[0];
+
+  int globalIndexSize = axisDestGlobalIndex_.size();
+  for (int idx = 0; idx < globalIndexSize; ++idx)
+  {
+    transMap[axisDestGlobalIndex_[idx]].push_back(axisDestGlobalSize_-axisDestGlobalIndex_[idx]-1);
+    transWeight[axisDestGlobalIndex_[idx]].push_back(1.0);
+  }
+
+  int niSrc   = axisSrc_->n.getValue();
+  int sizeSrc = axisSrc_->n_glo.getValue();
+  if (niSrc != sizeSrc) updateAxisValue();
+  else
+  {
+    for (int idx = 0; idx < sizeSrc; ++idx)
+    {
+      axisDest_->value(idx) = axisSrc_->value(sizeSrc-idx-1);
+    }
+  }
+}
+CATCH
+
+/*!
+  Update value on axis after inversing
+  After an axis is inversed, not only the data on it must be inversed but also the value
+*/
+void CAxisAlgorithmInverse::updateAxisValue()
+TRY
+{
+  CContext* context = CContext::getCurrent();
+  int clientRank = context->intraCommRank_;
+  int nbClient = context->intraCommSize_;
+
+  int niSrc     = axisSrc_->n.getValue();
+  int ibeginSrc = axisSrc_->begin.getValue();
+  int nSrc = axisSrc_->index.numElements();
+
+  CClientClientDHTInt::Index2VectorInfoTypeMap globalIndex2ProcRank;
+  for (int idx = 0; idx < nSrc; ++idx)
+  {
+    if ((axisSrc_->mask)(idx))
+    {
+      globalIndex2ProcRank[(axisSrc_->index)(idx)].resize(1);
+      globalIndex2ProcRank[(axisSrc_->index)(idx)][0] = clientRank;
+    }
+  }
+
+  typedef std::unordered_map<size_t, std::vector<double> > GlobalIndexMapFromSrcToDest;
+  GlobalIndexMapFromSrcToDest globalIndexMapFromSrcToDest;
+  TransformationIndexMap& transMap = this->transformationMapping_[0];
+  TransformationIndexMap::const_iterator itb = transMap.begin(), ite = transMap.end(), it;
+  CArray<size_t,1> globalSrcIndex(transMap.size());
+  int localIndex = 0;
+  for (it = itb; it != ite; ++it)
+  {
+    size_t srcIndex = it->second[0];
+    globalIndexMapFromSrcToDest[srcIndex].resize(1);
+    globalIndexMapFromSrcToDest[srcIndex][0] = it->first;
+    globalSrcIndex(localIndex) = srcIndex;
+    ++localIndex;
+  }
+
+  CClientClientDHTInt dhtIndexProcRank(globalIndex2ProcRank, context->intraComm_);
+  dhtIndexProcRank.computeIndexInfoMapping(globalSrcIndex);
+  CClientClientDHTInt::Index2VectorInfoTypeMap& computedGlobalIndexOnProc = dhtIndexProcRank.getInfoIndexMap();
+  std::unordered_map<int, std::vector<size_t> > globalSrcIndexSendToProc;
+  for (int idx = 0; idx < localIndex; ++idx)
+  {
+    size_t tmpIndex = globalSrcIndex(idx);
+    if (1 == computedGlobalIndexOnProc.count(tmpIndex))
+    {
+      std::vector<int>& tmpVec = computedGlobalIndexOnProc[tmpIndex];
+      globalSrcIndexSendToProc[tmpVec[0]].push_back(tmpIndex);
+    }
+  }
+
+  std::unordered_map<int, std::vector<size_t> >::const_iterator itbIndex = globalSrcIndexSendToProc.begin(), itIndex,
+                                                                  iteIndex = globalSrcIndexSendToProc.end();
+  std::map<int,int> sendRankSizeMap,recvRankSizeMap;
+  int connectedClient = globalSrcIndexSendToProc.size();
+  int* recvCount=new int[nbClient];
+  int* displ=new int[nbClient];
+  int* sendRankBuff=new int[connectedClient];
+  int* sendSizeBuff=new int[connectedClient];
+  int n = 0;
+  for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex, ++n)
+  {
+    sendRankBuff[n] = itIndex->first;
+    int sendSize = itIndex->second.size();
+    sendSizeBuff[n] = sendSize;
+    sendRankSizeMap[itIndex->first] = sendSize;
+  }
+  MPI_Allgather(&connectedClient,1,MPI_INT,recvCount,1,MPI_INT,context->intraComm_);
+
+  displ[0]=0 ;
+  for(int n=1;n<nbClient;n++) displ[n]=displ[n-1]+recvCount[n-1];
+  int recvSize=displ[nbClient-1]+recvCount[nbClient-1];
+  int* recvRankBuff=new int[recvSize];
+  int* recvSizeBuff=new int[recvSize];
+  MPI_Allgatherv(sendRankBuff,connectedClient,MPI_INT,recvRankBuff,recvCount,displ,MPI_INT,context->intraComm_);
+  MPI_Allgatherv(sendSizeBuff,connectedClient,MPI_INT,recvSizeBuff,recvCount,displ,MPI_INT,context->intraComm_);
+  for (int i = 0; i < nbClient; ++i)
+  {
+    int currentPos = displ[i];
+    for (int j = 0; j < recvCount[i]; ++j)
+      if (recvRankBuff[currentPos+j] == clientRank)
+      {
+        recvRankSizeMap[i] = recvSizeBuff[currentPos+j];
+      }
+  }
+
+  // Sending global index of grid source to corresponding process as well as the corresponding mask
+  std::vector<MPI_Request> requests;
+  std::vector<MPI_Status> status;
+  std::unordered_map<int, unsigned long* > recvGlobalIndexSrc;
+  std::unordered_map<int, double* > sendValueToDest;
+  for (std::map<int,int>::const_iterator itRecv = recvRankSizeMap.begin(); itRecv != recvRankSizeMap.end(); ++itRecv)
+  {
+    int recvRank = itRecv->first;
+    int recvSize = itRecv->second;
+    recvGlobalIndexSrc[recvRank] = new unsigned long [recvSize];
+    sendValueToDest[recvRank] = new double [recvSize];
+
+    requests.push_back(MPI_Request());
+    MPI_Irecv(recvGlobalIndexSrc[recvRank], recvSize, MPI_UNSIGNED_LONG, recvRank, 46, context->intraComm_, &requests.back());
+  }
+
+  std::unordered_map<int, unsigned long* > sendGlobalIndexSrc;
+  std::unordered_map<int, double* > recvValueFromSrc;
+  for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex)
+  {
+    int sendRank = itIndex->first;
+    int sendSize = sendRankSizeMap[sendRank];
+    const std::vector<size_t>& sendIndexMap = itIndex->second;
+    std::vector<size_t>::const_iterator itbSend = sendIndexMap.begin(), iteSend = sendIndexMap.end(), itSend;
+    sendGlobalIndexSrc[sendRank] = new unsigned long [sendSize];
+    recvValueFromSrc[sendRank] = new double [sendSize];
+    int countIndex = 0;
+    for (itSend = itbSend; itSend != iteSend; ++itSend)
+    {
+      sendGlobalIndexSrc[sendRank][countIndex] = *itSend;
+      ++countIndex;
+    }
+
+    // Send global index source and mask
+    requests.push_back(MPI_Request());
+    MPI_Isend(sendGlobalIndexSrc[sendRank], sendSize, MPI_UNSIGNED_LONG, sendRank, 46, context->intraComm_, &requests.back());
+  }
+
+  status.resize(requests.size());
+  MPI_Waitall(requests.size(), &requests[0], &status[0]);
+
+
+  std::vector<MPI_Request>().swap(requests);
+  std::vector<MPI_Status>().swap(status);
+
+  // Okie, on destination side, we will wait for information of masked index of source
+  for (std::map<int,int>::const_iterator itSend = sendRankSizeMap.begin(); itSend != sendRankSizeMap.end(); ++itSend)
+  {
+    int recvRank = itSend->first;
+    int recvSize = itSend->second;
+
+    requests.push_back(MPI_Request());
+    MPI_Irecv(recvValueFromSrc[recvRank], recvSize, MPI_DOUBLE, recvRank, 48, context->intraComm_, &requests.back());
+  }
+
+  for (std::map<int,int>::const_iterator itRecv = recvRankSizeMap.begin(); itRecv != recvRankSizeMap.end(); ++itRecv)
+  {
+    int recvRank = itRecv->first;
+    int recvSize = itRecv->second;
+    double* sendValue = sendValueToDest[recvRank];
+    unsigned long* recvIndexSrc = recvGlobalIndexSrc[recvRank];
+    int realSendSize = 0;
+    for (int idx = 0; idx < recvSize; ++idx)
+    {
+      size_t globalIndex = *(recvIndexSrc+idx);
+      int localIndex = globalIndex - ibeginSrc;
+      *(sendValue + idx) = axisSrc_->value(localIndex);
+    }
+    // Okie, now inform the destination which source index are masked
+    requests.push_back(MPI_Request());
+    MPI_Isend(sendValueToDest[recvRank], recvSize, MPI_DOUBLE, recvRank, 48, context->intraComm_, &requests.back());
+  }
+  status.resize(requests.size());
+  MPI_Waitall(requests.size(), &requests[0], &status[0]);
+
+
+  size_t nGloAxisDest = axisDest_->n_glo.getValue() - 1;
+  for (std::map<int,int>::const_iterator itSend = sendRankSizeMap.begin(); itSend != sendRankSizeMap.end(); ++itSend)
+  {
+    int recvRank = itSend->first;
+    int recvSize = itSend->second;
+
+    double* recvValue = recvValueFromSrc[recvRank];
+    unsigned long* recvIndex = sendGlobalIndexSrc[recvRank];
+    for (int idx = 0; idx < recvSize; ++idx)
+    {
+      size_t globalIndex = *(recvIndex+idx);
+      int localIndex = ((nGloAxisDest-globalIndex) - axisDest_->begin);
+      axisDest_->value(localIndex) = *(recvValue + idx);
+    }
+  }
+
+  delete [] recvCount;
+  delete [] displ;
+  delete [] sendRankBuff;
+  delete [] recvRankBuff;
+  delete [] sendSizeBuff;
+  delete [] recvSizeBuff;
+
+  std::unordered_map<int, double* >::const_iterator itChar;
+  for (itChar = sendValueToDest.begin(); itChar != sendValueToDest.end(); ++itChar)
+    delete [] itChar->second;
+  for (itChar = recvValueFromSrc.begin(); itChar != recvValueFromSrc.end(); ++itChar)
+    delete [] itChar->second;
+  std::unordered_map<int, unsigned long* >::const_iterator itLong;
+  for (itLong = sendGlobalIndexSrc.begin(); itLong != sendGlobalIndexSrc.end(); ++itLong)
+    delete [] itLong->second;
+  for (itLong = recvGlobalIndexSrc.begin(); itLong != recvGlobalIndexSrc.end(); ++itLong)
+    delete [] itLong->second;
+}
+CATCH
+
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_inverse.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_inverse.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_inverse.hpp	(revision 1985)
@@ -0,0 +1,53 @@
+/*!
+   \file axis_algorithm_inverse.hpp
+   \author Ha NGUYEN
+   \since 14 May 2015
+   \date 09 June 2015
+
+   \brief Algorithm for inversing an axis..
+ */
+#ifndef __XIOS_AXIS_ALGORITHM_INVERSE_HPP__
+#define __XIOS_AXIS_ALGORITHM_INVERSE_HPP__
+
+#include "axis_algorithm_transformation.hpp"
+#include "transformation.hpp"
+
+namespace xios {
+
+class CAxis;
+class CInverseAxis;
+
+/*!
+  \class CAxisAlgorithmInverse
+  Inversing an axis
+*/
+class CAxisAlgorithmInverse : public CAxisAlgorithmTransformation
+{
+public:
+  CAxisAlgorithmInverse(CAxis* axisDestination, CAxis* axisSource, CInverseAxis* inverseAxis);
+
+  virtual ~CAxisAlgorithmInverse() {}
+
+  static bool registerTrans();
+
+protected:
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs);
+
+private:
+  void updateAxisValue();
+
+private:
+  static CGenericAlgorithmTransformation* create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                CTransformation<CAxis>* transformation,
+                                                int elementPositionInGrid,
+                                                std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridDst2DomainPosition);
+  static bool dummyRegistered_;
+};
+
+}
+#endif // __XIOS_AXIS_ALGORITHM_INVERSE_HPP__
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_reduce_axis.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_reduce_axis.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_reduce_axis.cpp	(revision 1985)
@@ -0,0 +1,127 @@
+/*!
+   \file axis_algorithm_reduce_axis.cpp
+   \author Ha NGUYEN
+   \since 23 June 2016
+   \date 23 June 2016
+
+   \brief Algorithm for reduce a axis to an axis
+ */
+#include "axis_algorithm_reduce_axis.hpp"
+#include "reduce_axis_to_axis.hpp"
+#include "axis.hpp"
+#include "grid.hpp"
+#include "grid_transformation_factory_impl.hpp"
+#include "reduction.hpp"
+
+namespace xios {
+CGenericAlgorithmTransformation* CAxisAlgorithmReduceAxis::create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                                   CTransformation<CAxis>* transformation,
+                                                                   int elementPositionInGrid,
+                                                                   std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                                   std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                                   std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                                   std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                                   std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                                   std::map<int, int>& elementPositionInGridDst2DomainPosition)
+TRY
+{
+  std::vector<CAxis*> axisListDestP = gridDst->getAxis();
+  std::vector<CAxis*> axisListSrcP  = gridSrc->getAxis();
+
+  CReduceAxisToAxis* reduceAxis = dynamic_cast<CReduceAxisToAxis*> (transformation);
+  int axisDstIndex   = elementPositionInGridDst2AxisPosition[elementPositionInGrid];
+  int axisSrcIndex = elementPositionInGridSrc2AxisPosition[elementPositionInGrid];
+
+  return (new CAxisAlgorithmReduceAxis(axisListDestP[axisDstIndex], axisListSrcP[axisSrcIndex], reduceAxis));
+}
+CATCH
+
+bool CAxisAlgorithmReduceAxis::dummyRegistered_ = CAxisAlgorithmReduceAxis::registerTrans();
+bool CAxisAlgorithmReduceAxis::registerTrans()
+TRY
+{
+  return CGridTransformationFactory<CAxis>::registerTransformation(TRANS_REDUCE_AXIS_TO_AXIS, create);
+}
+CATCH
+
+
+CAxisAlgorithmReduceAxis::CAxisAlgorithmReduceAxis(CAxis* axisDestination, CAxis* axisSource, CReduceAxisToAxis* algo)
+ : CAxisAlgorithmTransformation(axisDestination, axisSource), reduction_(0)
+TRY
+{
+  eliminateRedondantSrc_= false ;
+  algo->checkValid(axisDestination, axisSource);
+  StdString op;
+  switch (algo->operation)
+  {
+    case CReduceAxisToAxis::operation_attr::sum:
+      op = "sum";
+      break;
+    case CReduceAxisToAxis::operation_attr::min:
+      op = "min";
+      break;
+    case CReduceAxisToAxis::operation_attr::max:
+      op = "max";
+      break;
+    case CReduceAxisToAxis::operation_attr::average:
+      op = "average";
+      break;
+    default:
+        ERROR("CAxisAlgorithmReduceAxis::CAxisAlgorithmReduceAxis(CAxis* axisDestination, CAxis* axisSource, CReduceAxisToAxis* algo)",
+         << "Operation is wrongly defined. Supported operations: sum, min, max, average." << std::endl
+         << "Axis source " <<axisSource->getId() << std::endl
+         << "Axis destination " << axisDestination->getId());
+
+  }
+
+  reduction_ = CReductionAlgorithm::createOperation(CReductionAlgorithm::ReductionOperations[op]);
+}
+CATCH
+
+void CAxisAlgorithmReduceAxis::apply(const std::vector<std::pair<int,double> >& localIndex,
+                                       const double* dataInput,
+                                       CArray<double,1>& dataOut,
+                                       std::vector<bool>& flagInitial,                     
+                                       bool ignoreMissingValue, bool firstPass)
+TRY
+{
+  reduction_->apply(localIndex, dataInput, dataOut, flagInitial, ignoreMissingValue, firstPass);
+}
+CATCH
+
+void CAxisAlgorithmReduceAxis::updateData(CArray<double,1>& dataOut)
+TRY
+{
+  reduction_->updateData(dataOut);
+}
+CATCH
+
+CAxisAlgorithmReduceAxis::~CAxisAlgorithmReduceAxis()
+TRY
+{
+  if (0 != reduction_) delete reduction_;
+}
+CATCH
+
+void CAxisAlgorithmReduceAxis::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+TRY
+{
+  this->transformationMapping_.resize(1);
+  this->transformationWeight_.resize(1);
+
+  TransformationIndexMap& transMap = this->transformationMapping_[0];
+  TransformationWeightMap& transWeight = this->transformationWeight_[0];
+  CArray<int,1>& axisDstIndex = axisDest_->index;
+  int nbAxisIdx = axisDstIndex.numElements();
+  for (int idxAxis = 0; idxAxis < nbAxisIdx; ++idxAxis)
+  {
+    int globalAxisIdx = axisDstIndex(idxAxis);
+    transMap[globalAxisIdx].resize(1);
+    transWeight[globalAxisIdx].resize(1);
+    transMap[globalAxisIdx][0]=globalAxisIdx ;      
+    transWeight[globalAxisIdx][0] = 1.0 ;      
+  }
+}
+CATCH
+
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_reduce_axis.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_reduce_axis.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_reduce_axis.hpp	(revision 1985)
@@ -0,0 +1,63 @@
+/*!
+   \file axis_algorithm_reduce_axis.hpp
+   \author Ha NGUYEN
+   \since 23 June 2016
+   \date 23 June 2016
+
+   \brief Algorithm for reduce an axis to a axis
+ */
+#ifndef __XIOS_AXIS_ALGORITHM_REDUCE_AXIS_HPP__
+#define __XIOS_AXIS_ALGORITHM_REDUCE_AXIS_HPP__
+
+#include "axis_algorithm_transformation.hpp"
+#include "transformation.hpp"
+
+namespace xios {
+
+class CAxis;
+class CReduceAxisToAxis;
+class CReductionAlgorithm;
+
+/*!
+  \class CAxisAlgorithmReduceAxis
+  Reduce a axis to an axis
+*/
+class CAxisAlgorithmReduceAxis : public CAxisAlgorithmTransformation
+{
+public:
+  CAxisAlgorithmReduceAxis(CAxis* axisDestination, CAxis* axisSource, CReduceAxisToAxis* algo);
+
+  virtual void apply(const std::vector<std::pair<int,double> >& localIndex,
+                     const double* dataInput,
+                     CArray<double,1>& dataOut,
+                     std::vector<bool>& flagInitial,                     
+                     bool ignoreMissingValue, bool firstPass);
+
+  virtual void updateData(CArray<double,1>& dataOut);
+  
+  virtual ~CAxisAlgorithmReduceAxis();
+
+  static bool registerTrans();
+protected:
+  
+protected:
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs); 
+  
+protected:
+  CReductionAlgorithm* reduction_;
+
+private:
+  static CGenericAlgorithmTransformation* create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                CTransformation<CAxis>* transformation,
+                                                int elementPositionInGrid,
+                                                std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridDst2DomainPosition);
+  static bool dummyRegistered_;
+};
+
+}
+#endif // __XIOS_AXIS_ALGORITHM_REDUCE_AXIS_HPP__
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_reduce_domain.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_reduce_domain.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_reduce_domain.cpp	(revision 1985)
@@ -0,0 +1,192 @@
+/*!
+   \file axis_algorithm_reduce_domain.cpp
+   \author Ha NGUYEN
+   \since 23 June 2016
+   \date 23 June 2016
+
+   \brief Algorithm for reduce a domain to an axis
+ */
+#include "axis_algorithm_reduce_domain.hpp"
+#include "reduce_domain_to_axis.hpp"
+#include "axis.hpp"
+#include "domain.hpp"
+#include "grid.hpp"
+#include "grid_transformation_factory_impl.hpp"
+#include "reduction.hpp"
+
+namespace xios {
+CGenericAlgorithmTransformation* CAxisAlgorithmReduceDomain::create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                                   CTransformation<CAxis>* transformation,
+                                                                   int elementPositionInGrid,
+                                                                   std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                                   std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                                   std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                                   std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                                   std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                                   std::map<int, int>& elementPositionInGridDst2DomainPosition)
+TRY
+{
+  std::vector<CAxis*> axisListDestP = gridDst->getAxis();
+  std::vector<CDomain*> domainListSrcP = gridSrc->getDomains();
+
+  CReduceDomainToAxis* reduceDomain = dynamic_cast<CReduceDomainToAxis*> (transformation);
+  int axisDstIndex   = elementPositionInGridDst2AxisPosition[elementPositionInGrid];
+  int domainSrcIndex = elementPositionInGridSrc2DomainPosition[elementPositionInGrid];
+
+  return (new CAxisAlgorithmReduceDomain(axisListDestP[axisDstIndex], domainListSrcP[domainSrcIndex], reduceDomain));
+}
+CATCH
+
+bool CAxisAlgorithmReduceDomain::dummyRegistered_ = CAxisAlgorithmReduceDomain::registerTrans();
+bool CAxisAlgorithmReduceDomain::registerTrans()
+TRY
+{
+  return CGridTransformationFactory<CAxis>::registerTransformation(TRANS_REDUCE_DOMAIN_TO_AXIS, create);
+}
+CATCH
+
+
+CAxisAlgorithmReduceDomain::CAxisAlgorithmReduceDomain(CAxis* axisDestination, CDomain* domainSource, CReduceDomainToAxis* algo)
+ : CAxisAlgorithmTransformation(axisDestination, domainSource), reduction_(0)
+TRY
+{
+  algo->checkValid(axisDestination, domainSource);
+  StdString op;
+  switch (algo->operation)
+  {
+    case CReduceDomainToAxis::operation_attr::sum:
+      op = "sum";
+      break;
+    case CReduceDomainToAxis::operation_attr::min:
+      op = "min";
+      break;
+    case CReduceDomainToAxis::operation_attr::max:
+      op = "max";
+      break;
+    case CReduceDomainToAxis::operation_attr::average:
+      op = "average";
+      break;
+    default:
+        ERROR("CAxisAlgorithmReduceDomain::CAxisAlgorithmReduceDomain(CAxis* axisDestination, CDomain* domainSource, CReduceDomainToAxis* algo)",
+         << "Operation is wrongly defined. Supported operations: sum, min, max, average." << std::endl
+         << "Domain source " <<domainSource->getId() << std::endl
+         << "Axis destination " << axisDestination->getId());
+
+  }
+
+  dir_ = (CReduceDomainToAxis::direction_attr::iDir == algo->direction)  ? iDir : jDir;
+  reduction_ = CReductionAlgorithm::createOperation(CReductionAlgorithm::ReductionOperations[op]);
+  local = algo->local ;
+}
+CATCH
+
+void CAxisAlgorithmReduceDomain::apply(const std::vector<std::pair<int,double> >& localIndex,
+                                       const double* dataInput,
+                                       CArray<double,1>& dataOut,
+                                       std::vector<bool>& flagInitial,                     
+                                       bool ignoreMissingValue, bool firstPass)
+TRY
+{
+  reduction_->apply(localIndex, dataInput, dataOut, flagInitial, ignoreMissingValue, firstPass);
+}
+CATCH
+
+void CAxisAlgorithmReduceDomain::updateData(CArray<double,1>& dataOut)
+TRY
+{
+  reduction_->updateData(dataOut);
+}
+CATCH
+
+CAxisAlgorithmReduceDomain::~CAxisAlgorithmReduceDomain()
+TRY
+{
+  if (0 != reduction_) delete reduction_;
+}
+CATCH
+
+void CAxisAlgorithmReduceDomain::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+TRY
+{
+  this->transformationMapping_.resize(1);
+  this->transformationWeight_.resize(1);
+
+  TransformationIndexMap& transMap = this->transformationMapping_[0];
+  TransformationWeightMap& transWeight = this->transformationWeight_[0];
+
+  CArray<int,1>& axisDstIndex = axisDest_->index;
+  int ni_glo = domainSrc_->ni_glo, nj_glo = domainSrc_->nj_glo;
+  if (iDir == dir_)
+  {
+    if (local)
+    {
+      const CArray<int, 1>& i_index = domainSrc_-> i_index.getValue() ;
+      const CArray<int, 1>& j_index = domainSrc_-> j_index.getValue() ;
+      const CArray<bool,1>& localMask = domainSrc_-> localMask ;
+      int nbDomainIdx = i_index.numElements();
+      
+      for (int idxDomain = 0; idxDomain < nbDomainIdx; ++idxDomain)
+      {
+        if (localMask(idxDomain))
+        { 
+          transMap[j_index(idxDomain)].push_back(j_index(idxDomain)* ni_glo + i_index(idxDomain));
+          transWeight[j_index(idxDomain)].push_back(1.0) ;
+        }
+      }
+    }
+    else
+    {
+      int nbAxisIdx = axisDstIndex.numElements();
+      for (int idxAxis = 0; idxAxis < nbAxisIdx; ++idxAxis)
+      {
+        int globalAxisIdx = axisDstIndex(idxAxis);
+        transMap[globalAxisIdx].resize(ni_glo);
+        transWeight[globalAxisIdx].resize(ni_glo);
+        for (int idx = 0; idx < ni_glo; ++idx)
+        {
+          transMap[globalAxisIdx][idx] = globalAxisIdx * ni_glo + idx;
+          transWeight[globalAxisIdx][idx] = 1.0;
+        }
+      }
+    }
+  }
+  else if (jDir == dir_)
+  {
+    int nbAxisIdx = axisDstIndex.numElements();
+    if (local)
+    {
+      const CArray<int, 1>& i_index = domainSrc_-> i_index.getValue() ;
+      const CArray<int, 1>& j_index = domainSrc_-> j_index.getValue() ;
+      const CArray<bool,1>& localMask = domainSrc_-> localMask ;
+      int nbDomainIdx = i_index.numElements();
+      
+      for (int idxDomain = 0; idxDomain < nbDomainIdx; ++idxDomain)
+      {
+        if (localMask(idxDomain))
+        { 
+          transMap[i_index(idxDomain)].push_back(j_index(idxDomain)* ni_glo + i_index(idxDomain));
+          transWeight[i_index(idxDomain)].push_back(1.0) ;
+        }
+      }
+    }
+    else
+    {
+      for (int idxAxis = 0; idxAxis < nbAxisIdx; ++idxAxis)
+      {
+        int globalAxisIdx = axisDstIndex(idxAxis);
+        transMap[globalAxisIdx].resize(nj_glo);
+        transWeight[globalAxisIdx].resize(nj_glo);
+        for (int idx = 0; idx < nj_glo; ++idx)
+        {
+          transMap[globalAxisIdx][idx] = globalAxisIdx + ni_glo*idx;
+          transWeight[globalAxisIdx][idx] = 1.0;
+        }
+      }
+    }
+  }
+  else
+  {}
+}
+CATCH
+
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_reduce_domain.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_reduce_domain.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_reduce_domain.hpp	(revision 1985)
@@ -0,0 +1,72 @@
+/*!
+   \file axis_algorithm_reduce_domain.hpp
+   \author Ha NGUYEN
+   \since 23 June 2016
+   \date 23 June 2016
+
+   \brief Algorithm for reduce an axis to a axis
+ */
+#ifndef __XIOS_AXIS_ALGORITHM_REDUCE_DOMAIN_HPP__
+#define __XIOS_AXIS_ALGORITHM_REDUCE_DOMAIN_HPP__
+
+#include "axis_algorithm_transformation.hpp"
+#include "transformation.hpp"
+
+namespace xios {
+
+class CAxis;
+class CDomain;
+class CReduceDomainToAxis;
+class CReductionAlgorithm;
+
+/*!
+  \class CAxisAlgorithmReduceDomain
+  Reduce a domain to an axis
+*/
+class CAxisAlgorithmReduceDomain : public CAxisAlgorithmTransformation
+{
+public:
+  CAxisAlgorithmReduceDomain(CAxis* axisDestination, CDomain* domainSource, CReduceDomainToAxis* algo);
+
+  virtual void apply(const std::vector<std::pair<int,double> >& localIndex,
+                     const double* dataInput,
+                     CArray<double,1>& dataOut,
+                     std::vector<bool>& flagInitial,                     
+                     bool ignoreMissingValue, bool firstPass);
+
+  virtual void updateData(CArray<double,1>& dataOut);
+  
+  virtual ~CAxisAlgorithmReduceDomain();
+
+  static bool registerTrans();
+protected:
+  enum ReduceDirection {
+    undefined = 0,
+    iDir = 1,
+    jDir = 2
+  };
+  
+  ReduceDirection dir_;
+  bool local ;
+  
+protected:
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs); 
+  
+protected:
+  CReductionAlgorithm* reduction_;
+
+private:
+  static CGenericAlgorithmTransformation* create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                CTransformation<CAxis>* transformation,
+                                                int elementPositionInGrid,
+                                                std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridDst2DomainPosition);
+  static bool dummyRegistered_;
+};
+
+}
+#endif // __XIOS_AXIS_ALGORITHM_REDUCE_DOMAIN_HPP__
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_temporal_splitting.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_temporal_splitting.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_temporal_splitting.cpp	(revision 1985)
@@ -0,0 +1,81 @@
+/*!
+   \file axis_algorithm_temporal_splitting.cpp
+
+   \brief Algorithm to split scalar into  axis by temporal accumulation
+ */
+#include "axis_algorithm_temporal_splitting.hpp"
+#include "temporal_splitting.hpp"
+#include "axis.hpp"
+#include "scalar.hpp"
+#include "grid.hpp"
+#include "grid_transformation_factory_impl.hpp"
+
+namespace xios {
+CGenericAlgorithmTransformation* CAxisAlgorithmTemporalSplitting::create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                                     CTransformation<CAxis>* transformation,
+                                                                     int elementPositionInGrid,
+                                                                     std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                                     std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                                     std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2DomainPosition)
+TRY
+{
+  std::vector<CAxis*> axisListDestP = gridDst->getAxis();
+  std::vector<CScalar*> scalarListSrcP = gridSrc->getScalars();
+
+  CTemporalSplitting* temporalSplitting = dynamic_cast<CTemporalSplitting*> (transformation);
+  int axisDstIndex = elementPositionInGridDst2AxisPosition[elementPositionInGrid];
+  int scalarSrcIndex = elementPositionInGridSrc2ScalarPosition[elementPositionInGrid];
+
+  return (new CAxisAlgorithmTemporalSplitting(axisListDestP[axisDstIndex], scalarListSrcP[scalarSrcIndex], temporalSplitting));
+}
+CATCH
+
+bool CAxisAlgorithmTemporalSplitting::dummyRegistered_ = CAxisAlgorithmTemporalSplitting::registerTrans();
+bool CAxisAlgorithmTemporalSplitting::registerTrans()
+TRY
+{
+  return CGridTransformationFactory<CAxis>::registerTransformation(TRANS_TEMPORAL_SPLITTING, create);
+}
+CATCH
+
+CAxisAlgorithmTemporalSplitting::CAxisAlgorithmTemporalSplitting(CAxis* axisDestination, CScalar* scalarSource, CTemporalSplitting* algo)
+ : CAxisAlgorithmTransformation(axisDestination, scalarSource)
+{
+
+}
+
+
+CAxisAlgorithmTemporalSplitting::~CAxisAlgorithmTemporalSplitting()
+{
+}
+
+void CAxisAlgorithmTemporalSplitting::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+TRY
+{
+  this->transformationMapping_.resize(1);
+  this->transformationWeight_.resize(1);
+
+  TransformationIndexMap& transMap = this->transformationMapping_[0];
+  TransformationWeightMap& transWeight = this->transformationWeight_[0];
+
+  CArray<int,1>& axisDstIndex = axisDest_->index;
+
+  int nbAxisIdx = axisDstIndex.numElements();
+  for (int idxAxis = 0; idxAxis < nbAxisIdx; ++idxAxis)
+  {
+    int globalAxisIdx = axisDstIndex(idxAxis);
+    if (idxAxis==0)
+    {
+      transMap[globalAxisIdx].resize(1);
+      transWeight[globalAxisIdx].resize(1);
+      transMap[globalAxisIdx][0] = 0 ;
+      transWeight[globalAxisIdx][0] = 1.0 ;
+    }
+  }
+}
+CATCH
+
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_temporal_splitting.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_temporal_splitting.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_temporal_splitting.hpp	(revision 1985)
@@ -0,0 +1,53 @@
+/*!
+   \file axis_algorithm_reduce_domain.hpp
+   \author Ha NGUYEN
+   \since 23 June 2016
+   \date 23 June 2016
+
+   \brief Algorithm for reduce an axis to a axis
+ */
+#ifndef __XIOS_AXIS_ALGORITHM_TEMPORAL_SPLITTING_HPP__
+#define __XIOS_AXIS_ALGORITHM_TEMPORAL_SPLITTING_HPP__
+
+#include "axis_algorithm_transformation.hpp"
+#include "transformation.hpp"
+
+namespace xios {
+
+class CAxis;
+class CScalar;
+class CTemporalSplitting;
+
+
+
+/*!
+  \class CAxisAlgorithmExtractDomain
+  Extract a domain to an axis
+*/
+class CAxisAlgorithmTemporalSplitting : public CAxisAlgorithmTransformation
+{
+public:
+  CAxisAlgorithmTemporalSplitting(CAxis* axisDestination, CScalar* scalarSource, CTemporalSplitting* algo);
+
+  virtual ~CAxisAlgorithmTemporalSplitting();
+
+  static bool registerTrans();
+
+protected:
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs);
+
+private:
+  static CGenericAlgorithmTransformation* create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                CTransformation<CAxis>* transformation,
+                                                int elementPositionInGrid,
+                                                std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridDst2DomainPosition);
+  static bool dummyRegistered_;
+};
+
+}
+#endif // __XIOS_AXIS_ALGORITHM_TEMPORAL_SPLITTING_HPP__
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_transformation.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_transformation.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_transformation.cpp	(revision 1985)
@@ -0,0 +1,135 @@
+/*!
+   \file axis_algorithm_transformation.hpp
+   \author Ha NGUYEN
+   \since 14 May 2015
+   \date 29 June 2015
+
+   \brief Interface for all axis transformation algorithms.
+ */
+
+#include "axis_algorithm_transformation.hpp"
+#include "axis_algorithm_inverse.hpp"
+#include "axis_algorithm_zoom.hpp"
+#include "context.hpp"
+#include "context_client.hpp"
+#include "client_client_dht_template.hpp"
+#include "axis.hpp"
+
+namespace xios {
+
+CAxisAlgorithmTransformation::CAxisAlgorithmTransformation(CAxis* axisDestination, CAxis* axisSource)
+ : CGenericAlgorithmTransformation(), axisDest_(axisDestination), axisSrc_(axisSource), domainSrc_(0),scalarSrc_(0)
+TRY
+{
+  axisDestGlobalSize_ = axisDestination->n_glo.getValue();
+  int niDest = axisDestination->n.getValue();
+  int ibeginDest = axisDestination->begin.getValue();
+
+  for (int idx = 0; idx < niDest; ++idx)
+    if ((axisDestination->mask)(idx)) axisDestGlobalIndex_.push_back(ibeginDest+idx);
+}
+CATCH
+
+CAxisAlgorithmTransformation::CAxisAlgorithmTransformation(CAxis* axisDestination, CDomain* domainSource)
+ : CGenericAlgorithmTransformation(), axisDest_(axisDestination), axisSrc_(0), domainSrc_(domainSource),scalarSrc_(0)
+TRY
+{
+  axisDestGlobalSize_ = axisDestination->n_glo.getValue();
+  int niDest = axisDestination->n.getValue();
+  int ibeginDest = axisDestination->begin.getValue();
+
+  for (int idx = 0; idx < niDest; ++idx)
+    if ((axisDestination->mask)(idx)) axisDestGlobalIndex_.push_back(ibeginDest+idx);
+}
+CATCH
+
+CAxisAlgorithmTransformation::CAxisAlgorithmTransformation(CAxis* axisDestination, CScalar* scalarSource)
+ : CGenericAlgorithmTransformation(), axisDest_(axisDestination), axisSrc_(0), domainSrc_(0), scalarSrc_(scalarSource)
+TRY
+{
+  axisDestGlobalSize_ = axisDestination->n_glo.getValue();
+  int niDest = axisDestination->n.getValue();
+  int ibeginDest = axisDestination->begin.getValue();
+
+  for (int idx = 0; idx < niDest; ++idx)
+    if ((axisDestination->mask)(idx)) axisDestGlobalIndex_.push_back(ibeginDest+idx);
+}
+CATCH
+
+CAxisAlgorithmTransformation::~CAxisAlgorithmTransformation()
+{
+}
+
+void CAxisAlgorithmTransformation::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+{
+}
+
+/*!
+  Compute global index of axis on different processes
+  \param [in] globalAxisIndex global index of axis source
+  \param [out] globalAxisIndexOnProc processes which contain the corresponding global index of axis source
+*/
+void CAxisAlgorithmTransformation::computeExchangeGlobalIndex(const CArray<size_t,1>& globalAxisIndex,
+                                                              int elementType,
+                                                              CClientClientDHTInt::Index2VectorInfoTypeMap& globalAxisIndexOnProc)
+TRY
+{
+  CContext* context = CContext::getCurrent();
+  int clientRank = context->intraCommRank_;
+  int clientSize = context->intraCommSize_;
+
+  size_t globalIndex;
+  int nIndexSize = 0;
+  CArray<bool,1>* ptLocalMask ;
+  CArray<bool,1> scalarMask(1) ;
+  scalarMask=true ;
+  
+  if (2 == elementType) 
+  {
+    nIndexSize = domainSrc_->i_index.numElements();
+    ptLocalMask=&(domainSrc_->localMask) ;
+  }
+  else if (1 == elementType)
+  {
+     nIndexSize = axisSrc_->index.numElements();
+     ptLocalMask=&(axisSrc_->mask) ;
+  }
+  else
+  {
+     nIndexSize=1  ; //  scalar
+     ptLocalMask=&scalarMask ;
+  }
+  CArray<bool,1>& localMask=*ptLocalMask ;
+  
+  CClientClientDHTInt::Index2VectorInfoTypeMap globalIndex2ProcRank;
+  globalIndex2ProcRank.rehash(std::ceil(nIndexSize/globalIndex2ProcRank.max_load_factor()));
+  for (int idx = 0; idx < nIndexSize; ++idx)
+  {
+    if (localMask(idx))
+    {
+      if (2 == elementType) // domain
+      {
+       globalIndex = domainSrc_->i_index(idx) + domainSrc_->j_index(idx) * domainSrc_->ni_glo;
+      }
+      else if (1 == elementType) // axis
+      {
+        globalIndex = axisSrc_->index(idx);
+      }
+      else // scalar
+      {
+        globalIndex = 0;
+      }
+     
+      globalIndex2ProcRank[globalIndex].resize(1);
+      globalIndex2ProcRank[globalIndex][0] = clientRank;
+    }
+  }
+
+  CClientClientDHTInt dhtIndexProcRank(globalIndex2ProcRank, context->intraComm_);
+  dhtIndexProcRank.computeIndexInfoMapping(globalAxisIndex);
+
+  globalAxisIndexOnProc = dhtIndexProcRank.getInfoIndexMap();
+}
+CATCH
+
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_transformation.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_transformation.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_transformation.hpp	(revision 1985)
@@ -0,0 +1,62 @@
+/*!
+   \file axis_algorithm_transformation.hpp
+   \author Ha NGUYEN
+   \since 14 May 2015
+   \date 29 June 2015
+
+   \brief Interface for all axis transformation algorithms.
+ */
+#ifndef __XIOS_AXIS_ALGORITHM_TRANSFORMATION_HPP__
+#define __XIOS_AXIS_ALGORITHM_TRANSFORMATION_HPP__
+
+#include "generic_algorithm_transformation.hpp"
+#include "grid_transformation_factory_impl.hpp"
+
+namespace xios {
+
+class CAxis;
+class CDomain;
+class CScalar;
+/*!
+  \class CAxisAlgorithmTransformation
+  Algorithms for axis.
+*/
+class CAxisAlgorithmTransformation : public virtual CGenericAlgorithmTransformation
+{
+public:
+  CAxisAlgorithmTransformation(CAxis* axisDestination, CAxis* axisSource);
+  CAxisAlgorithmTransformation(CAxis* axisDestination, CDomain* domainSource);
+  CAxisAlgorithmTransformation(CAxis* axisDestination, CScalar* scalarSource);
+
+  virtual ~CAxisAlgorithmTransformation();
+
+protected:
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs);
+
+  virtual void computeExchangeGlobalIndex(const CArray<size_t,1>& globalAxisIndex,
+                                          int elementType,
+                                          CClientClientDHTInt::Index2VectorInfoTypeMap& globalDomainIndexOnProc);
+
+protected:
+  //! Global index of an axis on grid destination
+  std::vector<int> axisDestGlobalIndex_;
+
+  //! Size of
+  int axisDestGlobalSize_;
+
+    //! Axis on grid destination
+  CAxis* axisDest_;
+
+  //! Axis on grid source
+  CAxis* axisSrc_;
+
+  //! Domain on grid source
+  CDomain* domainSrc_;
+
+  //! Scalar on grid source
+
+  CScalar* scalarSrc_;
+};
+
+}
+#endif // __XIOS_AXIS_ALGORITHM_TRANSFORMATION_HPP__
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_zoom.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_zoom.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_zoom.cpp	(revision 1985)
@@ -0,0 +1,122 @@
+/*!
+   \file axis_algorithm_zoom.cpp
+   \brief Algorithm for zooming an axis.
+ */
+#include "axis_algorithm_zoom.hpp"
+#include "axis.hpp"
+#include "grid.hpp"
+#include "grid_transformation_factory_impl.hpp"
+#include "zoom_axis.hpp"
+
+namespace xios {
+CGenericAlgorithmTransformation* CAxisAlgorithmZoom::create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                           CTransformation<CAxis>* transformation,
+                                                           int elementPositionInGrid,
+                                                           std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                           std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                           std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                           std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                           std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                           std::map<int, int>& elementPositionInGridDst2DomainPosition)
+TRY
+{
+  std::vector<CAxis*> axisListDestP = gridDst->getAxis();
+  std::vector<CAxis*> axisListSrcP  = gridSrc->getAxis();
+
+  CZoomAxis* zoomAxis = dynamic_cast<CZoomAxis*> (transformation);
+  int axisDstIndex = elementPositionInGridDst2AxisPosition[elementPositionInGrid];
+  int axisSrcIndex = elementPositionInGridSrc2AxisPosition[elementPositionInGrid];
+
+  return (new CAxisAlgorithmZoom(axisListDestP[axisDstIndex], axisListSrcP[axisSrcIndex], zoomAxis));
+}
+CATCH
+
+bool CAxisAlgorithmZoom::dummyRegistered_ = CAxisAlgorithmZoom::registerTrans();
+bool CAxisAlgorithmZoom::registerTrans()
+TRY
+{
+  return CGridTransformationFactory<CAxis>::registerTransformation(TRANS_ZOOM_AXIS, create);
+}
+CATCH
+
+CAxisAlgorithmZoom::CAxisAlgorithmZoom(CAxis* axisDestination, CAxis* axisSource, CZoomAxis* zoomAxis)
+: CAxisAlgorithmTransformation(axisDestination, axisSource)
+TRY
+{
+  zoomAxis->checkValid(axisSource);
+  zoomBegin_ = zoomAxis->begin.getValue();
+  zoomN_  = zoomAxis->n.getValue();
+  zoomEnd_   = zoomBegin_ + zoomN_ - 1;
+
+  if (zoomN_ > axisSource->n_glo.getValue())
+  {
+    ERROR("CAxisAlgorithmZoom::CAxisAlgorithmZoom(CAxis* axisDestination, CAxis* axisSource, CZoomAxis* zoomAxis)",
+           << "Zoom size is greater than global size of source axis"
+           << "Global size of source axis " <<axisSource->getId() << " is " << axisSource->n_glo.getValue()  << std::endl
+           << "Zoom size is " << zoomN_ );
+  }
+
+  int idxSrc, nDest = 0, beginDestLoc, beginDestGlo = 0 ;
+  int indGloDest, indGloSrc, iSrc;
+  for (int i = 0; i < axisSrc_->n.getValue(); i++)
+  {
+    idxSrc = axisSrc_->index(i);
+    if ((idxSrc >= zoomBegin_) && (idxSrc <= zoomEnd_))
+    {
+      if (nDest == 0) beginDestLoc = i;
+      ++nDest;
+    }
+  }
+  beginDestGlo = beginDestLoc + axisSrc_->begin - zoomBegin_;
+  axisDest_->n_glo.setValue(zoomN_);
+  axisDest_->n.setValue(nDest);
+  axisDest_->begin.setValue(beginDestGlo);
+  axisDest_->index.resize(nDest);
+
+  axisDest_->data_n.setValue(nDest);
+  axisDest_->data_begin.setValue(0);
+  axisDest_->data_index.resize(nDest);
+
+  axisDest_->mask.resize(nDest);
+  if (axisSrc_->hasValue) axisDest_->value.resize(nDest);
+  if (axisSrc_->hasLabel) axisDest_->label.resize(nDest);
+  if (axisSrc_->hasBounds) axisDest_->bounds.resize(2,nDest);
+
+  this->transformationMapping_.resize(1);
+  this->transformationWeight_.resize(1);
+  TransformationIndexMap& transMap = this->transformationMapping_[0];
+  TransformationWeightMap& transWeight = this->transformationWeight_[0];
+
+  for (int iDest = 0; iDest < nDest; iDest++)
+  {
+    iSrc = iDest + beginDestLoc;
+    axisDest_->index(iDest) = iDest + beginDestGlo;
+    axisDest_->data_index(iDest) = axisSrc_->data_index(iSrc) - beginDestLoc;
+    axisDest_->mask(iDest) = axisSrc_->mask(iSrc);
+
+    if (axisSrc_->hasValue)
+      axisDest_->value(iDest) = axisSrc_->value(iSrc);
+    if (axisSrc_->hasLabel)
+      axisDest_->label(iDest) = axisSrc_->label(iSrc);
+    if (axisSrc_->hasBounds)
+    {
+      axisDest_->bounds(0,iDest) = axisSrc_->bounds(0,iSrc);
+      axisDest_->bounds(1,iDest) = axisSrc_->bounds(1,iSrc);
+    }
+    indGloDest = axisDest_->index(iDest);
+    indGloSrc = axisSrc_->index(iSrc);
+    transMap[indGloDest].push_back(indGloSrc);
+    transWeight[indGloDest].push_back(1.0);
+
+  }
+}
+CATCH
+
+/*!
+  Compute the index mapping between domain on grid source and one on grid destination
+*/
+void CAxisAlgorithmZoom::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+{
+}
+
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_zoom.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_zoom.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_zoom.hpp	(revision 1985)
@@ -0,0 +1,61 @@
+/*!
+   \file axis_algorithm_zoom.hpp
+   \brief Algorithm for zooming an axis.
+ */
+#ifndef __XIOS_AXIS_ALGORITHM_ZOOM_HPP__
+#define __XIOS_AXIS_ALGORITHM_ZOOM_HPP__
+
+#include "axis_algorithm_transformation.hpp"
+#include "transformation.hpp"
+
+namespace xios {
+class CAxis;
+class CZoomAxis;
+
+/*!
+  \class CAxisAlgorithmZoom
+  Implementing zoom on axis
+  A zoomed region can be considered as region that isn't masked.
+  Only this zoomed region is zoomed to write on Netcdf.
+*/
+class CAxisAlgorithmZoom : public CAxisAlgorithmTransformation
+{
+public:
+  CAxisAlgorithmZoom(CAxis* axisDestination, CAxis* axisSource, CZoomAxis* zoomAxis);
+
+  virtual ~CAxisAlgorithmZoom() {}
+
+  static bool registerTrans();
+protected:
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs);
+
+private:
+
+private:
+  //! Global zoom begin on axis
+  StdSize zoomBegin_;
+
+  //! Global zoom end on axis
+  StdSize zoomEnd_;
+
+  //! Global zoom size on axis
+  StdSize zoomN_;
+
+  std::vector<int> zoomIndex_;
+
+private:
+
+  static CGenericAlgorithmTransformation* create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                CTransformation<CAxis>* transformation,
+                                                int elementPositionInGrid,
+                                                std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridDst2DomainPosition);
+  static bool dummyRegistered_;
+};
+
+}
+#endif // __XIOS_AXIS_ALGORITHM_ZOOM_HPP__
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_compute_connectivity.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_compute_connectivity.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_compute_connectivity.cpp	(revision 1985)
@@ -0,0 +1,112 @@
+/*!
+   \file domain_algorithm_compute_connectivity.cpp
+   \author Ha NGUYEN
+   \since 15 Jul 2016
+   \date 15 Jul 2016
+
+   \brief Algorithm for compute_connectivity on an domain.
+ */
+#include "domain_algorithm_compute_connectivity.hpp"
+#include "compute_connectivity_domain.hpp"
+#include "mesh.hpp"
+#include "domain.hpp"
+#include "grid.hpp"
+#include "grid_transformation_factory_impl.hpp"
+
+namespace xios {
+CGenericAlgorithmTransformation* CDomainAlgorithmComputeConnectivity::create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                                     CTransformation<CDomain>* transformation,
+                                                                     int elementPositionInGrid,
+                                                                     std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                                     std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                                     std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2DomainPosition)
+TRY
+{
+  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
+  std::vector<CDomain*> domainListSrcP  = gridSrc->getDomains();
+
+  CComputeConnectivityDomain* compute_connectivityDomain = dynamic_cast<CComputeConnectivityDomain*> (transformation);
+  int domainDstIndex = elementPositionInGridDst2DomainPosition[elementPositionInGrid];
+  int domainSrcIndex = elementPositionInGridSrc2DomainPosition[elementPositionInGrid];
+
+  return (new CDomainAlgorithmComputeConnectivity(domainListDestP[domainDstIndex], domainListSrcP[domainSrcIndex], compute_connectivityDomain));
+}
+CATCH
+
+bool CDomainAlgorithmComputeConnectivity::dummyRegistered_ = CDomainAlgorithmComputeConnectivity::registerTrans();
+bool CDomainAlgorithmComputeConnectivity::registerTrans()
+TRY
+{
+  return CGridTransformationFactory<CDomain>::registerTransformation(TRANS_COMPUTE_CONNECTIVITY_DOMAIN, create);
+}
+CATCH
+
+CDomainAlgorithmComputeConnectivity::CDomainAlgorithmComputeConnectivity(CDomain* domainDestination, CDomain* domainSource,
+                                                                         CComputeConnectivityDomain* compute_connectivityDomain)
+: CDomainAlgorithmTransformation(domainDestination, domainSource)
+TRY
+{
+  this->type_ = (ELEMENT_NO_MODIFICATION_WITHOUT_DATA);
+  compute_connectivityDomain->checkValid(domainDestination);
+  int& nbNeighborMax = compute_connectivityDomain->n_neighbor_max;
+  CArray<int,1>& nbNeighbor = compute_connectivityDomain->n_neighbor;
+  CArray<int,2>& localNeighbors = compute_connectivityDomain->local_neighbor;
+  int type = 1; // Edge type
+  switch (compute_connectivityDomain->type)
+  {
+    case CComputeConnectivityDomain::type_attr::node :
+      type = 0;
+      break;
+    case CComputeConnectivityDomain::type_attr::edge :
+      type = 1;
+      break;
+    default:
+      break;
+  }
+
+  computeLocalConnectivity(type, domainDestination, nbNeighborMax, nbNeighbor, localNeighbors);
+}
+CATCH
+
+/*!
+ *  Compute local connectivity of a domain
+ *  \param[in] type type of connectivity (node or edge)
+ *  \param[in] domain domain on which we calculate local connectivity
+ *  \param[in/out] nbConnectivityMax maximum number of neighbor a cell of domain has
+ *  \param[in/out] nbConnectivity number of neighbor a cell has
+ *  \param[in/out] localConnectivity localConnectivity local index of neighbor of a cell
+ */
+void CDomainAlgorithmComputeConnectivity::computeLocalConnectivity(int type,
+                                                                  CDomain* domain,
+                                                                  int& nbConnectivityMax,
+                                                                  CArray<int,1>& nbConnectivity,
+                                                                  CArray<int,2>& localConnectivity)
+TRY
+{
+
+  CMesh mesh;
+
+  CArray<double,2>& bounds_lon = domain->bounds_lon_1d;
+  CArray<double,2>& bounds_lat = domain->bounds_lat_1d;
+  int ncell = bounds_lon.shape()[1];
+  CArray<int,1> localIndex(ncell);
+  for (int idx = 0; idx <ncell; ++idx) localIndex(idx) = idx;
+
+  mesh.getLocalNghbFaces(type, localIndex, bounds_lon, bounds_lat, localConnectivity, nbConnectivity);
+  nbConnectivityMax = 0;
+  for (int idx =0; idx < nbConnectivity.numElements(); ++idx)
+    if (nbConnectivityMax < nbConnectivity(idx)) nbConnectivityMax = nbConnectivity(idx);
+}
+CATCH
+
+/*!
+  Compute the index mapping between domain on grid source and one on grid destination
+*/
+void CDomainAlgorithmComputeConnectivity::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+{
+}
+
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_compute_connectivity.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_compute_connectivity.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_compute_connectivity.hpp	(revision 1985)
@@ -0,0 +1,54 @@
+/*!
+   \file domain_algorithm_compute_connectivity.hpp
+   \author Ha NGUYEN
+   \since 03 June 2015
+   \date 12 June 2015
+
+   \brief Algorithm for compute_connectivitying on an domain.
+ */
+#ifndef __XIOS_DOMAIN_ALGORITHM_COMPUTE_CONNECTIVITY_HPP__
+#define __XIOS_DOMAIN_ALGORITHM_COMPUTE_CONNECTIVITY_HPP__
+
+#include "domain_algorithm_transformation.hpp"
+#include "transformation.hpp"
+
+namespace xios {
+class CDomain;
+class CComputeConnectivityDomain;
+
+/*!
+  \class CDomainAlgorithmComputeConnectivity
+*/
+class CDomainAlgorithmComputeConnectivity : public CDomainAlgorithmTransformation
+{
+public:
+  CDomainAlgorithmComputeConnectivity(CDomain* domainDestination, CDomain* domainSource, CComputeConnectivityDomain* compute_connectivityDomain);
+
+  virtual ~CDomainAlgorithmComputeConnectivity() {}
+
+  static bool registerTrans();
+protected:
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs);
+
+protected:
+  void computeLocalConnectivity(int type,
+                                CDomain* domain,
+                                int& nbConnectivityMax,
+                                CArray<int,1>& nbConnectivity,
+                                CArray<int,2>& localConnectivity);
+
+private:
+  static CGenericAlgorithmTransformation* create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                CTransformation<CDomain>* transformation,
+                                                int elementPositionInGrid,
+                                                std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridDst2DomainPosition);
+  static bool dummyRegistered_;
+};
+
+}
+#endif // __XIOS_DOMAIN_ALGORITHM_COMPUTE_CONNECTIVITY_HPP__
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_expand.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_expand.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_expand.cpp	(revision 1985)
@@ -0,0 +1,714 @@
+/*!
+   \file domain_algorithm_expand.cpp
+   \author Ha NGUYEN
+   \since 08 Aug 2016
+   \date 19 Sep 2016
+
+   \brief Algorithm for expanding an domain.
+ */
+#include "domain_algorithm_expand.hpp"
+#include "expand_domain.hpp"
+#include "mesh.hpp"
+#include "domain.hpp"
+#include "grid.hpp"
+#include "grid_transformation_factory_impl.hpp"
+#include "context.hpp"
+#include "context_client.hpp"
+
+namespace xios {
+CGenericAlgorithmTransformation* CDomainAlgorithmExpand::create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                               CTransformation<CDomain>* transformation,
+                                                               int elementPositionInGrid,
+                                                               std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                               std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                               std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                               std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                               std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                               std::map<int, int>& elementPositionInGridDst2DomainPosition)
+TRY
+{
+  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
+  std::vector<CDomain*> domainListSrcP  = gridSrc->getDomains();
+
+  CExpandDomain* expandDomain = dynamic_cast<CExpandDomain*> (transformation);
+  int domainDstIndex = elementPositionInGridDst2DomainPosition[elementPositionInGrid];
+  int domainSrcIndex = elementPositionInGridSrc2DomainPosition[elementPositionInGrid];
+
+  return (new CDomainAlgorithmExpand(domainListDestP[domainDstIndex], domainListSrcP[domainSrcIndex], expandDomain));
+}
+CATCH
+
+bool CDomainAlgorithmExpand::dummyRegistered_ = CDomainAlgorithmExpand::registerTrans();
+bool CDomainAlgorithmExpand::registerTrans()
+TRY
+{
+  return CGridTransformationFactory<CDomain>::registerTransformation(TRANS_EXPAND_DOMAIN, create);
+}
+CATCH
+
+CDomainAlgorithmExpand::CDomainAlgorithmExpand(CDomain* domainDestination,
+                                               CDomain* domainSource,
+                                               CExpandDomain* expandDomain)
+: CDomainAlgorithmTransformation(domainDestination, domainSource),
+  isXPeriodic_(false), isYPeriodic_(false)
+TRY
+{
+  if (domainDestination == domainSource)
+  {
+    ERROR("CDomainAlgorithmExpand::CDomainAlgorithmExpand(CDomain* domainDestination,CDomain* domainSource, CExpandDomain* expandDomain)",
+           << "Domain source and domain destination are the same. Please make sure domain destination refers to domain source" << std::endl
+           << "Domain source " <<domainSource->getId() << std::endl
+           << "Domain destination " <<domainDestination->getId() << std::endl);
+  }
+
+  this->type_ = (ELEMENT_MODIFICATION_WITH_DATA);
+  // Make sure domain source have all valid attributes
+  // domainSource->checkAllAttributes();
+  expandDomain->checkValid(domainDestination);
+  if (!expandDomain->i_periodic.isEmpty()) isXPeriodic_ = expandDomain->i_periodic;
+  if (!expandDomain->j_periodic.isEmpty()) isYPeriodic_ = expandDomain->j_periodic;
+  
+  switch (expandDomain->type)
+  {
+    case CExpandDomain::type_attr::node :
+      expandDomainNodeConnectivity(domainDestination,
+                                   domainSource);
+      break;
+    case CExpandDomain::type_attr::edge :
+      expandDomainEdgeConnectivity(domainDestination,
+                                   domainSource);
+      break;
+    default:
+      break;
+  }
+}
+CATCH
+
+/*!
+ *  Expand domain with edge-type neighbor
+ *  \param[in/out] domainDestination domain destination and will be modified
+ *  \param[in] domainSource domain source
+ */
+void CDomainAlgorithmExpand::expandDomainEdgeConnectivity(CDomain* domainDestination,
+                                                          CDomain* domainSource)
+TRY
+{
+  CContext* context = CContext::getCurrent();
+ 
+  int type = 1; // For edge
+  CMesh mesh;
+  CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d;
+  CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d;
+  CArray<int,2> neighborsSrc;
+  switch (domainSource->type) {
+   case CDomain::type_attr::unstructured:      
+      mesh.getGlobalNghbFaces(type, context->intraComm_, domainSource->i_index, bounds_lon_src, bounds_lat_src, neighborsSrc);
+      updateUnstructuredDomainAttributes(domainDestination, domainSource, neighborsSrc);
+      break;
+   default:
+      updateRectilinearDomainAttributes(domainDestination, domainSource, neighborsSrc);
+      break;      
+  }  
+}
+CATCH
+
+/*!
+ *  Expand domain with node-type neighbor
+ *  \param[in/out] domainDestination domain destination and will be modified
+ *  \param[in] domainSource domain source
+ */
+void CDomainAlgorithmExpand::expandDomainNodeConnectivity(CDomain* domainDestination,
+                                                          CDomain* domainSource)
+TRY
+{
+  CContext* context = CContext::getCurrent();
+
+  int type = 1; // For edge
+  CMesh mesh;
+  CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d;
+  CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d;
+  CArray<int,2> neighborsSrc;
+  switch (domainSource->type) {
+   case CDomain::type_attr::unstructured:      
+      mesh.getGlobalNghbFaces(type, context->intraComm_, domainSource->i_index, bounds_lon_src, bounds_lat_src, neighborsSrc);
+      updateUnstructuredDomainAttributes(domainDestination, domainSource, neighborsSrc);
+      break;
+   default:
+      updateRectilinearDomainAttributes(domainDestination, domainSource, neighborsSrc);
+      break;      
+  }  
+}
+CATCH
+
+/*!
+ *  Extend rectilinear or curvilinear domain destination and update its attributes
+ *  Suppose that domain destination and domain source have the same values for all attributes (by inheritance)
+ *  \param [in/out] domainDestination domain destination
+ *  \param [in] domainSource domain source
+ *  \param [in] neighborsDomainSrc neighbor of domain source. For now, we don't need it for rectilinear
+ */
+void CDomainAlgorithmExpand::updateRectilinearDomainAttributes(CDomain* domainDestination,
+                                                               CDomain* domainSource,
+                                                               CArray<int,2>& neighborsDomainSrc)
+TRY
+{
+  int index, globalIndex, idx;
+  int iindexDst, jindexDst, globIndexDst;
+  int iindexSrc, jindexSrc, globIndexSrc;
+  CContext* context = CContext::getCurrent();
+
+  // First of all, "copy" all attributes of domain source to domain destination
+  StdString domainDstRef = (!domainDestination->domain_ref.isEmpty()) ? domainDestination->domain_ref.getValue()
+                                                                      : "";
+  if (domainDstRef != domainSource->getId())
+  {
+    domainDestination->domain_ref.setValue(domainSource->getId());
+    domainDestination->solveRefInheritance(true);
+  }
+
+  if (domainDstRef.empty()) domainDestination->domain_ref.reset();
+  else domainDestination->domain_ref.setValue(domainDstRef);
+
+  
+  // Here are attributes of source need tranfering
+  int niGloSrc = domainSource->ni_glo;
+  int njGloSrc = domainSource->nj_glo;
+  int niSrc = domainSource->ni, ibegin = domainSource->ibegin;
+  int njSrc = domainSource->nj, jbegin = domainSource->jbegin;
+  int dataDimSrc = domainSource->data_dim;
+  CArray<bool,1>& mask_1d_src = domainSource->domainMask;
+  CArray<int,1>& i_index_src = domainSource->i_index;
+  CArray<int,1>& j_index_src = domainSource->j_index;
+  CArray<int,1>& data_i_index_src = domainSource->data_i_index;
+  CArray<int,1>& data_j_index_src = domainSource->data_j_index;
+  int data_i_begin_src = domainSource->data_ibegin;
+  int data_j_begin_src = domainSource->data_jbegin;
+  CArray<double,1>& lon_src = domainSource->lonvalue;
+  CArray<double,1>& lat_src = domainSource->latvalue;
+
+  // We need to generate boundary for longitude and latitude
+  if (domainSource->bounds_lon_1d.isEmpty() || domainSource->bounds_lat_1d.isEmpty())
+  {
+    CArray<double,1> lon = lon_src(Range(0,niSrc-1));
+    CArray<double,1> lat = lat_src(Range(0,lat_src.numElements()-niSrc,niSrc));
+    CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d;
+    CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d;
+    domainSource->fillInRectilinearBoundLonLat(lon_src, lat_src, bounds_lon_src, bounds_lat_src);
+  }
+
+
+  CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d;
+  CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d;
+
+  int nVertex       = bounds_lon_src.shape()[0];
+  int oldNbLocal = i_index_src.numElements();
+  // Calculate ni, nj by using i_index and j_index
+  int niSrcByIndex = max(i_index_src) - min(i_index_src) + 1;
+  int njSrcByIndex = max(j_index_src) - min(j_index_src) + 1;  
+  int dataIindexBoundSrc = (1 == dataDimSrc) ? (niSrcByIndex * njSrcByIndex) : niSrcByIndex;
+  int dataJindexBoundSrc = (1 == dataDimSrc) ? (niSrcByIndex * njSrcByIndex) : njSrcByIndex;
+
+  // Uncompress data_i_index, data_j_index
+  CArray<int,1> data_i_index_src_full(oldNbLocal);
+  CArray<int,1> data_j_index_src_full(oldNbLocal);
+  int nbUnMaskedPointOnLocalDomain = 0;
+  data_i_index_src_full = -1; // Suppose all values are masked
+  data_j_index_src_full = -1; // Suppose all values are masked
+  for (idx = 0; idx < data_i_index_src.numElements(); ++idx)
+  {
+    int dataIidx = data_i_index_src(idx) + data_i_begin_src;
+    int dataJidx = data_j_index_src(idx) + data_j_begin_src;
+    if ((0 <= dataIidx) && (dataIidx < dataIindexBoundSrc) &&
+        (0 <= dataJidx) && (dataJidx < dataJindexBoundSrc))
+    {
+      data_i_index_src_full(nbUnMaskedPointOnLocalDomain) = dataIidx;
+      data_j_index_src_full(nbUnMaskedPointOnLocalDomain) = dataJidx;
+      ++nbUnMaskedPointOnLocalDomain;
+    }
+  }
+
+  // Expand domain destination, not only local but also global
+  int niGloDst = niGloSrc + 2;
+  int njGloDst = njGloSrc + 2;
+  int niDst = niSrc + 2;
+  int njDst = njSrc + 2;
+  domainDestination->ni_glo.setValue(niGloDst);
+  domainDestination->nj_glo.setValue(njGloDst);
+  domainDestination->ni.setValue(niDst);
+  domainDestination->nj.setValue(njDst);
+
+  CArray<bool,1>& mask_1d_dst = domainDestination->domainMask;
+  CArray<int,1>& i_index_dst  = domainDestination->i_index;
+  CArray<int,1>& j_index_dst  = domainDestination->j_index;  
+  CArray<int,1>& data_i_index_dst  = domainDestination->data_i_index;
+  CArray<int,1>& data_j_index_dst  = domainDestination->data_j_index;
+  
+  // Make sure that we use only lonvalue_client, latvalue_client
+  if (!domainDestination->lonvalue_1d.isEmpty()) domainDestination->lonvalue_1d.reset();
+  if (!domainDestination->latvalue_1d.isEmpty()) domainDestination->latvalue_1d.reset();
+  if (!domainDestination->lonvalue_2d.isEmpty()) domainDestination->lonvalue_2d.reset();
+  if (!domainDestination->latvalue_2d.isEmpty()) domainDestination->latvalue_2d.reset();
+
+  // Recalculate i_index, j_index of extended domain
+  // Should be enough for common case, but if we have arbitrary distribution?
+  int newNbLocalDst = niDst * njDst;     
+
+  mask_1d_dst.resize(newNbLocalDst);
+  i_index_dst.resize(newNbLocalDst);
+  j_index_dst.resize(newNbLocalDst);
+  CArray<int,1> data_i_index_dst_full(newNbLocalDst);
+  CArray<int,1> data_j_index_dst_full(newNbLocalDst);
+
+  if (newNbLocalDst==0)
+  {
+    domainDestination->lonvalue.resize(newNbLocalDst);
+    domainDestination->latvalue.resize(newNbLocalDst);
+    domainDestination->bounds_lon_1d.resize(nVertex, newNbLocalDst);
+    domainDestination->bounds_lat_1d.resize(nVertex, newNbLocalDst);
+  }
+  else 
+  {
+    domainDestination->lonvalue.resizeAndPreserve(newNbLocalDst);
+    domainDestination->latvalue.resizeAndPreserve(newNbLocalDst);
+    domainDestination->bounds_lon_1d.resizeAndPreserve(nVertex, newNbLocalDst);
+    domainDestination->bounds_lat_1d.resizeAndPreserve(nVertex, newNbLocalDst);
+  }
+  CArray<double,1>& lon_dst   = domainDestination->lonvalue;
+  CArray<double,1>& lat_dst   = domainDestination->latvalue;
+  CArray<double,2>& bounds_lon_dst = domainDestination->bounds_lon_1d;
+  CArray<double,2>& bounds_lat_dst = domainDestination->bounds_lat_1d;
+
+  // Update i_index, j_index
+  for (int j = 0; j < njDst; ++j)
+    for (int i = 0; i < niDst; ++i)
+    {
+      idx = j * niDst + i; 
+      i_index_dst(idx) = i + ibegin;
+      j_index_dst(idx) = j + jbegin;
+    }
+ 
+
+  // 1. Fill in array relating to global index (i_index, j_index, transmap, etc, ...)
+  // Global index mapping between destination and source
+  this->transformationMapping_.resize(1);
+  this->transformationWeight_.resize(1);
+  TransformationIndexMap& transMap = this->transformationMapping_[0];
+  TransformationWeightMap& transWeight = this->transformationWeight_[0];
+
+  transMap.rehash(std::ceil(newNbLocalDst/transMap.max_load_factor()));
+  transWeight.rehash(std::ceil(newNbLocalDst/transWeight.max_load_factor()));
+  
+  // Index mapping for local domain
+  // Mapping global index of expanded domain into original one 
+  // (Representing global index of expanded domain in form of global index of original one)
+  CArray<size_t,1> globalIndexSrcOnDstDomain(newNbLocalDst); 
+  for (idx = 0; idx < newNbLocalDst; ++idx)
+  {
+    iindexDst = i_index_dst(idx);
+    jindexDst = j_index_dst(idx);
+    globIndexDst = jindexDst * niGloDst + iindexDst;
+    globIndexSrc = (((jindexDst-1)+njGloSrc) % njGloSrc) * niGloSrc + (((iindexDst-1)+niGloSrc) % niGloSrc) ;
+    globalIndexSrcOnDstDomain(idx) = globIndexSrc;
+
+    transMap[globIndexDst].push_back(globIndexSrc);
+    transWeight[globIndexDst].push_back(1.0); 
+  }
+
+  // 2. Exchange local info among domains (lon,lat,bounds,mask,etc,...)
+  CClientClientDHTDouble::Index2VectorInfoTypeMap localData;
+  localData.rehash(std::ceil(oldNbLocal/localData.max_load_factor()));
+
+  // Information exchanged among domains (attention to their order), number in parentheses presents size of data
+  // lon(1) + lat(1) + bounds_lon(nVertex) + bounds_lat(nVertex) + mask(1) + data_i_index(1)
+  int dataPackageSize =  1 + 1 + // lon + lat
+                         nVertex + nVertex + //bounds_lon + bounds_lat
+                         1 + // mask_1d_dst;
+                         1 + 1; // data_i_index + data_j_index
+  // Initialize database
+  for (int idx = 0; idx < oldNbLocal; ++idx)
+  {
+    index = i_index_src(idx) + j_index_src(idx) * niGloSrc;
+    localData[index].resize(dataPackageSize);
+    std::vector<double>& data = localData[index];
+
+    //Pack data
+    int dataIdx = 0;
+    data[dataIdx] = lon_src(idx);++dataIdx;
+    data[dataIdx] = lat_src(idx);++dataIdx;
+    for (int i = 0; i < nVertex; ++i)
+    {
+      data[dataIdx] = bounds_lon_src(i,idx); ++dataIdx;
+    }
+    for (int i = 0; i < nVertex; ++i)
+    {
+      data[dataIdx] = bounds_lat_src(i,idx); ++dataIdx;
+    }
+    data[dataIdx] = mask_1d_src(idx) ? 1.0 : -1.0; ++dataIdx;
+    data[dataIdx] = data_i_index_src_full(idx);++dataIdx;
+    data[dataIdx] = data_j_index_src_full(idx);
+  }
+
+  CClientClientDHTDouble dhtData(localData,context->intraComm_);
+  dhtData.computeIndexInfoMapping(globalIndexSrcOnDstDomain);
+  CClientClientDHTDouble::Index2VectorInfoTypeMap& neighborData = dhtData.getInfoIndexMap();
+  CClientClientDHTDouble::Index2VectorInfoTypeMap::iterator ite = neighborData.end(), it;
+
+  // Ok get all data for destination
+  // If domain is not periodic, then we mask all extended part.
+  int nbUnMaskedPointOnExtendedPart = 0, remainder = 0, dataIIndex, dataJIndex;
+  size_t nIdx;
+  double maskValue = 1.0;
+  for (index = 0; index < newNbLocalDst; ++index)
+  {
+     nIdx = globalIndexSrcOnDstDomain(index);
+     it = neighborData.find(nIdx);
+     if (ite != it)
+     {
+        std::vector<double>& data = it->second;
+        // Unpack data
+        int dataIdx = 0;
+        lon_dst(index) = data[dataIdx]; ++dataIdx;
+        lat_dst(index) = data[dataIdx]; ++dataIdx;
+        for (int i = 0; i < nVertex; ++i)
+        {
+          bounds_lon_dst(i,index) = data[dataIdx]; ++dataIdx;
+        }
+        for (int i = 0; i < nVertex; ++i)
+        {
+          bounds_lat_dst(i,index) = data[dataIdx]; ++dataIdx;
+        }
+        
+        // Check whether we have x periodic. If we don't, we should mask all point at 0 and niGloDst-1
+        maskValue = data[dataIdx];
+        if (!isXPeriodic_) 
+        {
+          remainder = i_index_dst(index) % (niGloDst-1);
+          if (0 == remainder) 
+          {
+            maskValue = -1.0;
+          }
+        }
+
+        if (!isYPeriodic_) 
+        {
+          remainder = j_index_dst(index) % (njGloDst-1);
+          if (0 == remainder) 
+          {
+            maskValue = -1.0;
+          }
+        }
+
+        mask_1d_dst(index) = (1.0 == maskValue) ? true : false; ++dataIdx;
+
+        dataIIndex = (int) data[dataIdx];
+        if (!isXPeriodic_) 
+        {
+          remainder = i_index_dst(index) % (niGloDst-1);
+          if (0 == remainder) 
+          {
+            dataIIndex = -1;
+          }
+        }
+        data_i_index_dst_full(index) = dataIIndex; ++dataIdx;
+        
+        dataJIndex = (int) data[dataIdx];
+        if (!isYPeriodic_) 
+        {
+          remainder = j_index_dst(index) % (njGloDst-1);
+          if (0 == remainder) 
+          {
+            dataJIndex = -1;
+          }
+        }        
+        data_j_index_dst_full(index) = dataJIndex;
+
+        if ((0 <= data_i_index_dst_full(index)) && (0 <= data_j_index_dst_full(index)))
+        {
+          ++nbUnMaskedPointOnExtendedPart;
+        }
+     }
+  }
+
+  
+  // Finally, update data_i_index, data_j_index
+  int dataDstDim = domainDestination->data_dim;
+  data_i_index_dst.resize(nbUnMaskedPointOnExtendedPart);
+  data_j_index_dst.resize(nbUnMaskedPointOnExtendedPart); 
+  int count = 0; 
+  for (idx = 0; idx < newNbLocalDst; ++idx)
+  {
+    dataIIndex = data_i_index_dst_full(idx);
+    dataJIndex = data_j_index_dst_full(idx);
+    if ((0 <= dataIIndex) && (0 <= dataJIndex))
+    {
+      data_i_index_dst(count) = (1 == dataDstDim) ? idx : i_index_dst(idx) - i_index_dst(0);
+      data_j_index_dst(count) = (1 == dataDstDim) ? 0   : j_index_dst(idx) - j_index_dst(0);
+      ++count;
+    }
+  }
+
+  // Update data_ni, data_nj
+  
+  domainDestination->data_ni.setValue((1==dataDstDim) ? niDst * njDst : niDst);
+  domainDestination->data_nj.setValue((1==dataDstDim) ? niDst * njDst : njDst);  
+  domainDestination->data_ibegin.setValue(0);
+  domainDestination->data_jbegin.setValue(0);
+
+  // Update longitude and latitude 
+  if (niSrc == domainSource->lonvalue_1d.numElements() && njSrc == domainSource->latvalue_1d.numElements()) // Ok, we have rectilinear here
+  {
+     domainDestination->lonvalue_1d.resize(niDst);
+     domainDestination->lonvalue_1d = lon_dst(Range(0,niDst-1));
+     domainDestination->latvalue_1d.resize(njDst);
+     domainDestination->latvalue_1d = lat_dst(Range(0,lat_dst.numElements()-niDst,niDst));
+  }
+  else // It should be curvilinear
+  {
+     domainDestination->lonvalue_1d.resize(lon_dst.numElements());
+     domainDestination->lonvalue_1d = lon_dst;
+     domainDestination->latvalue_1d.resize(lat_dst.numElements());
+     domainDestination->latvalue_1d = (lat_dst);
+  }
+   domainDestination->mask_1d.resize(domainDestination->domainMask.numElements()) ;
+   domainDestination->mask_1d=domainDestination->domainMask ;
+   domainDestination->computeLocalMask() ;
+}
+CATCH
+
+/*!
+ *  Extend domain destination and update its attributes
+ *  Suppose that domain destination and domain source have the same values for all attributes (by inheritance)
+ *  \param [in/out] domainDestination domain destination
+ *  \param [in] domainSource domain source
+ *  \param [in] neighborsDomainSrc domain extended part
+ */
+void CDomainAlgorithmExpand::updateUnstructuredDomainAttributes(CDomain* domainDestination,
+                                                                CDomain* domainSource,
+                                                                CArray<int,2>& neighborsDomainSrc)
+TRY
+{
+
+  CContext* context = CContext::getCurrent();
+
+  // First of all, "copy" all attributes of domain source to domain destination
+  StdString domainDstRef = (!domainDestination->domain_ref.isEmpty()) ? domainDestination->domain_ref.getValue()
+                                                                      : "";
+  if (domainDstRef != domainSource->getId())
+  {
+    domainDestination->domain_ref.setValue(domainSource->getId());
+    domainDestination->solveRefInheritance(true);
+  }
+
+  if (domainDstRef.empty()) domainDestination->domain_ref.reset();
+  else domainDestination->domain_ref.setValue(domainDstRef);
+
+  // Now extend domain destination
+  int niGlob = domainSource->ni_glo;
+  CArray<bool,1>& mask_1d_src = domainSource->domainMask;
+  CArray<int,1>& i_index_src = domainSource->i_index;
+  CArray<double,1>& lon_src = domainSource->lonvalue_1d;
+  CArray<double,1>& lat_src = domainSource->latvalue_1d;
+  CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d;
+  CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d;
+  CArray<int,1>& data_i_index_src = domainSource->data_i_index;
+
+  int oldNbLocal = i_index_src.numElements(), index, globalIndex;
+  // Uncompress data_i_index
+  CArray<int,1> data_i_index_src_full(oldNbLocal);
+  int nbUnMaskedPointOnLocalDomain = 0;
+  data_i_index_src_full = -1; // Suppose all values are masked
+  for (int idx = 0; idx < data_i_index_src.numElements(); ++idx)
+  {
+    int dataIdx = data_i_index_src(idx);
+    if ((0 <= dataIdx) && (dataIdx < oldNbLocal))
+    {
+      data_i_index_src_full(nbUnMaskedPointOnLocalDomain) = dataIdx;
+      ++nbUnMaskedPointOnLocalDomain;
+    }
+  }
+
+  CArray<bool,1>& mask_1d_dst = domainDestination->domainMask;
+  CArray<int,1>& i_index_dst = domainDestination->i_index;
+  CArray<int,1>& j_index_dst = domainDestination->j_index;
+  CArray<double,1>& lon_dst = domainDestination->lonvalue_1d;
+  CArray<double,1>& lat_dst = domainDestination->latvalue_1d;
+  CArray<double,2>& bounds_lon_dst = domainDestination->bounds_lon_1d;
+  CArray<double,2>& bounds_lat_dst = domainDestination->bounds_lat_1d;
+  CArray<int,1>& data_i_index_dst = domainDestination->data_i_index;
+  CArray<int,1>& data_j_index_dst = domainDestination->data_j_index;
+
+  // Resize all array-like attributes of domain destination
+  int nbNeighbor    = neighborsDomainSrc.shape()[1];
+  int newNbLocalDst = nbNeighbor + oldNbLocal;
+  int nVertex       = bounds_lon_dst.shape()[0];
+
+  if (newNbLocalDst==0)
+  {
+    mask_1d_dst.resize(newNbLocalDst);
+    i_index_dst.resize(newNbLocalDst);
+    j_index_dst.resize(newNbLocalDst);
+    lon_dst.resize(newNbLocalDst);
+    lat_dst.resize(newNbLocalDst);
+    bounds_lon_dst.resize(nVertex, newNbLocalDst);
+    bounds_lat_dst.resize(nVertex, newNbLocalDst);
+  }
+  else
+  {
+    mask_1d_dst.resizeAndPreserve(newNbLocalDst);
+    i_index_dst.resizeAndPreserve(newNbLocalDst);
+    j_index_dst.resizeAndPreserve(newNbLocalDst);
+    lon_dst.resizeAndPreserve(newNbLocalDst);
+    lat_dst.resizeAndPreserve(newNbLocalDst);
+    bounds_lon_dst.resizeAndPreserve(nVertex, newNbLocalDst);
+    bounds_lat_dst.resizeAndPreserve(nVertex, newNbLocalDst);
+  }
+  CArray<int,1> data_i_index_dst_full(newNbLocalDst);
+  data_i_index_dst_full(Range(0,oldNbLocal-1)) = data_i_index_src_full;
+  data_i_index_dst_full(Range(oldNbLocal,newNbLocalDst-1)) = -1;
+
+  // 1. Fill in array relating to global index (i_index, j_index, transmap, etc, ...)
+  // Global index mapping between destination and source
+  this->transformationMapping_.resize(1);
+  this->transformationWeight_.resize(1);
+  TransformationIndexMap& transMap = this->transformationMapping_[0];
+  TransformationWeightMap& transWeight = this->transformationWeight_[0];
+
+  transMap.rehash(std::ceil(newNbLocalDst/transMap.max_load_factor()));
+  transWeight.rehash(std::ceil(newNbLocalDst/transWeight.max_load_factor()));
+  // First, index mapping for local domain
+  for (int idx = 0; idx < oldNbLocal; ++idx)
+  {
+    index = i_index_dst(idx);
+    transMap[index].push_back(index);
+    transWeight[index].push_back(1.0);
+  }
+  // Then, index mapping for extended part
+  for (int idx = 0; idx < nbNeighbor; ++idx)
+  {
+    index = idx + oldNbLocal;
+    globalIndex = neighborsDomainSrc(0,idx);
+    i_index_dst(index) = globalIndex;
+    j_index_dst(index) = 0;
+    transMap[globalIndex].push_back(globalIndex);
+    transWeight[globalIndex].push_back(1.0);
+  }
+
+  // 2. Exchange local info among domains (lon,lat,bounds,mask,etc,...)
+  CClientClientDHTDouble::Index2VectorInfoTypeMap localData;
+  localData.rehash(std::ceil(oldNbLocal/localData.max_load_factor()));
+  // Information exchanged among domains (attention to their order), number in parentheses presents size of data
+  // lon(1) + lat(1) + bounds_lon(nVertex) + bounds_lat(nVertex) + mask(1) + data_i_index(1)
+  int dataPackageSize =  1 + 1 + // lon + lat
+                         nVertex + nVertex + //bounds_lon + bounds_lat
+                         1 + // mask_1d_dst;
+                         1; // data_i_index
+  // Initialize database
+  for (int idx = 0; idx < oldNbLocal; ++idx)
+  {
+    index = i_index_src(idx);
+    localData[index].resize(dataPackageSize);
+    std::vector<double>& data = localData[index];
+
+    //Pack data
+    int dataIdx = 0;
+    data[dataIdx] = lon_src(idx);++dataIdx;
+    data[dataIdx] = lat_src(idx);++dataIdx;
+    for (int i = 0; i < nVertex; ++i)
+    {
+      data[dataIdx] = bounds_lon_src(i,idx); ++dataIdx;
+    }
+    for (int i = 0; i < nVertex; ++i)
+    {
+      data[dataIdx] = bounds_lat_src(i,idx); ++dataIdx;
+    }
+    data[dataIdx] = mask_1d_src(idx) ? 1.0 : -1.0; ++dataIdx;
+    data[dataIdx] = data_i_index_src_full(idx);
+  }
+
+  CClientClientDHTDouble dhtData(localData, context->intraComm_);
+  CArray<size_t,1> neighborInd(nbNeighbor);
+  for (int idx = 0; idx < nbNeighbor; ++idx)
+    neighborInd(idx) = neighborsDomainSrc(0,idx);
+
+  // Compute local data on other domains
+  dhtData.computeIndexInfoMapping(neighborInd);
+  CClientClientDHTDouble::Index2VectorInfoTypeMap& neighborData = dhtData.getInfoIndexMap();
+  CClientClientDHTDouble::Index2VectorInfoTypeMap::iterator ite = neighborData.end(), it;
+  // Ok get neighbor data
+  size_t nIdx;
+  int nbUnMaskedPointOnExtendedPart = 0;
+  for (int idx = 0; idx < nbNeighbor; ++idx)
+  {
+    nIdx  = neighborInd(idx);
+    it = neighborData.find(nIdx);
+    if (ite != it)
+    {
+      index = idx + oldNbLocal;
+      std::vector<double>& data = it->second;
+      // Unpack data
+      int dataIdx = 0;
+      lon_dst(index) = data[dataIdx]; ++dataIdx;
+      lat_dst(index) = data[dataIdx]; ++dataIdx;
+      for (int i = 0; i < nVertex; ++i)
+      {
+        bounds_lon_dst(i,index) = data[dataIdx]; ++dataIdx;
+      }
+      for (int i = 0; i < nVertex; ++i)
+      {
+        bounds_lat_dst(i,index) = data[dataIdx]; ++dataIdx;
+      }
+      mask_1d_dst(index) = (1.0 == data[dataIdx]) ? true : false; ++dataIdx;
+      data_i_index_dst_full(index) = (int)(data[dataIdx]);
+      if (0 <= data_i_index_dst_full(index))
+      {
+        data_i_index_dst_full(index) = index;
+        ++nbUnMaskedPointOnExtendedPart;
+      }
+    }
+  }
+
+
+  // Finally, update data_i_index
+  int nbUnMaskedPointOnNewDstDomain = (nbUnMaskedPointOnExtendedPart + nbUnMaskedPointOnLocalDomain);
+  int count = 0, dataIdx;
+  for (int idx = 0; idx < newNbLocalDst; ++idx)
+  {
+    dataIdx = data_i_index_dst_full(idx);
+    if ((0 <= dataIdx))
+    {
+      ++count;
+    }
+  }
+
+  data_i_index_dst.resize(count);
+  data_j_index_dst.resize(count);
+  data_j_index_dst = 0;
+
+  count = 0;
+  for (int idx = 0; idx < newNbLocalDst; ++idx)
+  {
+    dataIdx = data_i_index_dst_full(idx);
+    if ((0 <= dataIdx))
+    {
+      data_i_index_dst(count) = dataIdx;
+      ++count;
+    }
+  }
+
+  // Update ni
+  domainDestination->ni.setValue(newNbLocalDst);
+  domainDestination->mask_1d.resize(domainDestination->domainMask.numElements()) ;
+  domainDestination->mask_1d=domainDestination->domainMask ;
+  domainDestination->computeLocalMask() ;
+}
+CATCH
+
+/*!
+  Compute the index mapping between domain on grid source and one on grid destination
+*/
+void CDomainAlgorithmExpand::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+{
+
+}
+
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_expand.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_expand.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_expand.hpp	(revision 1985)
@@ -0,0 +1,63 @@
+/*!
+   \file domain_algorithm_expand.hpp
+   \author Ha NGUYEN
+   \since 08 Aug 2016
+   \date 08 Aug 2016
+
+   \brief Algorithm for expanding an domain.
+ */
+#ifndef __XIOS_DOMAIN_ALGORITHM_EXPAND_HPP__
+#define __XIOS_DOMAIN_ALGORITHM_EXPAND_HPP__
+
+#include "domain_algorithm_transformation.hpp"
+#include "transformation.hpp"
+
+namespace xios {
+class CDomain;
+class CExpandDomain;
+
+/*!
+  \class CDomainAlgorithmExpand
+*/
+class CDomainAlgorithmExpand : public CDomainAlgorithmTransformation
+{
+public:
+  CDomainAlgorithmExpand(CDomain* domainDestination, CDomain* domainSource, CExpandDomain* expandDomain);
+
+  virtual ~CDomainAlgorithmExpand() {}
+
+  static bool registerTrans();
+
+protected:
+  bool isXPeriodic_; // Flag to determine the periodicity of expansion (only for rectilinear)
+  bool isYPeriodic_; // Flag to determine the periodicity of expansion (only for rectilinear)
+
+protected:
+  void expandDomainEdgeConnectivity(CDomain* domainDestination, CDomain* domainSource);
+  void expandDomainNodeConnectivity(CDomain* domainDestination, CDomain* domainSource);
+  void updateRectilinearDomainAttributes(CDomain* domainDestination,
+                                         CDomain* domainSource,
+                                         CArray<int,2>& neighborsDomainSrc);
+
+  void updateUnstructuredDomainAttributes(CDomain* domainDestination,
+                                          CDomain* domainSource,
+                                          CArray<int,2>& neighborsDomainSrc);
+
+protected:
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs);
+
+private:
+  static CGenericAlgorithmTransformation* create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                CTransformation<CDomain>* transformation,
+                                                int elementPositionInGrid,
+                                                std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridDst2DomainPosition);
+  static bool dummyRegistered_;
+};
+
+}
+#endif // __XIOS_DOMAIN_ALGORITHM_EXPAND_HPP__
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_extract.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_extract.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_extract.cpp	(revision 1985)
@@ -0,0 +1,252 @@
+#include "domain_algorithm_extract.hpp"
+#include "extract_domain.hpp"
+#include "domain.hpp"
+#include "grid.hpp"
+#include "grid_transformation_factory_impl.hpp"
+#include "attribute_template.hpp"
+
+namespace xios {
+CGenericAlgorithmTransformation* CDomainAlgorithmExtract::create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                             CTransformation<CDomain>* transformation,
+                                                             int elementPositionInGrid,
+                                                             std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                             std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                             std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                             std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                             std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                             std::map<int, int>& elementPositionInGridDst2DomainPosition)
+TRY
+{
+  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
+  std::vector<CDomain*> domainListSrcP  = gridSrc->getDomains();
+
+  CExtractDomain* extractDomain = dynamic_cast<CExtractDomain*> (transformation);
+  int domainDstIndex = elementPositionInGridDst2DomainPosition[elementPositionInGrid];
+  int domainSrcIndex = elementPositionInGridSrc2DomainPosition[elementPositionInGrid];
+
+  return (new CDomainAlgorithmExtract(domainListDestP[domainDstIndex], domainListSrcP[domainSrcIndex], extractDomain));
+}
+CATCH
+
+bool CDomainAlgorithmExtract::dummyRegistered_ = CDomainAlgorithmExtract::registerTrans();
+bool CDomainAlgorithmExtract::registerTrans()
+TRY
+{
+  return CGridTransformationFactory<CDomain>::registerTransformation(TRANS_EXTRACT_DOMAIN, create);
+}
+CATCH
+
+CDomainAlgorithmExtract::CDomainAlgorithmExtract(CDomain* domainDestination, CDomain* domainSource, CExtractDomain* extractDomain)
+: CDomainAlgorithmTransformation(domainDestination, domainSource)
+TRY
+{
+  extractDomain->checkValid(domainSource);
+  extractIBegin_ = extractDomain->ibegin.getValue();
+  extractJBegin_ = extractDomain->jbegin.getValue();
+
+  extractNi_  = extractDomain->ni.getValue();
+  extractNj_  = extractDomain->nj.getValue();
+
+  extractIEnd_ = extractIBegin_ + extractNi_ - 1;
+  extractJEnd_ = extractJBegin_ + extractNj_ - 1;
+
+  if (extractNi_ > domainSource->ni_glo.getValue())
+  {
+    ERROR("CDomainAlgorithmExtract::CDomainAlgorithmExtract(CDomain* domainDestination, CDomain* domainSource, CExtractDomain* extractDomain)",
+           << "Extract size is greater than size of domain source"
+           << "Size ni_glo of domain source " <<domainSource->getId() << " is " << domainSource->ni_glo.getValue()  << std::endl
+           << "Extract size is " << extractNi_ );
+  }
+
+  if (extractNj_ > domainSource->nj_glo.getValue())
+  {
+    ERROR("CDomainAlgorithmExtract::CDomainAlgorithmExtract(CDomain* domainDestination, CDomain* domainSource, CExtractDomain* extractDomain)",
+           << "Extract size is greater than size of domain source"
+           << "Size nj_glo of domain source " <<domainSource->getId() << " is " << domainSource->nj_glo.getValue()  << std::endl
+           << "Extract size is " << extractNj_ );
+  }
+
+  // Calculate the size of local domain
+  int ind, indLocSrc, indLocDest, iIdxSrc, jIdxSrc, destIBegin = -1, destJBegin = -1, niDest = 0, njDest = 0, ibeginDest, jbeginDest ;
+  int indGloDest, indGloSrc, niGloSrc = domainSrc_->ni_glo, iSrc, jSrc;
+  for (int j = 0; j < domainSrc_->nj.getValue(); j++)
+  {
+    for (int i = 0; i < domainSrc_->ni.getValue(); i++)
+    {
+      ind = j*domainSrc_->ni + i;
+      iIdxSrc = domainSrc_->i_index(ind);
+      if ((iIdxSrc >= extractIBegin_) && (iIdxSrc <= extractIEnd_))
+      {
+        jIdxSrc = domainSrc_->j_index(ind);
+        if ((jIdxSrc >= extractJBegin_) && (jIdxSrc <= extractJEnd_))
+        {
+          if ((niDest == 0) && (njDest == 0))
+          {
+            destIBegin = i;
+            destJBegin = j;
+          }
+          if (i == destIBegin) ++njDest;
+        }
+        if (j == destJBegin) ++niDest;
+
+      }
+    }
+  }
+  ibeginDest = destIBegin + domainSrc_->ibegin - extractIBegin_;
+  jbeginDest = destJBegin + domainSrc_->jbegin - extractJBegin_;
+  if (niDest==0) ibeginDest=0 ;
+  if (njDest==0) jbeginDest=0 ;
+  domainDest_->ni_glo.setValue(extractNi_);
+  domainDest_->nj_glo.setValue(extractNj_);
+  domainDest_->ni.setValue(niDest);
+  domainDest_->nj.setValue(njDest);
+  domainDest_->ibegin.setValue(ibeginDest);
+  domainDest_->jbegin.setValue(jbeginDest);
+  domainDest_->i_index.resize(niDest*njDest);
+  domainDest_->j_index.resize(niDest*njDest);
+
+  domainDest_->data_ni.setValue(niDest);
+  domainDest_->data_nj.setValue(njDest);
+  domainDest_->data_ibegin.setValue(0);  // local position
+  domainDest_->data_jbegin.setValue(0);  // local position
+  domainDest_->data_i_index.resize(niDest*njDest); // local position
+  domainDest_->data_j_index.resize(niDest*njDest); // local position
+
+  domainDest_->domainMask.resize(niDest*njDest);
+
+  if (!domainSrc_->mask_2d.isEmpty()) domainDest_->mask_2d.resize(niDest,njDest);
+  if (!domainSrc_->mask_1d.isEmpty()) domainDest_->mask_1d.resize(niDest*njDest);
+
+  if (!domainSrc_->lonvalue_1d.isEmpty())
+  {
+    if (domainDest_->type == CDomain::type_attr::rectilinear)
+    {
+      domainDest_->lonvalue_1d.resize(niDest);
+      domainDest_->latvalue_1d.resize(njDest);
+    }
+    else if (domainDest_->type == CDomain::type_attr::unstructured)
+    {
+      domainDest_->lonvalue_1d.resize(niDest);
+      domainDest_->latvalue_1d.resize(niDest);
+    }
+  }
+  else if (!domainSrc_->lonvalue_2d.isEmpty())
+  {
+    domainDest_->lonvalue_2d.resize(niDest,njDest);
+    domainDest_->latvalue_2d.resize(niDest,njDest);
+  }
+
+  if (domainSrc_->hasBounds)
+  {
+    if (!domainSrc_->bounds_lon_2d.isEmpty())
+    {
+      domainDest_->bounds_lon_2d.resize(domainDest_->nvertex, niDest, njDest);
+      domainDest_->bounds_lon_2d.resize(domainDest_->nvertex, niDest, njDest);
+    }
+    else if (!domainSrc_->bounds_lon_1d.isEmpty())
+    {
+      domainDest_->bounds_lon_1d.resize(domainDest_->nvertex, niDest);
+      domainDest_->bounds_lon_1d.resize(domainDest_->nvertex, niDest);
+    }
+  }
+  if (domainSrc_->hasArea) domainDest_->area.resize(niDest,njDest);
+
+  this->transformationMapping_.resize(1);
+  this->transformationWeight_.resize(1);
+  TransformationIndexMap& transMap = this->transformationMapping_[0];
+  TransformationWeightMap& transWeight = this->transformationWeight_[0];
+
+  for (int iDest = 0; iDest < niDest; iDest++)
+  {
+    iSrc = iDest + destIBegin;
+    for (int jDest = 0; jDest < njDest; jDest++)
+    {
+      jSrc = jDest + destJBegin;
+      ind = jSrc * domainSrc_->ni + iSrc;
+      iIdxSrc = domainSrc_->i_index(ind);
+      jIdxSrc = domainSrc_->j_index(ind);
+      indLocDest = jDest*niDest + iDest;
+      indGloDest = (jDest + jbeginDest)*extractNi_ + (iDest + ibeginDest);
+      indLocSrc = (jDest+destJBegin)*domainSrc_->ni + (iDest+destIBegin);
+      indGloSrc = (jIdxSrc )* niGloSrc + iIdxSrc;
+      domainDest_->i_index(indLocDest) = iDest + ibeginDest;                                             // i_index contains global positions
+      domainDest_->j_index(indLocDest) = jDest + jbeginDest;                                             // i_index contains global positions
+      domainDest_->data_i_index(indLocDest) = (domainSrc_->data_dim == 1) ? indLocDest : iDest;          // data_i_index contains local positions
+      domainDest_->data_j_index(indLocDest) = (domainSrc_->data_dim == 1) ? 0 :jDest;                    // data_i_index contains local positions
+      domainDest_->domainMask(indLocDest) = domainSrc_->domainMask(indLocSrc);
+
+      if (!domainSrc_->mask_2d.isEmpty())
+        domainDest_->mask_2d(iDest,jDest) = domainSrc_->mask_2d(iSrc,jSrc);
+
+      if (!domainSrc_->mask_1d.isEmpty())
+        domainDest_->mask_1d(indLocDest) = domainSrc_->mask_1d(indLocSrc);
+
+      if (domainSrc_->hasArea)
+        domainDest_->area(iDest,jDest) = domainSrc_->area(iSrc,jSrc);
+
+
+      if (domainSrc_->hasBounds)
+      {
+        if (!domainSrc_->bounds_lon_2d.isEmpty())
+        {
+          for (int n = 0; n < domainSrc_->nvertex; ++n)
+          {
+            domainDest_->bounds_lon_2d(n,iDest,jDest) = domainSrc_->bounds_lon_2d(n,iSrc,jSrc);
+            domainDest_->bounds_lat_2d(n,iDest,jDest) = domainSrc_->bounds_lat_2d(n,iSrc,jSrc);
+          }
+        }
+        else if (!domainSrc_->bounds_lon_1d.isEmpty())
+        {
+          for (int n = 0; n < domainSrc_->nvertex; ++n)
+          {
+            domainDest_->bounds_lon_1d(n,iDest) = domainSrc_->bounds_lon_1d(n,iSrc);
+            domainDest_->bounds_lat_1d(n,iDest) = domainSrc_->bounds_lat_1d(n,iSrc);
+          }
+        }
+      }
+
+      if (domainSrc_->hasLonLat)
+      {
+        if (domainDest_->type == CDomain::type_attr::rectilinear)
+        {
+          domainDest_->latvalue_1d(jDest) = domainSrc_->latvalue_1d(jSrc);
+        }
+        else if (domainDest_->type == CDomain::type_attr::curvilinear)
+        {
+          domainDest_->lonvalue_2d(iDest,jDest) = domainSrc_->lonvalue_2d(iSrc,jSrc);
+          domainDest_->latvalue_2d(iDest,jDest) = domainSrc_->latvalue_2d(iSrc,jSrc);
+        }
+      }
+
+      transMap[indGloDest].push_back(indGloSrc);
+      transWeight[indGloDest].push_back(1.0);
+
+    }
+    if (domainSrc_->hasLonLat)
+    {
+      if (domainDest_->type == CDomain::type_attr::unstructured)
+      {
+        domainDest_->lonvalue_1d(iDest) = domainSrc_->lonvalue_1d(iSrc);
+        domainDest_->latvalue_1d(iDest) = domainSrc_->latvalue_1d(iSrc);
+      }
+      else if (domainDest_->type == CDomain::type_attr::rectilinear)
+      {
+        domainDest_->lonvalue_1d(iDest) = domainSrc_->lonvalue_1d(iSrc);
+      }
+    }
+  }
+  
+  domainDestination->checkAttributes() ;
+  computeAlgorithm(domainSource->getLocalView(CElementView::WORKFLOW), domainDestination->getLocalView(CElementView::WORKFLOW)) ;
+}
+CATCH
+
+/*!
+  Compute the index mapping between domain on grid source and one on grid destination
+*/
+void CDomainAlgorithmExtract::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+{
+}
+
+
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_extract.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_extract.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_extract.hpp	(revision 1985)
@@ -0,0 +1,60 @@
+
+#ifndef __XIOS_DOMAIN_ALGORITHM_EXTRACT_HPP__
+#define __XIOS_DOMAIN_ALGORITHM_EXTRACT_HPP__
+
+#include "domain_algorithm_transformation.hpp"
+#include "transformation.hpp"
+
+namespace xios {
+
+class CDomain;
+class CExtractDomain;
+
+/*!
+  \class CDomainAlgorithmExtract
+  Implementing extract (alternative zoom) on domain
+*/
+class CDomainAlgorithmExtract : public CDomainAlgorithmTransformation
+{
+public:
+  CDomainAlgorithmExtract(CDomain* domainDestination, CDomain* domainSource, CExtractDomain* extractDomain);
+
+  virtual ~CDomainAlgorithmExtract() {}
+
+  static bool registerTrans();
+protected:
+  void updateDomainAttributes();
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs);
+
+private:
+  void updateExtract();
+
+private:
+  //! Global extract begin on domain
+  int extractIBegin_;
+  int extractJBegin_;
+
+  //! Global extract end on domain
+  int extractIEnd_;
+  int extractJEnd_;
+
+  //! Global extract size on domain
+  int extractNi_;
+  int extractNj_;
+
+private:
+
+  static CGenericAlgorithmTransformation* create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                CTransformation<CDomain>* transformation,
+                                                int elementPositionInGrid,
+                                                std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridDst2DomainPosition);
+  static bool dummyRegistered_;
+};
+
+}
+#endif // __XIOS_DOMAIN_ALGORITHM_EXTRACT_HPP__
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_generate_rectilinear.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_generate_rectilinear.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_generate_rectilinear.cpp	(revision 1985)
@@ -0,0 +1,141 @@
+/*!
+   \file domain_algorithm_generate_rectilinear.cpp
+   \author Ha NGUYEN
+   \since 31 Aug 2015
+   \date 31 Aug 2015
+
+   \brief Algorithm for automatic generation of rectilinear domain.
+ */
+#include "domain_algorithm_generate_rectilinear.hpp"
+#include "grid.hpp"
+#include "domain.hpp"
+#include "context.hpp"
+#include "context_client.hpp"
+#include "generate_rectilinear_domain.hpp"
+
+namespace xios {
+
+CDomainAlgorithmGenerateRectilinear::CDomainAlgorithmGenerateRectilinear(CDomain* domainDestination, CDomain* domainSource,
+                                                                         CGrid* gridDest, CGrid* gridSource,
+                                                                         CGenerateRectilinearDomain* genRectDomain)
+: CDomainAlgorithmTransformation(domainDestination, domainSource), nbDomainDistributedPart_(0)
+TRY
+{
+  type_ = ELEMENT_GENERATION;
+  genRectDomain->checkValid(domainDestination);
+  if (0 != gridSource) computeDistributionGridSource(gridSource);
+  else
+  {
+    computeDistributionGridDestination(gridDest);
+  }
+  fillInAttributesDomainDestination();
+}
+CATCH
+
+/*!
+  Compute the index mapping between domain on grid source and one on grid destination
+*/
+void CDomainAlgorithmGenerateRectilinear::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+{
+ /* Nothing to do */
+}
+
+/*!
+  Calculate the number of distributed parts on domain source
+*/
+void CDomainAlgorithmGenerateRectilinear::computeDistributionGridSource(CGrid* gridSrc)
+TRY
+{
+  CContext* context = CContext::getCurrent();
+  int clientSize = context->intraCommSize_ ;
+  int clientRank = context->intraCommRank_ ;
+  
+  std::vector<CDomain*> domListSrcP = gridSrc->getDomains();
+  std::vector<CAxis*> axisListSrcP = gridSrc->getAxis();
+
+  for (int i = 0; i < domListSrcP.size(); ++i) // support we have only domain, more than one, for now, dont know how to process
+  {
+    // First, find (roundly) distribution of associated axis (if any)
+    if (axisListSrcP.empty()) nbDomainDistributedPart_ = clientSize;
+    else
+    {
+      gridSrc->solveAxisRef(false);
+      int nbAxis = axisListSrcP.size();
+      std::vector<int> nbLocalAxis(nbAxis, 0);
+      for (int j = 0; j < nbAxis; ++j)
+      {
+        std::vector<int> globalAxisIndex(axisListSrcP[j]->n);
+        for (int idx = 0; idx < axisListSrcP[j]->n; ++idx)
+          globalAxisIndex[idx] = axisListSrcP[j]->begin + idx;
+        HashXIOS<int> hashFunc;
+        StdSize hashValue = hashFunc.hashVec(globalAxisIndex);
+        std::vector<StdSize> recvBuff(clientSize);
+        MPI_Gather(&hashValue, 1, MPI_UNSIGNED_LONG,
+                   &recvBuff[0], 1, MPI_UNSIGNED_LONG,
+                   0,
+                   context->intraComm_);
+        if (0 == clientRank)
+        {
+          std::set<StdSize> setTmp;
+          for (int k = 0; k < recvBuff.size(); ++k)
+          {
+            if (setTmp.end() == setTmp.find(recvBuff[k]))
+            {
+              ++nbLocalAxis[j];
+              setTmp.insert(recvBuff[k]);
+            }
+          }
+        }
+
+        MPI_Bcast(&nbLocalAxis[0], nbAxis, MPI_INT,
+                  0, context->intraComm_);
+      }
+
+      int nbAxisDistributedPart = 1;
+      for (int j = 0; j < nbAxis; ++j) nbAxisDistributedPart *= nbLocalAxis[j];
+      nbDomainDistributedPart_ = clientSize/nbAxisDistributedPart;
+    }
+  }
+}
+CATCH
+
+/*!
+  Compute the distribution of the domain destination by using available information provided by user such as n_distributed_partition of axis
+*/
+void CDomainAlgorithmGenerateRectilinear::computeDistributionGridDestination(CGrid* gridDest)
+TRY
+{
+  // For now, just suppose that the grid contains only one domain
+  std::vector<CAxis*> axisListDestP = gridDest->getAxis();
+  int nbPartition = 1, idx = 0;
+  for (int i = 0; i < gridDest->axis_domain_order.numElements(); ++i)
+  {
+    if (false == (gridDest->axis_domain_order)(i))
+    {
+      nbPartition *= (axisListDestP[idx]->n_distributed_partition.isEmpty()) ? 1: (axisListDestP[idx]->n_distributed_partition.getValue());
+      ++idx;
+    }
+  }
+
+  CContext* context = CContext::getCurrent();
+  int modPart = (context->intraCommSize_) % nbPartition;
+  if (0 != modPart)
+    ERROR("CDomainAlgorithmGenerateRectilinear::computeDistributionGridDestination(CGrid* gridDest)",
+       << "The grid " <<gridDest->getId() << " is not well-distributed. There is an incompatibility between distribution of axis and domain.");
+  nbDomainDistributedPart_ = context->intraCommSize_/nbPartition;
+
+}
+CATCH
+
+/*!
+  Fill in all necessary attributes of domain destination and their values
+*/
+void CDomainAlgorithmGenerateRectilinear::fillInAttributesDomainDestination()
+TRY
+{
+  if (!domainDest_->distributionAttributesHaveValue())
+    domainDest_->redistribute(nbDomainDistributedPart_);
+  domainDest_->fillInLonLat();
+}
+CATCH
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_generate_rectilinear.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_generate_rectilinear.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_generate_rectilinear.hpp	(revision 1985)
@@ -0,0 +1,50 @@
+/*!
+   \file domain_algorithm_generate_rectilinear.hpp
+   \author Ha NGUYEN
+   \since 31 Aug 2015
+   \date 01 Sep 2015
+
+   \brief Algorithm for automatic generation of rectilinear domain.
+ */
+#ifndef __XIOS_DOMAIN_ALGORITHM_GENERATE_RECTILINEAR_HPP__
+#define __XIOS_DOMAIN_ALGORITHM_GENERATE_RECTILINEAR_HPP__
+
+#include "domain_algorithm_transformation.hpp"
+
+namespace xios {
+
+class CGrid;
+class CDomain;
+class CGenerateRectilinearDomain;
+
+/*!
+  \class CDomainAlgorithmGenerateRectilinear
+  Generate a rectilinear or CURVILINEAR domain and fill in necessary its attributes automatically
+  A new rectilinear (or CURVILINEAR) domain will also be distributed automatically among the processes.
+  The number of processes is deduced from the distribution of the grid source.
+*/
+class CDomainAlgorithmGenerateRectilinear : public CDomainAlgorithmTransformation
+{
+public:
+  CDomainAlgorithmGenerateRectilinear(CDomain* domainDestination, CDomain* domainSource,
+                                      CGrid* gridDest, CGrid* gridSource,
+                                      CGenerateRectilinearDomain* zoomDomain);
+
+  virtual ~CDomainAlgorithmGenerateRectilinear() {}
+
+protected:
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs);
+
+private:
+  void computeDistributionGridSource(CGrid* gridSrc);
+  void computeDistributionGridDestination(CGrid* gridDest);
+  void fillInAttributesDomainDestination();
+
+private:
+  int nbDomainDistributedPart_; //! Number of local domain.
+
+  static bool dummyRegistered_;
+};
+
+}
+#endif // __XIOS_DOMAIN_ALGORITHM_GENERATE_RECTILINEAR_HPP__
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_interpolate.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_interpolate.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_interpolate.cpp	(revision 1985)
@@ -0,0 +1,1026 @@
+/*!
+   \file domain_algorithm_interpolate_from_file.cpp
+   \author Ha NGUYEN
+   \since 09 Jul 2015
+   \date 15 Sep 2015
+
+   \brief Algorithm for interpolation on a domain.
+ */
+#include "domain_algorithm_interpolate.hpp"
+#include <unordered_map>
+#include "context.hpp"
+#include "context_client.hpp"
+#include "distribution_client.hpp"
+#include "client_server_mapping_distributed.hpp"
+#include "netcdf.hpp"
+#include "mapper.hpp"
+#include "mpi_tag.hpp"
+#include "domain.hpp"
+#include "grid_transformation_factory_impl.hpp"
+#include "interpolate_domain.hpp"
+#include "grid.hpp"
+
+namespace xios {
+CGenericAlgorithmTransformation* CDomainAlgorithmInterpolate::create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                                     CTransformation<CDomain>* transformation,
+                                                                     int elementPositionInGrid,
+                                                                     std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                                     std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                                     std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2DomainPosition)
+TRY
+{
+  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
+  std::vector<CDomain*> domainListSrcP  = gridSrc->getDomains();
+
+  CInterpolateDomain* interpolateDomain = dynamic_cast<CInterpolateDomain*> (transformation);
+  int domainDstIndex = elementPositionInGridDst2DomainPosition[elementPositionInGrid];
+  int domainSrcIndex = elementPositionInGridSrc2DomainPosition[elementPositionInGrid];
+
+  return (new CDomainAlgorithmInterpolate(domainListDestP[domainDstIndex], domainListSrcP[domainSrcIndex], interpolateDomain));
+}
+CATCH
+
+bool CDomainAlgorithmInterpolate::dummyRegistered_ = CDomainAlgorithmInterpolate::registerTrans();
+bool CDomainAlgorithmInterpolate::registerTrans()
+TRY
+{
+  return CGridTransformationFactory<CDomain>::registerTransformation(TRANS_INTERPOLATE_DOMAIN, create);
+}
+CATCH
+
+CDomainAlgorithmInterpolate::CDomainAlgorithmInterpolate(CDomain* domainDestination, CDomain* domainSource, CInterpolateDomain* interpDomain)
+: CDomainAlgorithmTransformation(domainDestination, domainSource), interpDomain_(interpDomain), writeToFile_(false), readFromFile_(false)
+TRY
+{
+  CContext* context = CContext::getCurrent();
+  interpDomain_->checkValid(domainSource);
+
+  detectMissingValue = interpDomain_->detect_missing_value ;
+  renormalize = interpDomain_->renormalize ;
+  quantity = interpDomain_->quantity ;
+
+  if (interpDomain_->read_write_convention == CInterpolateDomain::read_write_convention_attr::fortran) fortranConvention=true ;
+  else fortranConvention=false ;
+  
+  fileToReadWrite_ = "xios_interpolation_weights_";
+
+  if (interpDomain_->weight_filename.isEmpty())
+  {
+    fileToReadWrite_ += context->getId() + "_" + 
+                    domainSource->getDomainOutputName() + "_" + 
+                    domainDestination->getDomainOutputName() + ".nc";    
+  }
+  else 
+    fileToReadWrite_ = interpDomain_->weight_filename;
+
+  ifstream f(fileToReadWrite_.c_str());  
+  switch (interpDomain_->mode)
+  {
+    case CInterpolateDomain::mode_attr::read:
+      readFromFile_ = true;      
+      break;
+    case CInterpolateDomain::mode_attr::compute:
+      readFromFile_ = false;
+      break;
+    case CInterpolateDomain::mode_attr::read_or_compute:      
+      if (!f.good())
+        readFromFile_ = false;
+      else
+        readFromFile_ = true;
+      break;
+    default:
+      break;
+  } 
+
+  writeToFile_ = interpDomain_->write_weight;  
+    
+}
+CATCH
+
+/*!
+  Compute remap with integrated remap calculation module
+*/
+void CDomainAlgorithmInterpolate::computeRemap()
+TRY
+{
+  using namespace sphereRemap;
+
+  CContext* context = CContext::getCurrent();
+  int clientRank = context->intraCommRank_;
+  int i, j, k, idx;
+  std::vector<double> srcPole(3,0), dstPole(3,0);
+  int orderInterp = interpDomain_->order.getValue();
+
+
+  const double poleValue = 90.0;
+  const int constNVertex = 4; // Value by default number of vertex for rectangular domain
+  int nVertexSrc, nVertexDest;
+  nVertexSrc = nVertexDest = constNVertex;
+
+  // First of all, try to retrieve the boundary values of domain source and domain destination
+  int localDomainSrcSize = domainSrc_->i_index.numElements();
+  int niSrc = domainSrc_->ni.getValue(), njSrc = domainSrc_->nj.getValue();
+  bool hasBoundSrc = domainSrc_->hasBounds;
+  if (hasBoundSrc) nVertexSrc = domainSrc_->nvertex.getValue();
+  CArray<double,2> boundsLonSrc(nVertexSrc,localDomainSrcSize);
+  CArray<double,2> boundsLatSrc(nVertexSrc,localDomainSrcSize);
+
+  if (domainSrc_->hasPole) srcPole[2] = 1;
+  if (hasBoundSrc)  // Suppose that domain source is curvilinear or unstructured
+  {
+    if (!domainSrc_->bounds_lon_2d.isEmpty())
+    {
+      for (j = 0; j < njSrc; ++j)
+        for (i = 0; i < niSrc; ++i)
+        {
+          k=j*niSrc+i;
+          for(int n=0;n<nVertexSrc;++n)
+          {
+            boundsLonSrc(n,k) = domainSrc_->bounds_lon_2d(n,i,j);
+            boundsLatSrc(n,k) = domainSrc_->bounds_lat_2d(n,i,j);
+          }
+        }
+    }
+    else
+    {
+      boundsLonSrc = domainSrc_->bounds_lon_1d;
+      boundsLatSrc = domainSrc_->bounds_lat_1d;
+    }
+  }
+  else // if domain source is rectilinear, not do anything now
+  {
+    CArray<double,1> lon_g ;
+    CArray<double,1> lat_g ;
+
+    if (!domainSrc_->lonvalue_1d.isEmpty() && !domainSrc_->latvalue_1d.isEmpty())
+    {
+      domainSrc_->AllgatherRectilinearLonLat(domainSrc_->lonvalue_1d,domainSrc_->latvalue_1d, lon_g,lat_g) ;
+    }
+    else if (! domainSrc_->latvalue_rectilinear_read_from_file.isEmpty() && ! domainSrc_->lonvalue_rectilinear_read_from_file.isEmpty() )
+    {
+      lat_g=domainSrc_->latvalue_rectilinear_read_from_file ;
+      lon_g=domainSrc_->lonvalue_rectilinear_read_from_file ;
+    }
+    else if (!domainSrc_->lon_start.isEmpty() && !domainSrc_->lon_end.isEmpty() &&
+             !domainSrc_->lat_start.isEmpty() && !domainSrc_->lat_end.isEmpty())
+    {
+      double step=(domainSrc_->lon_end-domainSrc_->lon_start)/domainSrc_->ni_glo ;
+      for (int i=0; i<domainSrc_->ni_glo; ++i) lon_g(i)=domainSrc_->lon_start+i*step ;
+      step=(domainSrc_->lat_end-domainSrc_->lat_start)/domainSrc_->nj_glo ;
+      for (int i=0; i<domainSrc_->ni_glo; ++i) lat_g(i)=domainSrc_->lat_start+i*step ;
+    }
+    else if (!domainSrc_->lonvalue.isEmpty() && !domainSrc_->latvalue.isEmpty())  // data can come possibly from coupler so only lonvalue and latvalue are defined ?
+    {
+      CArray<double,1> lon(niSrc), lat(njSrc) ;
+      for(int i=0; i< niSrc;++i) lon(i) = domainSrc_->lonvalue(i) ;
+      for(int j=0, n=0; j<njSrc; ++j, n+=niSrc) lat(j) = domainSrc_->latvalue(n) ;
+      domainSrc_->AllgatherRectilinearLonLat(lon, lat, lon_g, lat_g) ;
+    }
+    else ERROR("void CDomainAlgorithmInterpolate::computeRemap()",<<"Cannot compute bounds for rectilinear domain") ;
+
+    nVertexSrc = constNVertex;
+    domainSrc_->fillInRectilinearBoundLonLat(lon_g,lat_g, boundsLonSrc, boundsLatSrc);
+  }
+
+  std::map<int,std::vector<std::pair<int,double> > > interpMapValueNorthPole;
+  std::map<int,std::vector<std::pair<int,double> > > interpMapValueSouthPole;
+
+  int localDomainDestSize = domainDest_->i_index.numElements();
+  int niDest = domainDest_->ni.getValue(), njDest = domainDest_->nj.getValue();
+  bool hasBoundDest = domainDest_->hasBounds;
+  if (hasBoundDest) nVertexDest = domainDest_->nvertex.getValue();
+  CArray<double,2> boundsLonDest(nVertexDest,localDomainDestSize);
+  CArray<double,2> boundsLatDest(nVertexDest,localDomainDestSize);
+
+  if (domainDest_->hasPole) dstPole[2] = 1;
+  if (hasBoundDest)
+  {
+    if (!domainDest_->bounds_lon_2d.isEmpty())
+    {
+      for (j = 0; j < njDest; ++j)
+        for (i = 0; i < niDest; ++i)
+        {
+          k=j*niDest+i;
+          for(int n=0;n<nVertexDest;++n)
+          {
+            boundsLonDest(n,k) = domainDest_->bounds_lon_2d(n,i,j);
+            boundsLatDest(n,k) = domainDest_->bounds_lat_2d(n,i,j);
+          }
+        }
+    }
+    else
+    {
+      boundsLonDest = domainDest_->bounds_lon_1d;
+      boundsLatDest = domainDest_->bounds_lat_1d;
+    }
+  }
+  else
+  {
+    bool isNorthPole = false;
+    bool isSouthPole = false;
+
+    CArray<double,1> lon_g ;
+    CArray<double,1> lat_g ;
+
+    if (!domainDest_->lonvalue_1d.isEmpty() && !domainDest_->latvalue_1d.isEmpty())
+    {
+      domainDest_->AllgatherRectilinearLonLat(domainDest_->lonvalue_1d,domainDest_->latvalue_1d, lon_g,lat_g) ;
+    }
+    else if (! domainDest_->latvalue_rectilinear_read_from_file.isEmpty() && ! domainDest_->lonvalue_rectilinear_read_from_file.isEmpty() )
+    {
+      lat_g=domainDest_->latvalue_rectilinear_read_from_file ;
+      lon_g=domainDest_->lonvalue_rectilinear_read_from_file ;
+    }
+    else if (!domainDest_->lon_start.isEmpty() && !domainDest_->lon_end.isEmpty() &&
+             !domainDest_->lat_start.isEmpty() && !domainDest_->lat_end.isEmpty())
+    {
+      double step=(domainDest_->lon_end-domainDest_->lon_start)/domainDest_->ni_glo ;
+      for(int i=0; i<domainDest_->ni_glo; ++i) lon_g(i)=domainDest_->lon_start+i*step ;
+      step=(domainDest_->lat_end-domainDest_->lat_start)/domainDest_->nj_glo ;
+      for(int i=0; i<domainDest_->ni_glo; ++i) lat_g(i)=domainDest_->lat_start+i*step ;
+    }
+    else ERROR("void CDomainAlgorithmInterpolate::computeRemap()",<<"Cannot compute bounds for rectilinear domain") ;
+    
+    if (std::abs(poleValue - std::abs(lat_g(0))) < NumTraits<double>::epsilon()) isNorthPole = true;
+    if (std::abs(poleValue - std::abs(lat_g(domainDest_->nj_glo-1))) < NumTraits<double>::epsilon()) isSouthPole = true;
+
+
+
+
+    if (isNorthPole && (0 == domainDest_->jbegin.getValue()))
+    {
+      int ibegin = domainDest_->ibegin.getValue();
+      for (i = 0; i < niDest; ++i)
+      {
+        interpMapValueNorthPole[i+ibegin];
+      }
+    }
+
+    if (isSouthPole && (domainDest_->nj_glo.getValue() == (domainDest_->jbegin.getValue() + njDest)))
+    {
+      int ibegin = domainDest_->ibegin.getValue();
+      int njGlo = domainDest_->nj_glo.getValue();
+      int niGlo = domainDest_->ni_glo.getValue();
+      for (i = 0; i < niDest; ++i)
+      {
+        k = (njGlo - 1)*niGlo + i + ibegin;
+        interpMapValueSouthPole[k];
+      }
+    }
+
+    // Ok, fill in boundary values for rectangular domain
+    nVertexDest = constNVertex;
+    domainDest_->fillInRectilinearBoundLonLat(lon_g,lat_g, boundsLonDest, boundsLatDest);
+  }
+
+
+
+  // Ok, now use mapper to calculate
+  int nSrcLocal = domainSrc_->i_index.numElements();
+  int nDstLocal = domainDest_->i_index.numElements();
+  long int * globalSrc = new long int [nSrcLocal];
+  long int * globalDst = new long int [nDstLocal];
+
+  long int globalIndex;
+  int i_ind, j_ind;
+  for (int idx = 0; idx < nSrcLocal; ++idx)
+  {
+    i_ind=domainSrc_->i_index(idx) ;
+    j_ind=domainSrc_->j_index(idx) ;
+
+    globalIndex = i_ind + j_ind * domainSrc_->ni_glo;
+    globalSrc[idx] = globalIndex;
+  }
+
+  for (int idx = 0; idx < nDstLocal; ++idx)
+  {
+    i_ind=domainDest_->i_index(idx) ;
+    j_ind=domainDest_->j_index(idx) ;
+
+    globalIndex = i_ind + j_ind * domainDest_->ni_glo;
+    globalDst[idx] = globalIndex;
+  }
+
+
+  // Calculate weight index
+  Mapper mapper(context->intraComm_);
+  mapper.setVerbosity(PROGRESS) ;
+
+     
+  // supress masked data for the source
+  int nSrcLocalUnmasked = 0 ;
+  for (int idx=0 ; idx < nSrcLocal; idx++) if (domainSrc_->localMask(idx)) ++nSrcLocalUnmasked ;
+
+
+  CArray<double,2> boundsLonSrcUnmasked(nVertexSrc,nSrcLocalUnmasked);
+  CArray<double,2> boundsLatSrcUnmasked(nVertexSrc,nSrcLocalUnmasked);
+  CArray<double,1> areaSrcUnmasked(nSrcLocalUnmasked);
+  
+  long int * globalSrcUnmasked = new long int [nSrcLocalUnmasked];
+
+  nSrcLocalUnmasked=0 ;
+  bool hasSrcArea=domainSrc_->hasArea && !domainSrc_->radius.isEmpty() && !interpDomain_->use_area.isEmpty() && interpDomain_->use_area==true  ;
+  double srcAreaFactor ;
+  if (hasSrcArea) srcAreaFactor=1./(domainSrc_->radius*domainSrc_->radius) ;
+  
+  for (int idx=0 ; idx < nSrcLocal; idx++)
+  {
+    if (domainSrc_->localMask(idx))
+    {
+      for(int n=0;n<nVertexSrc;++n)
+      {
+        boundsLonSrcUnmasked(n,nSrcLocalUnmasked) = boundsLonSrc(n,idx) ;
+        boundsLatSrcUnmasked(n,nSrcLocalUnmasked) = boundsLatSrc(n,idx) ;
+      }
+      if (hasSrcArea) areaSrcUnmasked(nSrcLocalUnmasked) = domainSrc_->areavalue(idx)*srcAreaFactor ;
+      globalSrcUnmasked[nSrcLocalUnmasked]=globalSrc[idx] ;
+      ++nSrcLocalUnmasked ;
+    }
+  }
+ 
+
+  int nDstLocalUnmasked = 0 ;
+  for (int idx=0 ; idx < nDstLocal; idx++) if (domainDest_->localMask(idx)) ++nDstLocalUnmasked ;
+
+  CArray<double,2> boundsLonDestUnmasked(nVertexDest,nDstLocalUnmasked);
+  CArray<double,2> boundsLatDestUnmasked(nVertexDest,nDstLocalUnmasked);
+  CArray<double,1>   areaDstUnmasked(nDstLocalUnmasked);
+
+  long int * globalDstUnmasked = new long int [nDstLocalUnmasked];
+
+  nDstLocalUnmasked=0 ;
+  bool hasDstArea=domainDest_->hasArea && !domainDest_->radius.isEmpty() && !interpDomain_->use_area.isEmpty() && interpDomain_->use_area==true ;
+  double dstAreaFactor ;
+  if (hasDstArea) dstAreaFactor=1./(domainDest_->radius*domainDest_->radius) ;
+  for (int idx=0 ; idx < nDstLocal; idx++)
+  {
+    if (domainDest_->localMask(idx))
+    {
+      for(int n=0;n<nVertexDest;++n)
+      {
+        boundsLonDestUnmasked(n,nDstLocalUnmasked) = boundsLonDest(n,idx) ;
+        boundsLatDestUnmasked(n,nDstLocalUnmasked) = boundsLatDest(n,idx) ;
+      }
+      if (hasDstArea) areaDstUnmasked(nDstLocalUnmasked) = domainDest_->areavalue(idx)*dstAreaFactor ;
+      globalDstUnmasked[nDstLocalUnmasked]=globalDst[idx] ;
+      ++nDstLocalUnmasked ;
+    }
+  }
+
+  double* ptAreaSrcUnmasked = NULL ;
+  if (hasSrcArea) ptAreaSrcUnmasked=areaSrcUnmasked.dataFirst() ;
+
+  double* ptAreaDstUnmasked = NULL ;
+  if (hasDstArea) ptAreaDstUnmasked=areaDstUnmasked.dataFirst() ;
+
+  mapper.setSourceMesh(boundsLonSrcUnmasked.dataFirst(), boundsLatSrcUnmasked.dataFirst(), ptAreaSrcUnmasked, nVertexSrc, nSrcLocalUnmasked, &srcPole[0], globalSrcUnmasked);
+  mapper.setTargetMesh(boundsLonDestUnmasked.dataFirst(), boundsLatDestUnmasked.dataFirst(), ptAreaDstUnmasked, nVertexDest, nDstLocalUnmasked, &dstPole[0], globalDstUnmasked);
+
+  std::vector<double> timings = mapper.computeWeights(orderInterp,renormalize,quantity);
+
+  std::map<int,std::vector<std::pair<int,double> > > interpMapValue;
+  std::map<int,std::vector<std::pair<int,double> > >::const_iterator iteNorthPole = interpMapValueNorthPole.end(),
+                                                                     iteSouthPole = interpMapValueSouthPole.end();
+  for (int idx = 0;  idx < mapper.nWeights; ++idx)
+  {
+    interpMapValue[mapper.targetWeightId[idx]].push_back(make_pair(mapper.sourceWeightId[idx],mapper.remapMatrix[idx]));
+    if (iteNorthPole != interpMapValueNorthPole.find(mapper.targetWeightId[idx]))
+    {
+      interpMapValueNorthPole[mapper.targetWeightId[idx]].push_back(make_pair(mapper.sourceWeightId[idx],mapper.remapMatrix[idx]));
+    }
+
+    if (iteSouthPole != interpMapValueSouthPole.find(mapper.targetWeightId[idx]))
+    {
+      interpMapValueSouthPole[mapper.targetWeightId[idx]].push_back(make_pair(mapper.sourceWeightId[idx],mapper.remapMatrix[idx]));
+    }
+  }
+  int niGloDst = domainDest_->ni_glo.getValue();
+  processPole(interpMapValueNorthPole, niGloDst);
+  processPole(interpMapValueSouthPole, niGloDst);
+
+  if (!interpMapValueNorthPole.empty())
+  {
+     std::map<int,std::vector<std::pair<int,double> > >::iterator itNorthPole = interpMapValueNorthPole.begin();
+     for (; itNorthPole != iteNorthPole; ++itNorthPole)
+     {
+       if (!(itNorthPole->second.empty()))
+        itNorthPole->second.swap(interpMapValue[itNorthPole->first]);
+     }
+  }
+
+  if (!interpMapValueSouthPole.empty())
+  {
+     std::map<int,std::vector<std::pair<int,double> > >::iterator itSouthPole = interpMapValueSouthPole.begin();
+     for (; itSouthPole != iteSouthPole; ++itSouthPole)
+     {
+       if (!(itSouthPole->second.empty()))
+        itSouthPole->second.swap(interpMapValue[itSouthPole->first]);
+     }
+  }
+
+  if (writeToFile_ && !readFromFile_) writeRemapInfo(interpMapValue);
+//  exchangeRemapInfo(interpMapValue);
+  convertRemapInfo(interpMapValue) ;
+
+  delete [] globalSrc;
+  delete [] globalSrcUnmasked;
+  delete [] globalDst;
+  delete [] globalDstUnmasked;
+
+}
+CATCH
+
+void CDomainAlgorithmInterpolate::processPole(std::map<int,std::vector<std::pair<int,double> > >& interMapValuePole,
+                                              int nbGlobalPointOnPole)
+TRY
+{
+  CContext* context = CContext::getCurrent();
+  
+  MPI_Comm poleComme(MPI_COMM_NULL);
+  MPI_Comm_split(context->intraComm_, interMapValuePole.empty() ? MPI_UNDEFINED : 1, 0, &poleComme);
+  if (MPI_COMM_NULL != poleComme)
+  {
+    int nbClientPole;
+    MPI_Comm_size(poleComme, &nbClientPole);
+
+    std::map<int,std::vector<std::pair<int,double> > >::iterator itePole = interMapValuePole.end(), itPole,
+                                                                 itbPole = interMapValuePole.begin();
+
+    int nbWeight = 0;
+    for (itPole = itbPole; itPole != itePole; ++itPole)
+       nbWeight += itPole->second.size();
+
+    std::vector<int> recvCount(nbClientPole,0);
+    std::vector<int> displ(nbClientPole,0);
+    MPI_Allgather(&nbWeight,1,MPI_INT,&recvCount[0],1,MPI_INT,poleComme) ;
+
+    displ[0]=0;
+    for(int n=1;n<nbClientPole;++n) displ[n]=displ[n-1]+recvCount[n-1] ;
+    int recvSize=displ[nbClientPole-1]+recvCount[nbClientPole-1] ;
+
+    std::vector<int> sendSourceIndexBuff(nbWeight);
+    std::vector<double> sendSourceWeightBuff(nbWeight);
+    int k = 0;
+    for (itPole = itbPole; itPole != itePole; ++itPole)
+    {
+      for (int idx = 0; idx < itPole->second.size(); ++idx)
+      {
+        sendSourceIndexBuff[k] = (itPole->second)[idx].first;
+        sendSourceWeightBuff[k] = (itPole->second)[idx].second;
+        ++k;
+      }
+    }
+
+    std::vector<int> recvSourceIndexBuff(recvSize);
+    std::vector<double> recvSourceWeightBuff(recvSize);
+
+    // Gather all index and weight for pole
+    MPI_Allgatherv(&sendSourceIndexBuff[0],nbWeight,MPI_INT,&recvSourceIndexBuff[0],&recvCount[0],&displ[0],MPI_INT,poleComme);
+    MPI_Allgatherv(&sendSourceWeightBuff[0],nbWeight,MPI_DOUBLE,&recvSourceWeightBuff[0],&recvCount[0],&displ[0],MPI_DOUBLE,poleComme);
+
+    std::map<int,double> recvTemp;
+    for (int idx = 0; idx < recvSize; ++idx)
+    {
+      if (recvTemp.end() != recvTemp.find(recvSourceIndexBuff[idx]))
+        recvTemp[recvSourceIndexBuff[idx]] += recvSourceWeightBuff[idx]/nbGlobalPointOnPole;
+      else
+        recvTemp[recvSourceIndexBuff[idx]] = recvSourceWeightBuff[idx]/nbGlobalPointOnPole;
+    }
+
+    std::map<int,double>::const_iterator itRecvTemp, itbRecvTemp = recvTemp.begin(), iteRecvTemp = recvTemp.end();
+
+    for (itPole = itbPole; itPole != itePole; ++itPole)
+    {
+      itPole->second.clear();
+      for (itRecvTemp = itbRecvTemp; itRecvTemp != iteRecvTemp; ++itRecvTemp)
+          itPole->second.push_back(make_pair(itRecvTemp->first, itRecvTemp->second));
+    }
+  }
+
+}
+CATCH
+
+/*!
+  Compute the index mapping between domain on grid source and one on grid destination
+*/
+void CDomainAlgorithmInterpolate::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+TRY
+{
+  if (readFromFile_)  
+    readRemapInfo();
+  else
+  {
+    computeRemap(); 
+  }
+}
+CATCH
+
+void CDomainAlgorithmInterpolate::writeRemapInfo(std::map<int,std::vector<std::pair<int,double> > >& interpMapValue)
+TRY
+{  
+  writeInterpolationInfo(fileToReadWrite_, interpMapValue);
+}
+CATCH
+
+void CDomainAlgorithmInterpolate::readRemapInfo()
+TRY
+{  
+  std::map<int,std::vector<std::pair<int,double> > > interpMapValue;
+  readInterpolationInfo(fileToReadWrite_, interpMapValue);
+
+  exchangeRemapInfo(interpMapValue);
+}
+CATCH
+
+void CDomainAlgorithmInterpolate::convertRemapInfo(std::map<int,std::vector<std::pair<int,double> > >& interpMapValue)
+TRY
+{
+  CContext* context = CContext::getCurrent();
+  int clientRank = context->intraCommRank_;
+
+  this->transformationMapping_.resize(1);
+  this->transformationWeight_.resize(1);
+
+  TransformationIndexMap& transMap = this->transformationMapping_[0];
+  TransformationWeightMap& transWeight = this->transformationWeight_[0];
+
+  std::map<int,std::vector<std::pair<int,double> > >::const_iterator itb = interpMapValue.begin(), it,
+                                                                     ite = interpMapValue.end();
+  
+  for (it = itb; it != ite; ++it)
+  {    
+    const std::vector<std::pair<int,double> >& tmp = it->second;
+    for (int i = 0; i < tmp.size(); ++i)
+    {
+      transMap[it->first].push_back(tmp[i].first);
+      transWeight[it->first].push_back(tmp[i].second);
+    }      
+  }      
+}
+CATCH
+
+/*!
+  Read remap information from file then distribute it among clients
+*/
+void CDomainAlgorithmInterpolate::exchangeRemapInfo(std::map<int,std::vector<std::pair<int,double> > >& interpMapValue)
+TRY
+{
+  CContext* context = CContext::getCurrent();
+  int clientRank = context->intraCommRank_;
+  int nbClient = context-> intraCommSize_;
+
+  this->transformationMapping_.resize(1);
+  this->transformationWeight_.resize(1);
+
+  TransformationIndexMap& transMap = this->transformationMapping_[0];
+  TransformationWeightMap& transWeight = this->transformationWeight_[0];
+
+  std::unordered_map<size_t,int> globalIndexOfDomainDest;
+  int ni = domainDest_->ni.getValue();
+  int nj = domainDest_->nj.getValue();
+  int ni_glo = domainDest_->ni_glo.getValue();
+  size_t globalIndex;
+  int nIndexSize = domainDest_->i_index.numElements(), i_ind, j_ind;
+  for (int idx = 0; idx < nIndexSize; ++idx)
+  {
+    i_ind=domainDest_->i_index(idx) ;
+    j_ind=domainDest_->j_index(idx) ;
+
+    globalIndex = i_ind + j_ind * ni_glo;
+    globalIndexOfDomainDest[globalIndex] = clientRank;
+  }
+
+  CClientServerMappingDistributed domainIndexClientClientMapping(globalIndexOfDomainDest,
+                                                                 context->intraComm_,
+                                                                 true);
+  CArray<size_t,1> globalIndexInterp(interpMapValue.size());
+  std::map<int,std::vector<std::pair<int,double> > >::const_iterator itb = interpMapValue.begin(), it,
+                                                                     ite = interpMapValue.end();
+  size_t globalIndexCount = 0;
+  for (it = itb; it != ite; ++it)
+  {
+    globalIndexInterp(globalIndexCount) = it->first;
+    ++globalIndexCount;
+  }
+
+  domainIndexClientClientMapping.computeServerIndexMapping(globalIndexInterp, nbClient);
+  const CClientServerMapping::GlobalIndexMap& globalIndexInterpSendToClient = domainIndexClientClientMapping.getGlobalIndexOnServer();
+
+  //Inform each client number of index they will receive
+  int* sendBuff = new int[nbClient];
+  int* recvBuff = new int[nbClient];
+  for (int i = 0; i < nbClient; ++i)
+  {
+    sendBuff[i] = 0;
+    recvBuff[i] = 0;
+  }
+  int sendBuffSize = 0;
+  CClientServerMapping::GlobalIndexMap::const_iterator itbMap = globalIndexInterpSendToClient.begin(), itMap,
+                                                       iteMap = globalIndexInterpSendToClient.end();
+  for (itMap = itbMap; itMap != iteMap; ++itMap)
+  {
+    const std::vector<size_t>& tmp = itMap->second;
+    int sizeIndex = 0, mapSize = (itMap->second).size();
+    for (int idx = 0; idx < mapSize; ++idx)
+    {
+//      sizeIndex += interpMapValue.at((itMap->second)[idx]).size();
+      sizeIndex += (interpMapValue[(int)(itMap->second)[idx]]).size();
+    }
+    sendBuff[itMap->first] = sizeIndex;
+    sendBuffSize += sizeIndex;
+  }
+
+
+  MPI_Allreduce(sendBuff, recvBuff, nbClient, MPI_INT, MPI_SUM, context->intraComm_);
+
+  int* sendIndexDestBuff = new int [sendBuffSize];
+  int* sendIndexSrcBuff  = new int [sendBuffSize];
+  double* sendWeightBuff = new double [sendBuffSize];
+
+  std::vector<MPI_Request> sendRequest;
+
+  int sendOffSet = 0, l = 0;
+  for (itMap = itbMap; itMap != iteMap; ++itMap)
+  {
+    const std::vector<size_t>& indexToSend = itMap->second;
+    int mapSize = indexToSend.size();
+    int k = 0;
+    for (int idx = 0; idx < mapSize; ++idx)
+    {
+      std::vector<std::pair<int,double> >& interpMap = interpMapValue[(int)indexToSend[idx]]; //interpMapValue.at(indexToSend[idx]);
+      for (int i = 0; i < interpMap.size(); ++i)
+      {
+        sendIndexDestBuff[l] = indexToSend[idx];
+        sendIndexSrcBuff[l]  = interpMap[i].first;
+        sendWeightBuff[l]    = interpMap[i].second;
+        ++k;
+        ++l;
+      }
+    }
+
+    sendRequest.push_back(MPI_Request());
+    MPI_Isend(sendIndexDestBuff + sendOffSet,
+             k,
+             MPI_INT,
+             itMap->first,
+             MPI_DOMAIN_INTERPOLATION_DEST_INDEX,
+             context->intraComm_,
+             &sendRequest.back());
+    sendRequest.push_back(MPI_Request());
+    MPI_Isend(sendIndexSrcBuff + sendOffSet,
+             k,
+             MPI_INT,
+             itMap->first,
+             MPI_DOMAIN_INTERPOLATION_SRC_INDEX,
+             context->intraComm_,
+             &sendRequest.back());
+    sendRequest.push_back(MPI_Request());
+    MPI_Isend(sendWeightBuff + sendOffSet,
+             k,
+             MPI_DOUBLE,
+             itMap->first,
+             MPI_DOMAIN_INTERPOLATION_WEIGHT,
+             context->intraComm_,
+             &sendRequest.back());
+    sendOffSet += k;
+  }
+
+  int recvBuffSize = recvBuff[clientRank];
+  int* recvIndexDestBuff = new int [recvBuffSize];
+  int* recvIndexSrcBuff  = new int [recvBuffSize];
+  double* recvWeightBuff = new double [recvBuffSize];
+  int receivedSize = 0;
+  int clientSrcRank;
+  while (receivedSize < recvBuffSize)
+  {
+    MPI_Status recvStatus;
+    MPI_Recv((recvIndexDestBuff + receivedSize),
+             recvBuffSize,
+             MPI_INT,
+             MPI_ANY_SOURCE,
+             MPI_DOMAIN_INTERPOLATION_DEST_INDEX,
+             context->intraComm_,
+             &recvStatus);
+
+    int countBuff = 0;
+    MPI_Get_count(&recvStatus, MPI_INT, &countBuff);
+    clientSrcRank = recvStatus.MPI_SOURCE;
+
+    MPI_Recv((recvIndexSrcBuff + receivedSize),
+             recvBuffSize,
+             MPI_INT,
+             clientSrcRank,
+             MPI_DOMAIN_INTERPOLATION_SRC_INDEX,
+             context->intraComm_,
+             &recvStatus);
+
+    MPI_Recv((recvWeightBuff + receivedSize),
+             recvBuffSize,
+             MPI_DOUBLE,
+             clientSrcRank,
+             MPI_DOMAIN_INTERPOLATION_WEIGHT,
+             context->intraComm_,
+             &recvStatus);
+
+    for (int idx = 0; idx < countBuff; ++idx)
+    {
+      transMap[*(recvIndexDestBuff + receivedSize + idx)].push_back(*(recvIndexSrcBuff + receivedSize + idx));
+      transWeight[*(recvIndexDestBuff + receivedSize + idx)].push_back(*(recvWeightBuff + receivedSize + idx));
+    }
+    receivedSize += countBuff;
+  }
+
+  std::vector<MPI_Status> requestStatus(sendRequest.size());
+  MPI_Waitall(sendRequest.size(), &sendRequest[0], MPI_STATUS_IGNORE);
+
+  delete [] sendIndexDestBuff;
+  delete [] sendIndexSrcBuff;
+  delete [] sendWeightBuff;
+  delete [] recvIndexDestBuff;
+  delete [] recvIndexSrcBuff;
+  delete [] recvWeightBuff;
+  delete [] sendBuff;
+  delete [] recvBuff;
+}
+CATCH
+ 
+/*! Redefined some functions of CONetCDF4 to make use of them */
+CDomainAlgorithmInterpolate::WriteNetCdf::WriteNetCdf(const StdString& filename, const MPI_Comm comm)
+  : CNc4DataOutput(NULL, filename, false, false, true, comm, false, true) {}
+int CDomainAlgorithmInterpolate::WriteNetCdf::addDimensionWrite(const StdString& name, 
+                                                                const StdSize size)
+TRY
+{
+  return CONetCDF4::addDimension(name, size);  
+}
+CATCH
+
+int CDomainAlgorithmInterpolate::WriteNetCdf::addVariableWrite(const StdString& name, nc_type type,
+                                                               const std::vector<StdString>& dim)
+TRY
+{
+  return CONetCDF4::addVariable(name, type, dim);
+}
+CATCH
+
+void CDomainAlgorithmInterpolate::WriteNetCdf::endDefinition()
+TRY
+{
+  CONetCDF4::definition_end();
+}
+CATCH
+
+void CDomainAlgorithmInterpolate::WriteNetCdf::writeDataIndex(const CArray<int,1>& data, const StdString& name,
+                                                              bool collective, StdSize record,
+                                                              const std::vector<StdSize>* start,
+                                                              const std::vector<StdSize>* count)
+TRY
+{
+  CONetCDF4::writeData<int,1>(data, name, collective, record, start, count);
+}
+CATCH
+
+void CDomainAlgorithmInterpolate::WriteNetCdf::writeDataIndex(const CArray<double,1>& data, const StdString& name,
+                                                              bool collective, StdSize record,
+                                                              const std::vector<StdSize>* start,
+                                                              const std::vector<StdSize>* count)
+TRY
+{
+  CONetCDF4::writeData<double,1>(data, name, collective, record, start, count);
+}
+CATCH
+
+/*
+   Write interpolation weights into a file
+   \param [in] filename name of output file
+   \param interpMapValue mapping of global index of domain destination and domain source as well as the corresponding weight
+*/
+void CDomainAlgorithmInterpolate::writeInterpolationInfo(std::string& filename,
+                                                         std::map<int,std::vector<std::pair<int,double> > >& interpMapValue)
+TRY
+{
+  CContext* context = CContext::getCurrent();
+  
+  size_t n_src = domainSrc_->ni_glo * domainSrc_->nj_glo;
+  size_t n_dst = domainDest_->ni_glo * domainDest_->nj_glo;
+
+  long localNbWeight = 0;
+  long globalNbWeight;
+  long startIndex;
+  typedef std::map<int,std::vector<std::pair<int,double> > > IndexRemap;
+  IndexRemap::iterator itb = interpMapValue.begin(), it,
+                       ite = interpMapValue.end();
+  for (it = itb; it!=ite; ++it)
+  {
+    localNbWeight += (it->second).size();
+  }
+
+  CArray<int,1> src_idx(localNbWeight);
+  CArray<int,1> dst_idx(localNbWeight);
+  CArray<double,1> weights(localNbWeight);
+
+  int index = 0;
+  int indexOffset=0 ;
+  if (fortranConvention) indexOffset=1 ;
+  for (it = itb; it !=ite; ++it)
+  {
+    std::vector<std::pair<int,double> >& tmp = it->second;
+    for (int idx = 0; idx < tmp.size(); ++idx)
+    {
+      dst_idx(index) = it->first + indexOffset;
+      src_idx(index) = tmp[idx].first + indexOffset;
+      weights(index) = tmp[idx].second;
+      ++index;
+    }    
+  }
+
+  MPI_Allreduce(&localNbWeight, &globalNbWeight, 1, MPI_LONG, MPI_SUM, context->intraComm_);
+  MPI_Scan(&localNbWeight, &startIndex, 1, MPI_LONG, MPI_SUM, context->intraComm_);
+  
+  if (0 == globalNbWeight)
+  {
+    info << "There is no interpolation weights calculated between "
+         << "domain source: " << domainSrc_->getDomainOutputName()
+         << " and domain destination: " << domainDest_->getDomainOutputName()
+         << std::endl;
+    return;
+  }
+
+  std::vector<StdSize> start(1, startIndex - localNbWeight);
+  std::vector<StdSize> count(1, localNbWeight);
+  
+  WriteNetCdf netCdfWriter(filename, context->intraComm_);  
+
+  // Define some dimensions
+  netCdfWriter.addDimensionWrite("n_src", n_src);
+  netCdfWriter.addDimensionWrite("n_dst", n_dst);
+  netCdfWriter.addDimensionWrite("n_weight", globalNbWeight);
+  
+  std::vector<StdString> dims(1,"n_weight");
+
+  // Add some variables
+  netCdfWriter.addVariableWrite("src_idx", NC_INT, dims);
+  netCdfWriter.addVariableWrite("dst_idx", NC_INT, dims);
+  netCdfWriter.addVariableWrite("weight", NC_DOUBLE, dims);
+
+  // End of definition
+  netCdfWriter.endDefinition();
+
+  // // Write variables
+  if (0 != localNbWeight)
+  {
+    netCdfWriter.writeDataIndex(src_idx, "src_idx", false, 0, &start, &count);
+    netCdfWriter.writeDataIndex(dst_idx, "dst_idx", false, 0, &start, &count);
+    netCdfWriter.writeDataIndex(weights, "weight", false, 0, &start, &count);
+  }
+
+  netCdfWriter.closeFile();
+}
+CATCH
+
+/*!
+  Read interpolation information from a file
+  \param [in] filename interpolation file
+  \param [in/out] interpMapValue Mapping between (global) index of domain on grid destination and
+         corresponding global index of domain and associated weight value on grid source
+*/
+void CDomainAlgorithmInterpolate::readInterpolationInfo(std::string& filename,
+                                                        std::map<int,std::vector<std::pair<int,double> > >& interpMapValue)
+TRY
+{
+  int ncid ;
+  int weightDimId ;
+  size_t nbWeightGlo ;
+
+
+  CContext* context = CContext::getCurrent();
+  int clientRank = context->intraCommRank_;
+  int clientSize = context->intraCommSize_;
+
+
+  {
+    ifstream f(filename.c_str());
+    if (!f.good()) ERROR("void CDomainAlgorithmInterpolate::readInterpolationInfo",
+                      << "Attempt to read file weight :"  << filename << " which doesn't seem to exist." << std::endl
+                      << "Please check this file ");
+  }
+                  
+  nc_open(filename.c_str(),NC_NOWRITE, &ncid) ;
+  nc_inq_dimid(ncid,"n_weight",&weightDimId) ;
+  nc_inq_dimlen(ncid,weightDimId,&nbWeightGlo) ;
+
+  size_t nbWeight ;
+  size_t start ;
+  size_t div = nbWeightGlo/clientSize ;
+  size_t mod = nbWeightGlo%clientSize ;
+  if (clientRank < mod)
+  {
+    nbWeight=div+1 ;
+    start=clientRank*(div+1) ;
+  }
+  else
+  {
+    nbWeight=div ;
+    start= mod * (div+1) + (clientRank-mod) * div ;
+  }
+
+  double* weight=new double[nbWeight] ;
+  int weightId ;
+  nc_inq_varid (ncid, "weight", &weightId) ;
+  nc_get_vara_double(ncid, weightId, &start, &nbWeight, weight) ;
+
+  long* srcIndex=new long[nbWeight] ;
+  int srcIndexId ;
+  nc_inq_varid (ncid, "src_idx", &srcIndexId) ;
+  nc_get_vara_long(ncid, srcIndexId, &start, &nbWeight, srcIndex) ;
+
+  long* dstIndex=new long[nbWeight] ;
+  int dstIndexId ;
+  nc_inq_varid (ncid, "dst_idx", &dstIndexId) ;
+  nc_get_vara_long(ncid, dstIndexId, &start, &nbWeight, dstIndex) ;
+
+  int indexOffset=0 ;
+  if (fortranConvention) indexOffset=1 ;
+    for(size_t ind=0; ind<nbWeight;++ind)
+      interpMapValue[dstIndex[ind]-indexOffset].push_back(make_pair(srcIndex[ind]-indexOffset,weight[ind]));
+ }
+CATCH
+
+void CDomainAlgorithmInterpolate::apply(const std::vector<std::pair<int,double> >& localIndex,
+                                            const double* dataInput,
+                                            CArray<double,1>& dataOut,
+                                            std::vector<bool>& flagInitial,
+                                            bool ignoreMissingValue, bool firstPass  )
+TRY
+{
+  int nbLocalIndex = localIndex.size();   
+  double defaultValue = std::numeric_limits<double>::quiet_NaN();
+   
+  if (detectMissingValue)
+  {
+     if (firstPass && renormalize)
+     {
+       renormalizationFactor.resize(dataOut.numElements()) ;
+       renormalizationFactor=1 ;
+     }
+
+    if (firstPass)
+    {
+      allMissing.resize(dataOut.numElements()) ;
+      allMissing=true ;
+    }
+
+    for (int idx = 0; idx < nbLocalIndex; ++idx)
+    {
+      if (NumTraits<double>::isNan(*(dataInput + idx)))
+      {
+        allMissing(localIndex[idx].first) = allMissing(localIndex[idx].first) && true;
+        if (renormalize) renormalizationFactor(localIndex[idx].first)-=localIndex[idx].second ;
+      }
+      else
+      {
+        dataOut(localIndex[idx].first) += *(dataInput + idx) * localIndex[idx].second;
+        allMissing(localIndex[idx].first) = allMissing(localIndex[idx].first) && false; // Reset flag to indicate not all data source are nan
+      }
+    }
+
+  }
+  else
+  {
+    for (int idx = 0; idx < nbLocalIndex; ++idx)
+    {
+      dataOut(localIndex[idx].first) += *(dataInput + idx) * localIndex[idx].second;
+    }
+  }
+}
+CATCH
+
+void CDomainAlgorithmInterpolate::updateData(CArray<double,1>& dataOut)
+TRY
+{
+  if (detectMissingValue)
+  {
+    double defaultValue = std::numeric_limits<double>::quiet_NaN();
+    size_t nbIndex=dataOut.numElements() ; 
+
+    if (allMissing.numElements()>0)
+    {
+      for (int idx = 0; idx < nbIndex; ++idx)
+      {
+        if (allMissing(idx)) dataOut(idx) = defaultValue; // If all data source are nan then data destination must be nan
+      }
+    }
+    
+    if (renormalize)
+    {
+      if (renormalizationFactor.numElements()>0) dataOut/=renormalizationFactor ; // In some case, process doesn't received any data for interpolation (mask)
+                                                                                 // so renormalizationFactor is not initialized
+    }
+  }
+}
+CATCH
+
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_interpolate.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_interpolate.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_interpolate.hpp	(revision 1985)
@@ -0,0 +1,103 @@
+/*!
+   \file domain_algorithm_interpolate_from_file.hpp
+   \author Ha NGUYEN
+   \since 09 July 2015
+   \date 09 Sep 2015
+
+   \brief Algorithm for interpolation on a domain.
+ */
+#ifndef __XIOS_DOMAIN_ALGORITHM_INTERPOLATE_HPP__
+#define __XIOS_DOMAIN_ALGORITHM_INTERPOLATE_HPP__
+
+#include "domain_algorithm_transformation.hpp"
+#include "transformation.hpp"
+#include "nc4_data_output.hpp"
+
+namespace xios {
+
+class CDomain;
+class CInterpolateDomain;
+
+
+/*!
+  \class CDomainAlgorithmInterpolate
+  Reading interpolation from file then apply on a domain
+*/
+class CDomainAlgorithmInterpolate : public CDomainAlgorithmTransformation
+{
+public:
+  CDomainAlgorithmInterpolate(CDomain* domainDestination, CDomain* domainSource, CInterpolateDomain* interpDomain);
+
+  virtual ~CDomainAlgorithmInterpolate() {}
+
+  virtual void apply(const std::vector<std::pair<int,double> >& localIndex,
+                     const double* dataInput,
+                     CArray<double,1>& dataOut,
+                     std::vector<bool>& flagInitial,                     
+                     bool ignoreMissingValue, bool firstPass);
+  virtual void updateData(CArray<double,1>& dataOut);
+
+  static bool registerTrans();
+protected:
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs);
+
+private:
+  void readInterpolationInfo(std::string& filename, std::map<int,std::vector<std::pair<int,double> > >& interpMapValue);
+  void writeInterpolationInfo(std::string& filename, std::map<int,std::vector<std::pair<int,double> > >& interpMapValue);
+  void processPole(std::map<int,std::vector<std::pair<int,double> > >& interMapValuePole,
+                   int nbGlobalPointOnPole);
+  void computeRemap();
+  void readRemapInfo();
+  void writeRemapInfo(std::map<int,std::vector<std::pair<int,double> > >&);
+  void exchangeRemapInfo(std::map<int,std::vector<std::pair<int,double> > >& interpMapValue);
+  void convertRemapInfo(std::map<int,std::vector<std::pair<int,double> > >& interpMapValue);
+
+private:
+  CArray<double,1> renormalizationFactor ;
+  CArray<bool,1> allMissing ;
+  bool detectMissingValue ;
+  bool renormalize ;
+  bool quantity ;
+  
+  CInterpolateDomain* interpDomain_;
+  bool writeToFile_;
+  bool readFromFile_;
+  StdString fileToReadWrite_;
+  bool fortranConvention ;
+
+  // class WriteNetCdf;
+  class WriteNetCdf : public CNc4DataOutput
+  {
+  public:
+    WriteNetCdf(const StdString& filename, const MPI_Comm comm);
+    int addDimensionWrite(const StdString& name, const StdSize size = UNLIMITED_DIM);
+    int addVariableWrite(const StdString& name, nc_type type,
+                         const std::vector<StdString>& dim);
+    void endDefinition();
+    void writeDataIndex(const CArray<int,1>& data, const StdString& name,
+                        bool collective, StdSize record,
+                        const std::vector<StdSize>* start = NULL,
+                        const std::vector<StdSize>* count = NULL);
+    void writeDataIndex(const CArray<double,1>& data, const StdString& name,
+                        bool collective, StdSize record,
+                        const std::vector<StdSize>* start = NULL,
+                        const std::vector<StdSize>* count = NULL);
+  };
+
+
+private:
+
+  static CGenericAlgorithmTransformation* create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                CTransformation<CDomain>* transformation,
+                                                int elementPositionInGrid,
+                                                std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridDst2DomainPosition);
+  static bool dummyRegistered_;
+};
+
+}
+#endif // __XIOS_DOMAIN_ALGORITHM_INTERPOLATE_HPP__
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_reorder.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_reorder.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_reorder.cpp	(revision 1985)
@@ -0,0 +1,135 @@
+/*!
+   \file domain_algorithm_reorder.cpp
+   \brief Algorithm for reorder a domain.
+ */
+#include "domain_algorithm_reorder.hpp"
+#include "reorder_domain.hpp"
+#include "domain.hpp"
+#include "grid.hpp"
+#include "grid_transformation_factory_impl.hpp"
+
+namespace xios {
+CGenericAlgorithmTransformation* CDomainAlgorithmReorder::create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                             CTransformation<CDomain>* transformation,
+                                                             int elementPositionInGrid,
+                                                             std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                             std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                             std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                             std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                             std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                             std::map<int, int>& elementPositionInGridDst2DomainPosition)
+TRY
+{
+  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
+  std::vector<CDomain*> domainListSrcP  = gridSrc->getDomains();
+
+  CReorderDomain* reorderDomain = dynamic_cast<CReorderDomain*> (transformation);
+  int domainDstIndex = elementPositionInGridDst2DomainPosition[elementPositionInGrid];
+  int domainSrcIndex = elementPositionInGridSrc2DomainPosition[elementPositionInGrid];
+
+  return (new CDomainAlgorithmReorder(domainListDestP[domainDstIndex], domainListSrcP[domainSrcIndex], reorderDomain));
+}
+CATCH
+
+bool CDomainAlgorithmReorder::dummyRegistered_ = CDomainAlgorithmReorder::registerTrans();
+bool CDomainAlgorithmReorder::registerTrans()
+TRY
+{
+  return CGridTransformationFactory<CDomain>::registerTransformation(TRANS_REORDER_DOMAIN, create);
+}
+CATCH
+
+CDomainAlgorithmReorder::CDomainAlgorithmReorder(CDomain* domainDestination, CDomain* domainSource, CReorderDomain* reorderDomain)
+: CDomainAlgorithmTransformation(domainDestination, domainSource)
+TRY
+{
+  reorderDomain->checkValid(domainSource);
+  if (domainDestination->type !=  CDomain::type_attr::rectilinear)
+  {
+      ERROR("CDomainAlgorithmReorder::CDomainAlgorithmReorder(CDomain* domainDestination, CDomain* domainSource, CReorderDomain* reorderDomain)",
+           << "Domain destination is not rectilinear. This filter work only for rectilinear domain and destination domain with < id = "
+           <<domainDestination->getId() <<" > is of type "<<domainDestination->type<<std::endl);
+  }
+  
+  if (domainDestination == domainSource)
+  {
+    ERROR("CDomainAlgorithmReorder::CDomainAlgorithmReorder(CDomain* domainDestination, CDomain* domainSource, CReorderDomain* reorderDomain)",
+           << "Domain source and domain destination are the same. Please make sure domain destination refers to domain source" << std::endl
+           << "Domain source " <<domainSource->getId() << std::endl
+           << "Domain destination " <<domainDestination->getId() << std::endl);
+  }
+  this->type_ = (ELEMENT_MODIFICATION_WITHOUT_DATA);
+
+  if (!reorderDomain->invert_lat.isEmpty())
+  {
+    CArray<int,1>& j_index=domainDestination->j_index ;
+    int nglo = j_index.numElements() ;
+    int nj_glo =domainDestination->nj_glo ;
+    
+    for (size_t i = 0; i < nglo ; ++i)
+    {
+      j_index(i)=(nj_glo-1)-j_index(i) ;
+    }
+  }
+
+  if (!reorderDomain->shift_lon_fraction.isEmpty())
+  {
+    int ni_glo =domainDestination->ni_glo ;
+    int  offset = ni_glo*reorderDomain->shift_lon_fraction ;
+    CArray<int,1>& i_index=domainDestination->i_index ;
+    int nglo = i_index.numElements() ;
+    
+    for (size_t i = 0; i < nglo ; ++i)
+    {
+      i_index(i)=  (i_index(i)+offset+ni_glo)%ni_glo ;
+    }
+  }
+
+  if (!reorderDomain->min_lon.isEmpty() && !reorderDomain->max_lon.isEmpty())
+  {
+    double min_lon=reorderDomain->min_lon ;
+    double max_lon=reorderDomain->max_lon ;
+    double delta=max_lon-min_lon ;
+    
+    if (!domainDestination->lonvalue_1d.isEmpty() )
+    {
+      CArray<double,1>& lon=domainDestination->lonvalue_1d ;
+      for (int i=0;i<lon.numElements();++i)
+      {
+        while  (lon(i) > max_lon) lon(i)=lon(i)-delta ;
+        while  (lon(i) < min_lon) lon(i)=lon(i)+delta ;
+      }
+    }
+
+    if (!domainDestination->bounds_lon_1d.isEmpty() )
+    {
+      CArray<double,2>& bounds_lon=domainDestination->bounds_lon_1d ;
+      for (int i=0;i<bounds_lon.extent(0);++i)
+      {
+        while  (bounds_lon(0,i) > max_lon) bounds_lon(0,i)=bounds_lon(0,i)-delta ;
+        while  (bounds_lon(1,i) > max_lon) bounds_lon(1,i)=bounds_lon(1,i)-delta ;
+
+        while  (bounds_lon(0,i) < min_lon) bounds_lon(0,i)=bounds_lon(0,i)+delta ;
+        while  (bounds_lon(1,i) < min_lon) bounds_lon(1,i)=bounds_lon(1,i)+delta ;
+      }
+    }
+  }
+}
+CATCH
+
+/*!
+  Compute the index mapping between domain on grid source and one on grid destination
+*/
+void CDomainAlgorithmReorder::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+{
+/*
+  this->transformationMapping_.resize(1);
+  this->transformationWeight_.resize(1);
+
+  TransformationIndexMap& transMap = this->transformationMapping_[0];
+  TransformationWeightMap& transWeight = this->transformationWeight_[0];
+*/
+}
+
+
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_reorder.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_reorder.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_reorder.hpp	(revision 1985)
@@ -0,0 +1,47 @@
+/*!
+   \brief Algorithm for reordering domain.
+ */
+#ifndef __XIOS_DOMAIN_ALGORITHM_REORDER_HPP__
+#define __XIOS_DOMAIN_ALGORITHM_REORDER_HPP__
+
+#include "domain_algorithm_transformation.hpp"
+#include "transformation.hpp"
+
+namespace xios {
+
+class CDomain;
+class CReorderDomain;
+
+/*!
+  \class CDomainAlgorithmReorder
+  Reordering data on domain
+*/
+class CDomainAlgorithmReorder : public CDomainAlgorithmTransformation
+{
+public:
+  CDomainAlgorithmReorder(CDomain* domainDestination, CDomain* domainSource, CReorderDomain* reorderDomain);
+
+  virtual ~CDomainAlgorithmReorder() {}
+
+  static bool registerTrans();
+
+protected:
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs);
+
+
+private:
+
+  static CGenericAlgorithmTransformation* create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                CTransformation<CDomain>* transformation,
+                                                int elementPositionInGrid,
+                                                std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridDst2DomainPosition);
+  static bool dummyRegistered_;
+};
+
+}
+#endif // __XIOS_DOMAIN_ALGORITHM_REORDER_HPP__
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_transformation.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_transformation.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_transformation.cpp	(revision 1985)
@@ -0,0 +1,69 @@
+/*!
+   \file domain_algorithm_transformation.hpp
+   \author Ha NGUYEN
+   \since 02 Jul 2015
+   \date 02 Jul 2015
+
+   \brief Interface for all domain transformation algorithms.
+ */
+
+#include "domain_algorithm_transformation.hpp"
+#include "context.hpp"
+#include "context_client.hpp"
+#include "domain.hpp"
+
+namespace xios {
+
+CDomainAlgorithmTransformation::CDomainAlgorithmTransformation(CDomain* domainDestination, CDomain* domainSource)
+ : CGenericAlgorithmTransformation(), domainDest_(domainDestination), domainSrc_(domainSource)
+{
+}
+
+CDomainAlgorithmTransformation::~CDomainAlgorithmTransformation()
+{
+}
+
+void CDomainAlgorithmTransformation::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+{
+}
+
+/*!
+  Compute global index of domain on different processes
+  \param [in] globalDomainIndex global index of domain source
+  \param [out] globalDomainIndexOnProc processes which contain the corresponding global index of domain source
+*/
+void CDomainAlgorithmTransformation::computeExchangeGlobalIndex(const CArray<size_t,1>& globalDomainIndex,
+                                                                int elementType,
+                                                                CClientClientDHTInt::Index2VectorInfoTypeMap& globalDomainIndexOnProc)
+TRY
+{
+  CContext* context = CContext::getCurrent();
+  int clientRank = context->intraCommRank_;
+  int clientSize = context->intraCommSize_;
+
+  int niGlob = domainSrc_->ni_glo.getValue();
+  int njGlob = domainSrc_->nj_glo.getValue();
+  size_t globalIndex;
+  int nIndexSize = domainSrc_->i_index.numElements(), i_ind, j_ind;
+  CClientClientDHTInt::Index2VectorInfoTypeMap globalIndex2ProcRank;
+  globalIndex2ProcRank.rehash(std::ceil(nIndexSize/globalIndex2ProcRank.max_load_factor()));
+  CArray<bool,1>& localMask=domainSrc_->localMask ;
+  for (int idx = 0; idx < nIndexSize; ++idx)
+  {
+    if (localMask(idx))
+    {
+      i_ind=domainSrc_->i_index(idx) ;
+      j_ind=domainSrc_->j_index(idx) ;
+
+      globalIndex = i_ind + j_ind * niGlob;
+      globalIndex2ProcRank[globalIndex].push_back(clientRank);
+    }
+  }
+
+  CClientClientDHTInt dhtIndexProcRank(globalIndex2ProcRank, context->intraComm_);
+  dhtIndexProcRank.computeIndexInfoMapping(globalDomainIndex);
+  globalDomainIndexOnProc = dhtIndexProcRank.getInfoIndexMap();
+}
+CATCH
+
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_transformation.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_transformation.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_transformation.hpp	(revision 1985)
@@ -0,0 +1,45 @@
+/*!
+   \file domain_algorithm_transformation.hpp
+   \author Ha NGUYEN
+   \since 02 Jul 2015
+   \date 02 Jul 2015
+
+   \brief Interface for all domain transformation algorithms.
+ */
+#ifndef __XIOS_DOMAIN_ALGORITHM_TRANSFORMATION_HPP__
+#define __XIOS_DOMAIN_ALGORITHM_TRANSFORMATION_HPP__
+
+#include "generic_algorithm_transformation.hpp"
+#include "client_client_dht_template.hpp"
+
+namespace xios {
+
+class CDomain;
+/*!
+  \class CDomainAlgorithmTransformation
+  Algorithms for domain.
+*/
+class CDomainAlgorithmTransformation : public virtual CGenericAlgorithmTransformation
+{
+public:
+  CDomainAlgorithmTransformation(CDomain* domainDestination, CDomain* domainSource);
+
+  virtual ~CDomainAlgorithmTransformation();
+
+protected:
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >&);
+
+  virtual void computeExchangeGlobalIndex(const CArray<size_t,1>& globalDomainIndex,
+                                          int elementType,
+                                          CClientClientDHTInt::Index2VectorInfoTypeMap& globalDomainIndexOnProc);
+
+protected:
+    //! Domain on grid destination
+  CDomain* domainDest_;
+
+  //! Domain on grid source
+  CDomain* domainSrc_;
+};
+
+}
+#endif // __XIOS_DOMAIN_ALGORITHM_TRANSFORMATION_HPP__
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_zoom.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_zoom.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_zoom.cpp	(revision 1985)
@@ -0,0 +1,290 @@
+#include "domain_algorithm_zoom.hpp"
+#include "zoom_domain.hpp"
+#include "domain.hpp"
+#include "grid.hpp"
+#include "grid_transformation_factory_impl.hpp"
+#include "attribute_template.hpp"
+
+namespace xios {
+CGenericAlgorithmTransformation* CDomainAlgorithmZoom::create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                             CTransformation<CDomain>* transformation,
+                                                             int elementPositionInGrid,
+                                                             std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                             std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                             std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                             std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                             std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                             std::map<int, int>& elementPositionInGridDst2DomainPosition)
+TRY
+{
+  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
+  std::vector<CDomain*> domainListSrcP  = gridSrc->getDomains();
+
+  CZoomDomain* zoomDomain = dynamic_cast<CZoomDomain*> (transformation);
+  int domainDstIndex = elementPositionInGridDst2DomainPosition[elementPositionInGrid];
+  int domainSrcIndex = elementPositionInGridSrc2DomainPosition[elementPositionInGrid];
+
+  return (new CDomainAlgorithmZoom(domainListDestP[domainDstIndex], domainListSrcP[domainSrcIndex], zoomDomain));
+}
+CATCH
+
+bool CDomainAlgorithmZoom::dummyRegistered_ = CDomainAlgorithmZoom::registerTrans();
+bool CDomainAlgorithmZoom::registerTrans()
+TRY
+{
+  return CGridTransformationFactory<CDomain>::registerTransformation(TRANS_ZOOM_DOMAIN, create);
+}
+CATCH
+
+CDomainAlgorithmZoom::CDomainAlgorithmZoom(CDomain* domainDestination, CDomain* domainSource, CZoomDomain* zoomDomain)
+: CDomainAlgorithmTransformation(domainDestination, domainSource)
+TRY
+{
+  zoomDomain->checkValid(domainSource);
+  zoomIBegin_ = zoomDomain->ibegin.getValue();
+  zoomJBegin_ = zoomDomain->jbegin.getValue();
+
+  zoomNi_  = zoomDomain->ni.getValue();
+  zoomNj_  = zoomDomain->nj.getValue();
+
+  zoomIEnd_ = zoomIBegin_ + zoomNi_ - 1;
+  zoomJEnd_ = zoomJBegin_ + zoomNj_ - 1;
+
+  if (zoomNi_ > domainSource->ni_glo.getValue())
+  {
+    ERROR("CDomainAlgorithmZoom::CDomainAlgorithmZoom(CDomain* domainDestination, CDomain* domainSource, CZoomDomain* zoomDomain)",
+           << "Zoom size is greater than size of domain source"
+           << "Size ni_glo of domain source " <<domainSource->getId() << " is " << domainSource->ni_glo.getValue()  << std::endl
+           << "Zoom size is " << zoomNi_ );
+  }
+
+  if (zoomNj_ > domainSource->nj_glo.getValue())
+  {
+    ERROR("CDomainAlgorithmZoom::CDomainAlgorithmZoom(CDomain* domainDestination, CDomain* domainSource, CZoomDomain* zoomDomain)",
+           << "Zoom size is greater than size of domain source"
+           << "Size nj_glo of domain source " <<domainSource->getId() << " is " << domainSource->nj_glo.getValue()  << std::endl
+           << "Zoom size is " << zoomNj_ );
+  }
+
+  // Calculate the size of local domain
+  int ind, indLocSrc, indLocDest, iIdxSrc, jIdxSrc, destIBegin = -1, destJBegin = -1, niDest = 0, njDest = 0, ibeginDest, jbeginDest ;
+  int indGloDest, indGloSrc, niGloSrc = domainSrc_->ni_glo, iSrc, jSrc, nvertex = 0;
+  for (int j = 0; j < domainSrc_->nj.getValue(); j++)
+  {
+    for (int i = 0; i < domainSrc_->ni.getValue(); i++)
+    {
+      ind = j*domainSrc_->ni + i;
+      iIdxSrc = domainSrc_->i_index(ind);
+      if ((iIdxSrc >= zoomIBegin_) && (iIdxSrc <= zoomIEnd_))
+      {
+        jIdxSrc = domainSrc_->j_index(ind);
+        if ((jIdxSrc >= zoomJBegin_) && (jIdxSrc <= zoomJEnd_))
+        {
+          if ((niDest == 0) && (njDest == 0))
+          {
+            destIBegin = i;
+            destJBegin = j;
+          }
+          if (i == destIBegin) ++njDest;
+        }
+        if (j == destJBegin) ++niDest;
+
+      }
+    }
+  }
+  ibeginDest = destIBegin + domainSrc_->ibegin - zoomIBegin_;
+  jbeginDest = destJBegin + domainSrc_->jbegin - zoomJBegin_;
+  domainDest_->ni_glo.setValue(zoomNi_);
+  domainDest_->nj_glo.setValue(zoomNj_);
+  domainDest_->ni.setValue(niDest);
+  domainDest_->nj.setValue(njDest);
+  if ( (niDest==0) || (njDest==0))
+  {
+    domainDest_->ibegin.setValue(0);
+    domainDest_->jbegin.setValue(0);
+  }
+  else
+  {
+    domainDest_->ibegin.setValue(ibeginDest);
+    domainDest_->jbegin.setValue(jbeginDest);
+  }
+  domainDest_->i_index.resize(niDest*njDest);
+  domainDest_->j_index.resize(niDest*njDest);
+
+  domainDest_->data_ni.setValue(niDest);
+  domainDest_->data_nj.setValue(njDest);
+  domainDest_->data_ibegin.setValue(0);
+  domainDest_->data_jbegin.setValue(0);
+  domainDest_->data_i_index.resize(niDest*njDest);
+  domainDest_->data_j_index.resize(niDest*njDest);
+
+  domainDest_->domainMask.resize(niDest*njDest);
+
+  if (!domainSrc_->lonvalue_1d.isEmpty())
+  {
+    if (domainDest_->type == CDomain::type_attr::rectilinear)
+    {
+      domainDest_->lonvalue_1d.resize(niDest);
+      domainDest_->latvalue_1d.resize(njDest);
+    }
+    else if (domainDest_->type == CDomain::type_attr::unstructured)
+    {
+      domainDest_->lonvalue_1d.resize(niDest);
+      domainDest_->latvalue_1d.resize(niDest);
+    }
+    else if (domainDest_->type == CDomain::type_attr::curvilinear)
+    {
+      domainDest_->lonvalue_1d.resize(niDest*njDest);
+      domainDest_->latvalue_1d.resize(niDest*njDest);
+    }
+  }
+  else if (!domainSrc_->lonvalue_2d.isEmpty())
+  {
+    domainDest_->lonvalue_2d.resize(niDest,njDest);
+    domainDest_->latvalue_2d.resize(niDest,njDest);
+  }
+
+  if (domainSrc_->hasBounds)
+  {
+    nvertex = domainSrc_->nvertex;
+    domainDest_->nvertex.setValue(nvertex);
+    if (!domainSrc_->bounds_lon_1d.isEmpty())
+    {
+      if (domainDest_->type == CDomain::type_attr::rectilinear)
+      {
+        domainDest_->bounds_lon_1d.resize(nvertex, niDest);
+        domainDest_->bounds_lat_1d.resize(nvertex, njDest);
+      }
+      else if (domainDest_->type == CDomain::type_attr::unstructured)
+      {
+        domainDest_->bounds_lon_1d.resize(nvertex, niDest);
+        domainDest_->bounds_lat_1d.resize(nvertex, niDest);
+      }
+      else if (domainDest_->type == CDomain::type_attr::curvilinear)
+      {
+        domainDest_->bounds_lon_1d.resize(nvertex, niDest*njDest);
+        domainDest_->bounds_lat_1d.resize(nvertex, niDest*njDest);
+      }
+    }
+    else if (!domainSrc_->bounds_lon_2d.isEmpty())
+    {
+      domainDest_->bounds_lon_2d.resize(nvertex, niDest, njDest);
+      domainDest_->bounds_lat_2d.resize(nvertex, niDest, njDest);
+    }
+  }
+  if (domainSrc_->hasArea) domainDest_->area.resize(niDest,njDest);
+
+  this->transformationMapping_.resize(1);
+  this->transformationWeight_.resize(1);
+  TransformationIndexMap& transMap = this->transformationMapping_[0];
+  TransformationWeightMap& transWeight = this->transformationWeight_[0];
+
+  for (int iDest = 0; iDest < niDest; iDest++)
+  {
+    iSrc = iDest + destIBegin;
+    for (int jDest = 0; jDest < njDest; jDest++)
+    {
+      jSrc = jDest + destJBegin;
+      ind = jSrc * domainSrc_->ni + iSrc;
+      iIdxSrc = domainSrc_->i_index(ind);
+      jIdxSrc = domainSrc_->j_index(ind);
+      indLocDest = jDest*niDest + iDest;
+      indGloDest = (jDest + jbeginDest)*zoomNi_ + (iDest + ibeginDest);
+      indLocSrc = (jDest+destJBegin)*domainSrc_->ni + (iDest+destIBegin);
+      indGloSrc = (jIdxSrc )* niGloSrc + iIdxSrc;
+      domainDest_->i_index(indLocDest) = iDest + ibeginDest;                                             // i_index contains global positions
+      domainDest_->j_index(indLocDest) = jDest + jbeginDest;                                             // i_index contains global positions
+      domainDest_->data_i_index(indLocDest) = (domainSrc_->data_dim == 1) ? indLocDest : iDest;          // data_i_index contains local positions
+      domainDest_->data_j_index(indLocDest) = (domainSrc_->data_dim == 1) ? 0 :jDest;                    // data_i_index contains local positions
+      domainDest_->domainMask(indLocDest) = domainSrc_->domainMask(indLocSrc);
+
+      if (domainSrc_->hasArea)
+        domainDest_->area(iDest,jDest) = domainSrc_->area(iSrc,jSrc);
+
+      if (domainSrc_->hasLonLat)
+      {
+        if (!domainSrc_->latvalue_1d.isEmpty())
+        {
+          if (domainDest_->type == CDomain::type_attr::rectilinear)
+          {
+            domainDest_->latvalue_1d(jDest) = domainSrc_->latvalue_1d(jSrc);
+          }
+          else
+          {
+            domainDest_->lonvalue_1d(indLocDest) = domainSrc_->lonvalue_1d(ind);
+            domainDest_->latvalue_1d(indLocDest) = domainSrc_->latvalue_1d(ind);
+          }
+        }
+        else if (!domainSrc_->latvalue_2d.isEmpty())
+        {
+          domainDest_->lonvalue_2d(iDest,jDest) = domainSrc_->lonvalue_2d(iSrc,jSrc);
+          domainDest_->latvalue_2d(iDest,jDest) = domainSrc_->latvalue_2d(iSrc,jSrc);
+        }
+      }
+
+      if (domainSrc_->hasBounds)
+      {
+        if (!domainSrc_->bounds_lon_1d.isEmpty())
+        {
+          if (domainDest_->type == CDomain::type_attr::rectilinear)
+          {
+            for (int n = 0; n < nvertex; ++n)
+              domainDest_->bounds_lat_1d(n,jDest) = domainSrc_->bounds_lat_1d(n,jSrc);
+          }
+          else
+          {
+            for (int n = 0; n < nvertex; ++n)
+            {
+              domainDest_->bounds_lon_1d(n,indLocDest) = domainSrc_->bounds_lon_1d(n,ind);
+              domainDest_->bounds_lat_1d(n,indLocDest) = domainSrc_->bounds_lat_1d(n,ind);
+            }
+          }
+        }
+        else if (!domainSrc_->bounds_lon_2d.isEmpty())
+        {
+          for (int n = 0; n < nvertex; ++n)
+          {
+            domainDest_->bounds_lon_2d(n,iDest,jDest) = domainSrc_->bounds_lon_2d(n,iSrc,jSrc);
+            domainDest_->bounds_lat_2d(n,iDest,jDest) = domainSrc_->bounds_lat_2d(n,iSrc,jSrc);
+          }
+        }
+
+      }
+
+      transMap[indGloDest].push_back(indGloSrc);
+      transWeight[indGloDest].push_back(1.0);
+    }
+
+    if (domainSrc_->hasLonLat && !domainSrc_->latvalue_1d.isEmpty())
+    {
+      if (domainDest_->type == CDomain::type_attr::rectilinear)
+      {
+        domainDest_->lonvalue_1d(iDest) = domainSrc_->lonvalue_1d(iSrc);
+      }
+    }
+
+    if (domainSrc_->hasBounds && !domainSrc_->bounds_lon_1d.isEmpty())
+    {
+      if (domainDest_->type == CDomain::type_attr::rectilinear)
+      {
+        for (int n = 0; n < nvertex; ++n)
+          domainDest_->bounds_lon_1d(n,iDest) = domainSrc_->bounds_lon_1d(n,iSrc);
+      }
+    }
+  }
+  domainDest_->computeLocalMask();
+
+  computeAlgorithm(domainSource->getLocalView(CElementView::WORKFLOW), domainDestination->getLocalView(CElementView::WORKFLOW)) ;
+  
+}
+CATCH
+
+/*!
+  Compute the index mapping between domain on grid source and one on grid destination
+*/
+void CDomainAlgorithmZoom::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+{
+}
+
+
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_zoom.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_zoom.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_zoom.hpp	(revision 1985)
@@ -0,0 +1,60 @@
+
+#ifndef __XIOS_DOMAIN_ALGORITHM_ZOOM_HPP__
+#define __XIOS_DOMAIN_ALGORITHM_ZOOM_HPP__
+
+#include "domain_algorithm_transformation.hpp"
+#include "transformation.hpp"
+
+namespace xios {
+
+class CDomain;
+class CZoomDomain;
+
+/*!
+  \class CDomainAlgorithmZoom
+  Implementing zoom (alternative zoom) on domain
+*/
+class CDomainAlgorithmZoom : public CDomainAlgorithmTransformation
+{
+public:
+  CDomainAlgorithmZoom(CDomain* domainDestination, CDomain* domainSource, CZoomDomain* zoomDomain);
+
+  virtual ~CDomainAlgorithmZoom() {}
+
+  static bool registerTrans();
+protected:
+  void updateDomainAttributes();
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs);
+
+private:
+  void updateZoom();
+
+private:
+  //! Global zoom begin on domain
+  int zoomIBegin_;
+  int zoomJBegin_;
+
+  //! Global zoom end on domain
+  int zoomIEnd_;
+  int zoomJEnd_;
+
+  //! Global zoom size on domain
+  int zoomNi_;
+  int zoomNj_;
+
+private:
+
+  static CGenericAlgorithmTransformation* create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                CTransformation<CDomain>* transformation,
+                                                int elementPositionInGrid,
+                                                std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridDst2DomainPosition);
+  static bool dummyRegistered_;
+};
+
+}
+#endif // __XIOS_DOMAIN_ALGORITHM_ZOOM_HPP__
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/grid_transformation.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/grid_transformation.cpp	(revision 1984)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/grid_transformation.cpp	(revision 1985)
@@ -47,5 +47,5 @@
   for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation
   CGenericAlgorithmTransformation* algo = 0;
-  algo = CGridTransformationFactory<CScalar>::createTransformation(transType,
+  algo = CGridTransformationFactory<CScalar>::createTransformation(transType, false,
                                                                   gridDestination_,
                                                                   gridSource_,
@@ -78,5 +78,5 @@
 
   CGenericAlgorithmTransformation* algo = 0;
-  algo = CGridTransformationFactory<CAxis>::createTransformation(transType,
+  algo = CGridTransformationFactory<CAxis>::createTransformation(transType, false,
                                                                  gridDestination_,
                                                                  gridSource_,
@@ -109,5 +109,5 @@
 
   CGenericAlgorithmTransformation* algo = 0;
-  algo = CGridTransformationFactory<CDomain>::createTransformation(transType,
+  algo = CGridTransformationFactory<CDomain>::createTransformation(transType, false,
                                                                    gridDestination_,
                                                                    gridSource_,
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/grid_transformation_factory_impl.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/grid_transformation_factory_impl.hpp	(revision 1984)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/grid_transformation_factory_impl.hpp	(revision 1985)
@@ -33,5 +33,5 @@
   virtual ~CGridTransformationFactory() {}
 
-  static CGenericAlgorithmTransformation* createTransformation(ETranformationType transType,
+  static CGenericAlgorithmTransformation* createTransformation(ETranformationType transType, bool isSource,
                                                                CGrid* gridDst, CGrid* gridSrc,
                                                                CTransformation<T>* transformation,
@@ -45,5 +45,5 @@
 
 public:
-  typedef CGenericAlgorithmTransformation* (*CreateTransformationCallBack)(CGrid* gridDst, CGrid* gridSrc,
+  typedef CGenericAlgorithmTransformation* (*CreateTransformationCallBack)(bool isSource, CGrid* gridDst, CGrid* gridSrc,
                                                                            CTransformation<T>* transformation,
                                                                            int elementPositionInGrid,
@@ -68,5 +68,5 @@
 
 template<typename T>
-CGenericAlgorithmTransformation* CGridTransformationFactory<T>::createTransformation(ETranformationType transType,
+CGenericAlgorithmTransformation* CGridTransformationFactory<T>::createTransformation(ETranformationType transType, bool isSource,
                                                                                CGrid* gridDst, CGrid* gridSrc,
                                                                                CTransformation<T>* transformation,
@@ -86,5 +86,5 @@
            << "doesn't exist. Please define.");
   }
-  return (it->second)(gridDst, gridSrc, transformation, elementPositionInGrid,
+  return (it->second)(isSource, gridDst, gridSrc, transformation, elementPositionInGrid,
                       elementPositionInGridSrc2ScalarPosition,
                       elementPositionInGridSrc2AxisPosition,
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_extract_axis.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_extract_axis.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_extract_axis.cpp	(revision 1985)
@@ -0,0 +1,93 @@
+/*!
+   \file scalar_algorithm_extract_scalar.cpp
+   \author Ha NGUYEN
+   \since 23 June 2016
+   \date 23 June 2016
+
+   \brief Algorithm for extract an axis to a scalar
+ */
+#include "scalar_algorithm_extract_axis.hpp"
+#include "axis.hpp"
+#include "scalar.hpp"
+#include "extract_axis_to_scalar.hpp"
+#include "grid.hpp"
+#include "grid_transformation_factory_impl.hpp"
+
+#include "reduction.hpp"
+
+namespace xios {
+CGenericAlgorithmTransformation* CScalarAlgorithmExtractAxis::create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                                     CTransformation<CScalar>* transformation,
+                                                                     int elementPositionInGrid,
+                                                                     std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                                     std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                                     std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2DomainPosition)
+TRY
+{
+  std::vector<CScalar*> scalarListDestP = gridDst->getScalars();
+  std::vector<CAxis*> axisListSrcP  = gridSrc->getAxis();
+
+  CExtractAxisToScalar* extractAxis = dynamic_cast<CExtractAxisToScalar*> (transformation);
+  int scalarDstIndex = elementPositionInGridDst2ScalarPosition[elementPositionInGrid];
+  int axisSrcIndex = elementPositionInGridSrc2AxisPosition[elementPositionInGrid];
+
+  return (new CScalarAlgorithmExtractAxis(scalarListDestP[scalarDstIndex], axisListSrcP[axisSrcIndex], extractAxis));
+}
+CATCH
+
+bool CScalarAlgorithmExtractAxis::dummyRegistered_ = CScalarAlgorithmExtractAxis::registerTrans();
+bool CScalarAlgorithmExtractAxis::registerTrans()
+TRY
+{
+  return CGridTransformationFactory<CScalar>::registerTransformation(TRANS_EXTRACT_AXIS_TO_SCALAR, create);
+}
+CATCH
+
+CScalarAlgorithmExtractAxis::CScalarAlgorithmExtractAxis(CScalar* scalarDestination, CAxis* axisSource, CExtractAxisToScalar* algo)
+ : CScalarAlgorithmTransformation(scalarDestination, axisSource),
+   reduction_(0)
+TRY
+{
+  algo->checkValid(scalarDestination, axisSource);
+  StdString op = "extract";
+  pos_ = algo->position;
+  reduction_ = CReductionAlgorithm::createOperation(CReductionAlgorithm::ReductionOperations[op]);
+}
+CATCH
+
+void CScalarAlgorithmExtractAxis::apply(const std::vector<std::pair<int,double> >& localIndex,
+                                         const double* dataInput,
+                                         CArray<double,1>& dataOut,
+                                         std::vector<bool>& flagInitial,                     
+                                         bool ignoreMissingValue, bool firstPass)
+TRY
+{
+  reduction_->apply(localIndex, dataInput, dataOut, flagInitial, ignoreMissingValue, firstPass);
+}
+CATCH
+
+CScalarAlgorithmExtractAxis::~CScalarAlgorithmExtractAxis()
+TRY
+{
+  if (0 != reduction_) delete reduction_;
+}
+CATCH
+
+void CScalarAlgorithmExtractAxis::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+TRY
+{
+  this->transformationMapping_.resize(1);
+  this->transformationWeight_.resize(1);
+
+  TransformationIndexMap& transMap = this->transformationMapping_[0];
+  TransformationWeightMap& transWeight = this->transformationWeight_[0];
+
+  transMap[0].push_back(pos_);
+  transWeight[0].push_back(1.0);
+}
+CATCH
+
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_extract_axis.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_extract_axis.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_extract_axis.hpp	(revision 1985)
@@ -0,0 +1,64 @@
+/*!
+   \file scalar_algorithm_extract_scalar.hpp
+   \author Ha NGUYEN
+   \since 23 June 2016
+   \date 23 June 2016
+
+   \brief Algorithm for extract an axis to a scalar
+ */
+#ifndef __XIOS_SCALAR_ALGORITHM_EXTRACT_AXIS_HPP__
+#define __XIOS_SCALAR_ALGORITHM_EXTRACT_AXIS_HPP__
+
+#include "scalar_algorithm_transformation.hpp"
+#include "transformation.hpp"
+
+namespace xios {
+
+class CScalar;
+class CAxis;
+class CExtractAxisToScalar;
+class CReductionAlgorithm;
+
+/*!
+  \class CScalarAlgorithmExtractAxis
+  Extract a scalar from an axis
+*/
+class CScalarAlgorithmExtractAxis : public CScalarAlgorithmTransformation
+{
+public:
+  CScalarAlgorithmExtractAxis(CScalar* scalarDestination, CAxis* axisSource, CExtractAxisToScalar* algo);
+
+  virtual void apply(const std::vector<std::pair<int,double> >& localIndex,
+                     const double* dataInput,
+                     CArray<double,1>& dataOut,
+                     std::vector<bool>& flagInitial,                     
+                     bool ignoreMissingValue, bool firstPass);
+
+  virtual ~CScalarAlgorithmExtractAxis();
+
+  static bool registerTrans();
+
+protected:
+  int pos_;
+
+protected:
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs);
+
+protected:
+  CReductionAlgorithm* reduction_;
+
+private:
+  static CGenericAlgorithmTransformation* create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                CTransformation<CScalar>* transformation,
+                                                int elementPositionInGrid,
+                                                std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridDst2DomainPosition);
+  static bool dummyRegistered_;
+};
+
+}
+#endif // __XIOS_SCALAR_ALGORITHM_EXTRACT_AXIS_HPP__
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_reduce_axis.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_reduce_axis.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_reduce_axis.cpp	(revision 1985)
@@ -0,0 +1,134 @@
+/*!
+   \file scalar_algorithm_reduce_scalar.cpp
+   \author Ha NGUYEN
+   \since 23 June 2016
+   \date 23 June 2016
+
+   \brief Algorithm for reduce an axis to a scalar
+ */
+#include "scalar_algorithm_reduce_axis.hpp"
+#include "axis.hpp"
+#include "scalar.hpp"
+#include "reduce_axis_to_scalar.hpp"
+#include "grid.hpp"
+#include "grid_transformation_factory_impl.hpp"
+#include "reduction.hpp"
+
+#include "reduction.hpp"
+
+namespace xios {
+CGenericAlgorithmTransformation* CScalarAlgorithmReduceAxis::create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                                     CTransformation<CScalar>* transformation,
+                                                                     int elementPositionInGrid,
+                                                                     std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                                     std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                                     std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2DomainPosition)
+TRY
+{
+  std::vector<CScalar*> scalarListDestP = gridDst->getScalars();
+  std::vector<CAxis*> axisListSrcP  = gridSrc->getAxis();
+
+  CReduceAxisToScalar* reduceAxis = dynamic_cast<CReduceAxisToScalar*> (transformation);
+  int scalarDstIndex = elementPositionInGridDst2ScalarPosition[elementPositionInGrid];
+  int axisSrcIndex = elementPositionInGridSrc2AxisPosition[elementPositionInGrid];
+
+  return (new CScalarAlgorithmReduceAxis(scalarListDestP[scalarDstIndex], axisListSrcP[axisSrcIndex], reduceAxis));
+}
+CATCH
+
+bool CScalarAlgorithmReduceAxis::dummyRegistered_ = CScalarAlgorithmReduceAxis::registerTrans();
+bool CScalarAlgorithmReduceAxis::registerTrans()
+TRY
+{
+  return CGridTransformationFactory<CScalar>::registerTransformation(TRANS_REDUCE_AXIS_TO_SCALAR, create);
+}
+CATCH
+
+CScalarAlgorithmReduceAxis::CScalarAlgorithmReduceAxis(CScalar* scalarDestination, CAxis* axisSource, CReduceAxisToScalar* algo)
+ : CScalarAlgorithmTransformation(scalarDestination, axisSource),
+   reduction_(0)
+TRY
+{
+  if (algo->operation.isEmpty())
+    ERROR("CScalarAlgorithmReduceAxis::CScalarAlgorithmReduceAxis(CAxis* axisDestination, CAxis* axisSource, CReduceAxisToScalar* algo)",
+           << "Operation must be defined."
+           << "Axis source " <<axisSource->getId() << std::endl
+           << "Scalar destination " << scalarDestination->getId());
+  StdString op;
+  switch (algo->operation)
+  {
+    case CReduceAxisToScalar::operation_attr::sum:
+      op = "sum";
+      break;
+    case CReduceAxisToScalar::operation_attr::min:
+      op = "min";
+      break;
+    case CReduceAxisToScalar::operation_attr::max:
+      op = "max";
+      break;
+    case CReduceAxisToScalar::operation_attr::average:
+      op = "average";
+      break;
+    default:
+        ERROR("CScalarAlgorithmReduceAxis::CScalarAlgorithmReduceAxis(CScalar* scalarDestination, CAxis* axisSource, CReduceAxisToScalar* algo)",
+         << "Operation is wrongly defined. Supported operations: sum, min, max, average." << std::endl
+         << "Domain source " <<axisSource->getId() << std::endl
+         << "Scalar destination " << scalarDestination->getId());
+
+  }
+  
+  if (CReductionAlgorithm::ReductionOperations.end() == CReductionAlgorithm::ReductionOperations.find(op))
+    ERROR("CScalarAlgorithmReduceAxis::CScalarAlgorithmReduceAxis(CAxis* axisDestination, CAxis* axisSource, CReduceAxisToScalar* algo)",
+       << "Operation '" << op << "' not found. Please make sure to use a supported one"
+       << "Axis source " <<axisSource->getId() << std::endl
+       << "Scalar destination " << scalarDestination->getId());
+
+  reduction_ = CReductionAlgorithm::createOperation(CReductionAlgorithm::ReductionOperations[op]);
+}
+CATCH
+
+void CScalarAlgorithmReduceAxis::apply(const std::vector<std::pair<int,double> >& localIndex, const double* dataInput, CArray<double,1>& dataOut,
+                                         std::vector<bool>& flagInitial, bool ignoreMissingValue, bool firstPass)
+TRY
+{
+  reduction_->apply(localIndex, dataInput, dataOut, flagInitial, ignoreMissingValue, firstPass);
+}
+CATCH
+
+void CScalarAlgorithmReduceAxis::updateData(CArray<double,1>& dataOut)
+TRY
+{
+  reduction_->updateData(dataOut);
+}
+CATCH
+
+CScalarAlgorithmReduceAxis::~CScalarAlgorithmReduceAxis()
+TRY
+{
+  if (0 != reduction_) delete reduction_;
+}
+CATCH
+
+void CScalarAlgorithmReduceAxis::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+TRY
+{
+  this->transformationMapping_.resize(1);
+  this->transformationWeight_.resize(1);
+
+  TransformationIndexMap& transMap = this->transformationMapping_[0];
+  TransformationWeightMap& transWeight = this->transformationWeight_[0];
+
+  int globalIndexSize = axisSrc_-> n_glo;
+
+  for (int idx = 0; idx < globalIndexSize; ++idx)
+  {
+    transMap[0].push_back(idx);
+    transWeight[0].push_back(1.0);
+  }
+}
+CATCH
+
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_reduce_axis.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_reduce_axis.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_reduce_axis.hpp	(revision 1985)
@@ -0,0 +1,63 @@
+/*!
+   \file scalar_algorithm_reduce_scalar.hpp
+   \author Ha NGUYEN
+   \since 23 June 2016
+   \date 23 June 2016
+
+   \brief Algorithm for reduce an axis to a scalar
+ */
+#ifndef __XIOS_SCALAR_ALGORITHM_REDUCE_AXIS_HPP__
+#define __XIOS_SCALAR_ALGORITHM_REDUCE_AXIS_HPP__
+
+#include "scalar_algorithm_transformation.hpp"
+#include "transformation.hpp"
+
+namespace xios {
+
+class CScalar;
+class CAxis;
+class CReduceAxisToScalar;
+class CReductionAlgorithm;
+
+/*!
+  \class CScalarAlgorithmReduceAxis
+  Reducing an axis to a scalar
+*/
+class CScalarAlgorithmReduceAxis : public CScalarAlgorithmTransformation
+{
+public:
+  CScalarAlgorithmReduceAxis(CScalar* scalarDestination, CAxis* axisSource, CReduceAxisToScalar* algo);
+
+  virtual void apply(const std::vector<std::pair<int,double> >& localIndex,
+                     const double* dataInput,
+                     CArray<double,1>& dataOut,
+                     std::vector<bool>& flagInitial,
+                     bool ignoreMissingValue, bool firstPass);
+
+  virtual void updateData(CArray<double,1>& dataOut);
+  
+  virtual ~CScalarAlgorithmReduceAxis();
+
+  static bool registerTrans();
+protected:
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs);
+
+protected:
+  CReductionAlgorithm* reduction_;
+
+private:
+
+  static CGenericAlgorithmTransformation* create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                CTransformation<CScalar>* transformation,
+                                                int elementPositionInGrid,
+                                                std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridDst2DomainPosition);
+  static bool dummyRegistered_;
+};
+
+}
+#endif // __XIOS_SCALAR_ALGORITHM_REDUCE_AXIS_HPP__
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_reduce_domain.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_reduce_domain.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_reduce_domain.cpp	(revision 1985)
@@ -0,0 +1,155 @@
+/*!
+   \file scalar_algorithm_reduce_domain.cpp
+   \author Ha NGUYEN
+   \since 23 June 2016
+   \date 23 June 2016
+
+   \brief Algorithm for reduce a domain to a scalar
+ */
+#include "scalar_algorithm_reduce_domain.hpp"
+#include "domain.hpp"
+#include "scalar.hpp"
+#include "reduce_domain_to_scalar.hpp"
+#include "grid.hpp"
+#include "grid_transformation_factory_impl.hpp"
+
+#include "reduction.hpp"
+
+namespace xios {
+CGenericAlgorithmTransformation* CScalarAlgorithmReduceDomain::create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                                     CTransformation<CScalar>* transformation,
+                                                                     int elementPositionInGrid,
+                                                                     std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                                     std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                                     std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2DomainPosition)
+TRY
+{
+  std::vector<CScalar*> scalarListDestP = gridDst->getScalars();
+  std::vector<CDomain*> domainListSrcP  = gridSrc->getDomains();
+
+  CReduceDomainToScalar* reduceDomain = dynamic_cast<CReduceDomainToScalar*> (transformation);
+  int scalarDstIndex = elementPositionInGridDst2ScalarPosition[elementPositionInGrid];
+  int domainSrcIndex = elementPositionInGridSrc2DomainPosition[elementPositionInGrid];
+
+  return (new CScalarAlgorithmReduceDomain(scalarListDestP[scalarDstIndex], domainListSrcP[domainSrcIndex], reduceDomain));
+}
+CATCH
+
+bool CScalarAlgorithmReduceDomain::dummyRegistered_ = CScalarAlgorithmReduceDomain::registerTrans();
+bool CScalarAlgorithmReduceDomain::registerTrans()
+TRY
+{
+  return CGridTransformationFactory<CScalar>::registerTransformation(TRANS_REDUCE_DOMAIN_TO_SCALAR, create);
+}
+CATCH
+
+CScalarAlgorithmReduceDomain::CScalarAlgorithmReduceDomain(CScalar* scalarDestination, CDomain* domainSource, CReduceDomainToScalar* algo)
+ : CScalarAlgorithmTransformation(scalarDestination, domainSource),
+   reduction_(0)
+TRY
+{
+  algo->checkValid(scalarDestination, domainSource);
+  
+  StdString op;
+  switch (algo->operation)
+  {
+    case CReduceDomainToScalar::operation_attr::sum:
+      op = "sum";
+      break;
+    case CReduceDomainToScalar::operation_attr::min:
+      op = "min";
+      break;
+    case CReduceDomainToScalar::operation_attr::max:
+      op = "max";
+      break;
+    case CReduceDomainToScalar::operation_attr::average:
+      op = "average";
+      break;
+    default:
+        ERROR("CScalarAlgorithmReduceDomain::CScalarAlgorithmReduceDomain(CDomain* domainDestination, CDomain* domainSource, CReduceDomainToScalar* algo)",
+         << "Operation must be defined."
+         << "Domain source " <<domainSource->getId() << std::endl
+         << "Scalar destination " << scalarDestination->getId());
+
+  }
+  
+  if (CReductionAlgorithm::ReductionOperations.end() == CReductionAlgorithm::ReductionOperations.find(op))
+    ERROR("CScalarAlgorithmReduceDomain::CScalarAlgorithmReduceDomain(CDomain* domainDestination, CDomain* domainSource, CReduceDomainToScalar* algo)",
+       << "Operation '" << op << "' not found. Please make sure to use a supported one"
+       << "Domain source " <<domainSource->getId() << std::endl
+       << "Scalar destination " << scalarDestination->getId());
+
+  reduction_ = CReductionAlgorithm::createOperation(CReductionAlgorithm::ReductionOperations[op]);
+  local = algo->local ;
+}
+CATCH
+
+void CScalarAlgorithmReduceDomain::apply(const std::vector<std::pair<int,double> >& localIndex,
+                                         const double* dataInput,
+                                         CArray<double,1>& dataOut,
+                                         std::vector<bool>& flagInitial,                     
+                                         bool ignoreMissingValue, bool firstPass)
+TRY
+{
+  reduction_->apply(localIndex, dataInput, dataOut, flagInitial, ignoreMissingValue, firstPass);
+}
+CATCH
+
+void CScalarAlgorithmReduceDomain::updateData(CArray<double,1>& dataOut)
+TRY
+{
+  reduction_->updateData(dataOut);
+}
+CATCH
+
+CScalarAlgorithmReduceDomain::~CScalarAlgorithmReduceDomain()
+TRY
+{
+  if (0 != reduction_) delete reduction_;
+}
+CATCH
+
+void CScalarAlgorithmReduceDomain::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+TRY
+{
+  this->transformationMapping_.resize(1);
+  this->transformationWeight_.resize(1);
+
+  TransformationIndexMap& transMap = this->transformationMapping_[0];
+  TransformationWeightMap& transWeight = this->transformationWeight_[0];
+
+  int ni_glo = domainSrc_->ni_glo ;
+  int nj_glo = domainSrc_->nj_glo ;
+  int nbDomainIdx ;
+  
+  if (local)
+  {
+      const CArray<int, 1>& i_index = domainSrc_-> i_index.getValue() ;
+      const CArray<int, 1>& j_index = domainSrc_-> j_index.getValue() ;
+      const CArray<bool,1>& localMask = domainSrc_-> localMask ;
+      int nbDomainIdx = i_index.numElements();
+
+      for (int idxDomain = 0; idxDomain < nbDomainIdx; ++idxDomain)
+      {
+        if (localMask(idxDomain))
+        { 
+          transMap[0].push_back(j_index(idxDomain)* ni_glo + i_index(idxDomain));
+          transWeight[0].push_back(1.0) ;
+        }
+      }
+  }
+  else
+  {  
+    nbDomainIdx = ni_glo * nj_glo;
+    transMap[0].resize(nbDomainIdx);
+    transWeight[0].resize(nbDomainIdx, 1.0);
+    for (int idxDomain = 0; idxDomain < nbDomainIdx; ++idxDomain) transMap[0][idxDomain] = idxDomain;    
+  }
+  
+}
+CATCH
+
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_reduce_domain.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_reduce_domain.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_reduce_domain.hpp	(revision 1985)
@@ -0,0 +1,64 @@
+/*!
+   \file scalar_algorithm_reduce_scalar.hpp
+   \author Ha NGUYEN
+   \since 13 Oct 2016
+   \date 13 Oct 2016
+
+   \brief Algorithm for reduce an DOMAIN to a scalar
+ */
+#ifndef __XIOS_SCALAR_ALGORITHM_REDUCE_DOMAIN_HPP__
+#define __XIOS_SCALAR_ALGORITHM_REDUCE_DOMAIN_HPP__
+
+#include "scalar_algorithm_transformation.hpp"
+#include "transformation.hpp"
+
+namespace xios {
+
+class CScalar;
+class CDomain;
+class CReduceDomainToScalar;
+class CReductionAlgorithm;
+
+/*!
+  \class CScalarAlgorithmReduceDomain
+  Reducing an DOMAIN to a scalar
+*/
+class CScalarAlgorithmReduceDomain : public CScalarAlgorithmTransformation
+{
+public:
+  CScalarAlgorithmReduceDomain(CScalar* scalarDestination, CDomain* domainSource, CReduceDomainToScalar* algo);
+
+  virtual void apply(const std::vector<std::pair<int,double> >& localIndex,
+                     const double* dataInput,
+                     CArray<double,1>& dataOut,
+                     std::vector<bool>& flagInitial,                     
+                     bool ignoreMissingValue, bool firstPass);
+
+  virtual void updateData(CArray<double,1>& dataOut);
+  
+  virtual ~CScalarAlgorithmReduceDomain();
+
+  static bool registerTrans();
+protected:
+  bool local ;
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs);
+
+protected:
+  CReductionAlgorithm* reduction_;
+
+private:
+
+  static CGenericAlgorithmTransformation* create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                CTransformation<CScalar>* transformation,
+                                                int elementPositionInGrid,
+                                                std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2DOMAINPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridDst2DOMAINPosition,
+                                                std::map<int, int>& elementPositionInGridDst2DomainPosition);
+  static bool dummyRegistered_;
+};
+
+}
+#endif // __XIOS_SCALAR_ALGORITHM_REDUCE_DOMAIN_HPP__
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_reduce_scalar.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_reduce_scalar.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_reduce_scalar.cpp	(revision 1985)
@@ -0,0 +1,126 @@
+/*!
+   \file scalar_algorithm_reduce_scalar.cpp
+ 
+   \brief Algorithm for reduce an scalar to a scalar
+ */
+#include "scalar_algorithm_reduce_scalar.hpp"
+#include "scalar.hpp"
+#include "reduce_scalar_to_scalar.hpp"
+#include "grid.hpp"
+#include "grid_transformation_factory_impl.hpp"
+#include "reduction.hpp"
+
+
+namespace xios {
+CGenericAlgorithmTransformation* CScalarAlgorithmReduceScalar::create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                                     CTransformation<CScalar>* transformation,
+                                                                     int elementPositionInGrid,
+                                                                     std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                                     std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                                     std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                                     std::map<int, int>& elementPositionInGridDst2DomainPosition)
+TRY
+{
+  std::vector<CScalar*> scalarListDestP = gridDst->getScalars();
+  std::vector<CScalar*> scalarListSrcP  = gridSrc->getScalars();
+
+  CReduceScalarToScalar* reduceScalar = dynamic_cast<CReduceScalarToScalar*> (transformation);
+  int scalarDstIndex = elementPositionInGridDst2ScalarPosition[elementPositionInGrid];
+  int scalarSrcIndex = elementPositionInGridSrc2AxisPosition[elementPositionInGrid];
+
+  return (new CScalarAlgorithmReduceScalar(scalarListDestP[scalarDstIndex], scalarListSrcP[scalarSrcIndex], reduceScalar));
+}
+CATCH
+
+bool CScalarAlgorithmReduceScalar::dummyRegistered_ = CScalarAlgorithmReduceScalar::registerTrans();
+bool CScalarAlgorithmReduceScalar::registerTrans()
+TRY
+{
+  return CGridTransformationFactory<CScalar>::registerTransformation(TRANS_REDUCE_SCALAR_TO_SCALAR, create);
+}
+CATCH
+
+CScalarAlgorithmReduceScalar::CScalarAlgorithmReduceScalar(CScalar* scalarDestination, CScalar* scalarSource, CReduceScalarToScalar* algo)
+ : CScalarAlgorithmTransformation(scalarDestination, scalarSource),
+   reduction_(0)
+TRY
+{
+  eliminateRedondantSrc_= false ;
+  if (algo->operation.isEmpty())
+    ERROR("CScalarAlgorithmReduceScalar::CScalarAlgorithmReduceScalar(CScalar* scalarDestination, CScalar* scalarSource, CReduceScalarToScalar* algo)",
+           << "Operation must be defined."
+           << "Scalar source " <<scalarSource->getId() << std::endl
+           << "Scalar destination " << scalarDestination->getId());
+  StdString op;
+  switch (algo->operation)
+  {
+    case CReduceScalarToScalar::operation_attr::sum:
+      op = "sum";
+      break;
+    case CReduceScalarToScalar::operation_attr::min:
+      op = "min";
+      break;
+    case CReduceScalarToScalar::operation_attr::max:
+      op = "max";
+      break;
+    case CReduceScalarToScalar::operation_attr::average:
+      op = "average";
+      break;
+    default:
+        ERROR("CScalarAlgorithmReduceScalar::CScalarAlgorithmReduceScalar(CScalar* scalarDestination, CScalar* scalarSource, CReduceScalarToScalar* algo)",
+         << "Operation is wrongly defined. Supported operations: sum, min, max, average." << std::endl
+         << "Scalar source " <<scalarSource->getId() << std::endl
+         << "Scalar destination " << scalarDestination->getId());
+
+  }
+  
+  if (CReductionAlgorithm::ReductionOperations.end() == CReductionAlgorithm::ReductionOperations.find(op))
+    ERROR("CScalarAlgorithmReduceScalar::CScalarAlgorithmReduceScalar(CScalar* scalarDestination, CScalar* scalarSource, CReduceScalarToScalar* algo)",
+       << "Operation '" << op << "' not found. Please make sure to use a supported one"
+       << "Scalar source " <<scalarSource->getId() << std::endl
+       << "Scalar destination " << scalarDestination->getId());
+
+  reduction_ = CReductionAlgorithm::createOperation(CReductionAlgorithm::ReductionOperations[op]);
+}
+CATCH
+
+void CScalarAlgorithmReduceScalar::apply(const std::vector<std::pair<int,double> >& localIndex, const double* dataInput, CArray<double,1>& dataOut,
+                                         std::vector<bool>& flagInitial, bool ignoreMissingValue, bool firstPass)
+TRY
+{
+  reduction_->apply(localIndex, dataInput, dataOut, flagInitial, ignoreMissingValue, firstPass);
+}
+CATCH
+
+void CScalarAlgorithmReduceScalar::updateData(CArray<double,1>& dataOut)
+TRY
+{
+  reduction_->updateData(dataOut);
+}
+CATCH
+
+CScalarAlgorithmReduceScalar::~CScalarAlgorithmReduceScalar()
+TRY
+{
+  if (0 != reduction_) delete reduction_;
+}
+CATCH
+
+void CScalarAlgorithmReduceScalar::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+TRY
+{
+  this->transformationMapping_.resize(1);
+  this->transformationWeight_.resize(1);
+
+  TransformationIndexMap& transMap = this->transformationMapping_[0];
+  TransformationWeightMap& transWeight = this->transformationWeight_[0];
+
+  transMap[0].push_back(0);
+  transWeight[0].push_back(1.0);
+
+}
+CATCH
+
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_reduce_scalar.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_reduce_scalar.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_reduce_scalar.hpp	(revision 1985)
@@ -0,0 +1,58 @@
+/*!
+   \file scalar_algorithm_reduce_scalar.hpp
+   \brief Algorithm for reduce an scalar to a scalar
+ */
+#ifndef __XIOS_SCALAR_ALGORITHM_REDUCE_SCALAR_HPP__
+#define __XIOS_SCALAR_ALGORITHM_REDUCE_SCALAR_HPP__
+
+#include "scalar_algorithm_transformation.hpp"
+#include "transformation.hpp"
+
+namespace xios {
+
+class CScalar;
+class CReduceScalarToScalar;
+class CReductionAlgorithm;
+
+/*!
+  \class CScalarAlgorithmReduceScalar
+  Reducing an scalar to a scalar
+*/
+class CScalarAlgorithmReduceScalar : public CScalarAlgorithmTransformation
+{
+public:
+  CScalarAlgorithmReduceScalar(CScalar* scalarDestination, CScalar* scalarSource, CReduceScalarToScalar* algo);
+
+  virtual void apply(const std::vector<std::pair<int,double> >& localIndex,
+                     const double* dataInput,
+                     CArray<double,1>& dataOut,
+                     std::vector<bool>& flagInitial,
+                     bool ignoreMissingValue, bool firstPass);
+
+  virtual void updateData(CArray<double,1>& dataOut);
+  
+  virtual ~CScalarAlgorithmReduceScalar();
+
+  static bool registerTrans();
+protected:
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs);
+
+protected:
+  CReductionAlgorithm* reduction_;
+
+private:
+
+  static CGenericAlgorithmTransformation* create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
+                                                CTransformation<CScalar>* transformation,
+                                                int elementPositionInGrid,
+                                                std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridSrc2DomainPosition,
+                                                std::map<int, int>& elementPositionInGridDst2ScalarPosition,
+                                                std::map<int, int>& elementPositionInGridDst2AxisPosition,
+                                                std::map<int, int>& elementPositionInGridDst2DomainPosition);
+  static bool dummyRegistered_;
+};
+
+}
+#endif // __XIOS_SCALAR_ALGORITHM_REDUCE_SCALAR_HPP__
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_transformation.cpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_transformation.cpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_transformation.cpp	(revision 1985)
@@ -0,0 +1,124 @@
+/*!
+   \file scalar_algorithm_transformation.hpp
+   \author Ha NGUYEN
+   \since 23 June 2016
+   \date 23 Oct 2016
+
+   \brief Interface for all scalar transformation algorithms.
+ */
+
+#include "scalar_algorithm_transformation.hpp"
+#include "context.hpp"
+#include "context_client.hpp"
+#include "client_client_dht_template.hpp"
+#include "domain.hpp"
+#include "axis.hpp"
+#include "scalar.hpp"
+
+namespace xios {
+
+CScalarAlgorithmTransformation::CScalarAlgorithmTransformation(CScalar* scalarDestination, CScalar* scalarSource)
+ : CGenericAlgorithmTransformation(),
+   scalarDest_(scalarDestination),
+   scalarSrc_(scalarSource), axisSrc_(0), domainSrc_(0)
+{
+}
+
+CScalarAlgorithmTransformation::CScalarAlgorithmTransformation(CScalar* scalarDestination, CAxis* axisSource)
+ : CGenericAlgorithmTransformation(),
+   scalarDest_(scalarDestination),
+   scalarSrc_(0), axisSrc_(axisSource), domainSrc_(0)
+{
+}
+
+CScalarAlgorithmTransformation::CScalarAlgorithmTransformation(CScalar* scalarDestination, CDomain* domainSource)
+ : CGenericAlgorithmTransformation(),
+   scalarDest_(scalarDestination),
+   scalarSrc_(0), axisSrc_(0), domainSrc_(domainSource)
+{
+}
+
+CScalarAlgorithmTransformation::~CScalarAlgorithmTransformation()
+{
+}
+
+void CScalarAlgorithmTransformation::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
+{
+}
+
+/*!
+  Compute global index of element source on different processes on knowing the global index source of transformation
+
+  \param [in] globalIndexElementSource global index source of transformation (represented in the transformation mapping)
+  \param [in] elementSourceType type of element source
+  \param [out] globalIndexElementSourceOnProc processes which contain the corresponding global index of scalar source
+*/
+void CScalarAlgorithmTransformation::computeExchangeGlobalIndex(const CArray<size_t,1>& globalIndexElementSource,
+                                                                int elementSourceType,
+                                                                CClientClientDHTInt::Index2VectorInfoTypeMap& globalIndexElementSourceOnProc)
+TRY
+{
+  CContext* context = CContext::getCurrent();
+  int clientRank = context->intraCommRank_;
+  int clientSize = context->intraCommSize_;
+
+  if (2 == elementSourceType) // Source is a domain
+  {
+    size_t globalIndex;
+    int nIndexSize = domainSrc_->i_index.numElements();
+    CArray<int,1>& iIndex = domainSrc_->i_index;
+    CArray<int,1>& jIndex = domainSrc_->j_index;
+    int niGlo = domainSrc_->ni_glo;
+    CArray<bool,1>& localMask = domainSrc_ -> localMask ;
+    CClientClientDHTInt::Index2VectorInfoTypeMap globalIndex2ProcRank;
+    globalIndex2ProcRank.rehash(std::ceil(nIndexSize/globalIndex2ProcRank.max_load_factor()));
+    for (int idx = 0; idx < nIndexSize; ++idx)
+    {
+      if (localMask(idx))
+      {
+        globalIndex = jIndex(idx) * niGlo + iIndex(idx);
+        globalIndex2ProcRank[globalIndex].resize(1);
+        globalIndex2ProcRank[globalIndex][0] = clientRank;
+      }
+    }
+
+    CClientClientDHTInt dhtIndexProcRank(globalIndex2ProcRank, context->intraComm_);
+    dhtIndexProcRank.computeIndexInfoMapping(globalIndexElementSource);
+    globalIndexElementSourceOnProc = dhtIndexProcRank.getInfoIndexMap();
+  }
+  else if (1 == elementSourceType) // Source is an axis
+  {
+    size_t globalIndex;
+    int nIndexSize = axisSrc_->index.numElements();
+    CClientClientDHTInt::Index2VectorInfoTypeMap globalIndex2ProcRank;
+    globalIndex2ProcRank.rehash(std::ceil(nIndexSize/globalIndex2ProcRank.max_load_factor()));
+    CArray<bool,1>& localMask = axisSrc_ -> mask ;
+   
+    for (int idx = 0; idx < nIndexSize; ++idx)
+    {
+      if (localMask(idx))
+      {
+        globalIndex = axisSrc_->index(idx);
+        globalIndex2ProcRank[globalIndex].resize(1);
+        globalIndex2ProcRank[globalIndex][0] = clientRank;
+      }
+    }
+
+    CClientClientDHTInt dhtIndexProcRank(globalIndex2ProcRank, context->intraComm_);
+    dhtIndexProcRank.computeIndexInfoMapping(globalIndexElementSource);
+    globalIndexElementSourceOnProc = dhtIndexProcRank.getInfoIndexMap();
+  }
+  else // scalar
+  {
+    size_t globalIndex = 0;    
+    CClientClientDHTInt::Index2VectorInfoTypeMap globalIndex2ProcRank;
+    globalIndex2ProcRank[globalIndex].resize(1);
+    globalIndex2ProcRank[globalIndex][0] = clientRank;
+
+    CClientClientDHTInt dhtIndexProcRank(globalIndex2ProcRank, context->intraComm_);
+    dhtIndexProcRank.computeIndexInfoMapping(globalIndexElementSource);
+    globalIndexElementSourceOnProc = dhtIndexProcRank.getInfoIndexMap();
+  }
+}
+CATCH
+}
Index: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_transformation.hpp
===================================================================
--- XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_transformation.hpp	(revision 1985)
+++ XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/scalar_algorithm/scalar_algorithm_transformation.hpp	(revision 1985)
@@ -0,0 +1,59 @@
+/*!
+   \file scalar_algorithm_transformation.hpp
+   \author Ha NGUYEN
+   \since 23 June 2016
+   \date 23 June 2016
+
+   \brief Interface for all scalar transformation algorithms.
+ */
+#ifndef __XIOS_SCALAR_ALGORITHM_TRANSFORMATION_HPP__
+#define __XIOS_SCALAR_ALGORITHM_TRANSFORMATION_HPP__
+
+#include "generic_algorithm_transformation.hpp"
+
+namespace xios {
+
+class CScalar;
+class CAxis;
+class CDomain;
+
+/*!
+  \class CScalarAlgorithmTransformation
+  Algorithms for scalar.
+*/
+class CScalarAlgorithmTransformation : public virtual CGenericAlgorithmTransformation
+{
+public:
+  CScalarAlgorithmTransformation(CScalar* scalarDestination, CScalar* scalarSource);
+  CScalarAlgorithmTransformation(CScalar* scalarDestination, CAxis* axisSource);
+  CScalarAlgorithmTransformation(CScalar* scalarDestination, CDomain* domainSource);
+
+  virtual ~CScalarAlgorithmTransformation();
+
+protected:
+  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs);
+
+  virtual void computeExchangeGlobalIndex(const CArray<size_t,1>& globalScalarIndex,
+                                          int elementType,
+                                          CClientClientDHTInt::Index2VectorInfoTypeMap& globalScalarIndexOnProc);
+
+protected:
+  //! Global index of an scalar on grid destination
+  std::vector<int> scalarDestGlobalIndex_;
+
+  //! Size of
+  int scalarDestGlobalSize_;
+
+    //! Scalar on grid destination
+  CScalar* scalarDest_;
+
+  //! Scalar on grid source
+  CScalar* scalarSrc_;
+  //! Axis source
+  CAxis* axisSrc_;
+  //! Domain source
+  CDomain* domainSrc_;
+};
+
+}
+#endif // __XIOS_SCALAR_ALGORITHM_TRANSFORMATION_HPP__
