Index: XIOS/dev/branch_openmp/extern/remap/src/cputime.cpp
===================================================================
--- XIOS/dev/branch_openmp/extern/remap/src/cputime.cpp	(revision 1545)
+++ XIOS/dev/branch_openmp/extern/remap/src/cputime.cpp	(revision 1642)
@@ -1,4 +1,3 @@
 #include "mpi.hpp"
-using namespace ep_lib;
 
 namespace sphereRemap {
Index: XIOS/dev/branch_openmp/extern/remap/src/elt.hpp
===================================================================
--- XIOS/dev/branch_openmp/extern/remap/src/elt.hpp	(revision 1545)
+++ XIOS/dev/branch_openmp/extern/remap/src/elt.hpp	(revision 1642)
@@ -48,4 +48,5 @@
 	int n; /* number of vertices */
 	double area;
+  double given_area ;
 	Coord x; /* barycentre */
 };
@@ -80,4 +81,5 @@
 		n    = rhs.n;
 		area = rhs.area;
+		given_area = rhs.given_area;
 		x    = rhs.x;
 		val  = rhs.val;
Index: XIOS/dev/branch_openmp/extern/remap/src/gridRemap.cpp
===================================================================
--- XIOS/dev/branch_openmp/extern/remap/src/gridRemap.cpp	(revision 1545)
+++ XIOS/dev/branch_openmp/extern/remap/src/gridRemap.cpp	(revision 1642)
@@ -11,8 +11,5 @@
 
 CRemapGrid srcGrid;
-#pragma omp threadprivate(srcGrid)
-
 CRemapGrid tgtGrid;
-#pragma omp threadprivate(tgtGrid)
 
 }
Index: XIOS/dev/branch_openmp/extern/remap/src/gridRemap.hpp
===================================================================
--- XIOS/dev/branch_openmp/extern/remap/src/gridRemap.hpp	(revision 1545)
+++ XIOS/dev/branch_openmp/extern/remap/src/gridRemap.hpp	(revision 1642)
@@ -14,6 +14,6 @@
 Coord readPole(std::istream&);
 
-//extern CRemapGrid srcGrid;
-//extern CRemapGrid tgtGrid;
+extern CRemapGrid srcGrid;
+extern CRemapGrid tgtGrid;
 
 }
Index: XIOS/dev/branch_openmp/extern/remap/src/intersect.cpp
===================================================================
--- XIOS/dev/branch_openmp/extern/remap/src/intersect.cpp	(revision 1545)
+++ XIOS/dev/branch_openmp/extern/remap/src/intersect.cpp	(revision 1642)
@@ -14,10 +14,4 @@
 
 namespace sphereRemap {
-
-extern CRemapGrid srcGrid;
-#pragma omp threadprivate(srcGrid)
-
-extern CRemapGrid tgtGrid;
-#pragma omp threadprivate(tgtGrid)
 
 using namespace std;
Index: XIOS/dev/branch_openmp/extern/remap/src/intersection_ym.cpp
===================================================================
--- XIOS/dev/branch_openmp/extern/remap/src/intersection_ym.cpp	(revision 1545)
+++ XIOS/dev/branch_openmp/extern/remap/src/intersection_ym.cpp	(revision 1642)
@@ -8,4 +8,9 @@
 #include <stdlib.h>
 #include <limits>
+#include <array>
+#include <cstdint> 
+#include "earcut.hpp"
+#include <fstream> 
+
 
 #define epsilon 1e-3  // epsilon distance ratio over side lenght for approximate small circle by great circle
@@ -14,16 +19,15 @@
 namespace sphereRemap {
 
-extern CRemapGrid srcGrid;
-#pragma omp threadprivate(srcGrid)
-
-extern CRemapGrid tgtGrid;
-#pragma omp threadprivate(tgtGrid)
-
-
 using namespace std;
 using namespace ClipperLib ;
+        
 
 double intersect_ym(Elt *a, Elt *b)
 {
+
+  using N = uint32_t;
+  using Point = array<double, 2>;
+  vector<Point> vect_points;
+  vector< vector<Point> > polyline;
 
 // transform small circle into piece of great circle if necessary
@@ -52,4 +56,5 @@
   double cos_alpha;
 
+  /// vector<p2t::Point*> polyline;
   for(int n=0; n<na;n++)
   {
@@ -58,5 +63,28 @@
     a_gno[n].y=scalarprod(dstPolygon[n],Oy)/cos_alpha ;
     a_gno[n].z=scalarprod(dstPolygon[n],Oz)/cos_alpha ; // must be equal to 1
-  }
+
+    vect_points.push_back( array<double, 2>() );
+    vect_points[n][0] = a_gno[n].x;
+    vect_points[n][1] = a_gno[n].y;
+
+  }
+
+  polyline.push_back(vect_points);
+  vector<N> indices_a_gno = mapbox::earcut<N>(polyline);
+  
+  double area_a_gno=0 ;
+  for(int i=0;i<indices_a_gno.size()/3;++i)
+    {
+      Coord x0 = Ox * polyline[0][indices_a_gno[3*i]][0] + Oy* polyline[0][indices_a_gno[3*i]][1] + Oz ;
+      Coord x1 = Ox * polyline[0][indices_a_gno[3*i+1]][0] + Oy* polyline[0][indices_a_gno[3*i+1]][1] + Oz ;
+      Coord x2 = Ox * polyline[0][indices_a_gno[3*i+2]][0] + Oy* polyline[0][indices_a_gno[3*i+2]][1] + Oz ;
+      area_a_gno+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ;
+    }
+
+  vect_points.clear();
+  polyline.clear();
+  indices_a_gno.clear();
+
+ 
 
   for(int n=0; n<nb;n++)
@@ -66,6 +94,26 @@
     b_gno[n].y=scalarprod(srcPolygon[n],Oy)/cos_alpha ;
     b_gno[n].z=scalarprod(srcPolygon[n],Oz)/cos_alpha ; // must be equal to 1
-  }
-
+
+    vect_points.push_back( array<double, 2>() );
+    vect_points[n][0] = b_gno[n].x;
+    vect_points[n][1] = b_gno[n].y;
+  }
+
+
+  polyline.push_back(vect_points);
+  vector<N> indices_b_gno = mapbox::earcut<N>(polyline);
+
+  double area_b_gno=0 ;
+  for(int i=0;i<indices_b_gno.size()/3;++i)
+    {
+      Coord x0 = Ox * polyline[0][indices_b_gno[3*i]][0] + Oy* polyline[0][indices_b_gno[3*i]][1] + Oz ;
+      Coord x1 = Ox * polyline[0][indices_b_gno[3*i+1]][0] + Oy* polyline[0][indices_b_gno[3*i+1]][1] + Oz ;
+      Coord x2 = Ox * polyline[0][indices_b_gno[3*i+2]][0] + Oy* polyline[0][indices_b_gno[3*i+2]][1] + Oz ;
+      area_b_gno+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ;
+    }
+
+  vect_points.clear();
+  polyline.clear();
+  indices_b_gno.clear();
 
 
@@ -116,38 +164,70 @@
   clip.AddPaths(dst, ptClip, true);
   clip.Execute(ctIntersection, intersection);
-
+ 
   double area=0 ;
 
   for(int ni=0;ni<intersection.size(); ni++)
   {
-    // go back into real coordinate on the sphere
     Coord* intersectPolygon=new Coord[intersection[ni].size()] ;
     for(int n=0; n < intersection[ni].size(); n++)
     {
-      double x=intersection[ni][n].X/xscale+xoffset ;
-      double y=intersection[ni][n].Y/yscale+yoffset ;
-
-      intersectPolygon[n]=Ox*x+Oy*y+Oz ;
-      intersectPolygon[n]=intersectPolygon[n]*(1./norm(intersectPolygon[n])) ;
-    }
-
-// remove redondants vertex
-    int nv=0 ;
+      intersectPolygon[n].x=intersection[ni][n].X/xscale+xoffset ;
+      intersectPolygon[n].y=intersection[ni][n].Y/yscale+yoffset ;
+    }
+    
+
+    int nv=0;
+
     for(int n=0; n < intersection[ni].size(); n++)
     {
-      if (norm(intersectPolygon[n]-intersectPolygon[(n+1)%intersection[ni].size()])>fusion_vertex)
+       double dx=intersectPolygon[n].x-intersectPolygon[(n+1)%intersection[ni].size()].x ;
+       double dy=intersectPolygon[n].y-intersectPolygon[(n+1)%intersection[ni].size()].y ;
+     
+       if (dx*dx+dy*dy>fusion_vertex*fusion_vertex)
+       {
+          intersectPolygon[nv]=intersectPolygon[n] ;
+
+	  vect_points.push_back( array<double, 2>() );
+          vect_points[nv][0] = intersectPolygon[n].x;
+	  vect_points[nv][1] = intersectPolygon[n].y;
+
+	  nv++ ;
+       }
+      
+
+    }
+
+    polyline.push_back(vect_points);
+    vect_points.clear();
+
+    if (nv>2)
+    {
+ 
+      vector<N> indices = mapbox::earcut<N>(polyline);
+
+      double area2=0 ;
+      for(int i=0;i<indices.size()/3;++i)
       {
-        intersectPolygon[nv]=intersectPolygon[n] ;
-        nv++ ;
+	  Coord x0 = Ox * polyline[0][indices[3*i]][0] + Oy* polyline[0][indices[3*i]][1] + Oz ;
+	  Coord x1 = Ox * polyline[0][indices[3*i+1]][0] + Oy* polyline[0][indices[3*i+1]][1] + Oz ;
+	  Coord x2 = Ox * polyline[0][indices[3*i+2]][0] + Oy* polyline[0][indices[3*i+2]][1] + Oz ;
+	  area2+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ;
       }
-    }
-
-
-    if (nv>2)
-    {
+
+      polyline.clear();
+
+      for(int n=0; n < nv; n++)
+      {
+        intersectPolygon[n] = Ox*intersectPolygon[n].x+Oy*intersectPolygon[n].y+Oz;
+        intersectPolygon[n] = intersectPolygon[n]*(1./norm(intersectPolygon[n])) ;
+      }
+
+
 //     assign intersection to source and destination polygons
        Polyg *is = new Polyg;
        is->x = exact_barycentre(intersectPolygon,nv);
-       is->area = polygonarea(intersectPolygon,nv) ;
+//       is->area = polygonarea(intersectPolygon,nv) ;
+       is->area = area2 ;
+
 //        if (is->area < 1e-12) cout<<"Small intersection : "<<is->area<<endl ;
        if (is->area==0.) delete is ;
@@ -168,5 +248,8 @@
   delete[] b_gno ;
   return area ;
+
 }
+
+
 
 void createGreatCirclePolygon(const Elt& element, const Coord& pole, vector<Coord>& coordinates)
Index: XIOS/dev/branch_openmp/extern/remap/src/libmapper.cpp
===================================================================
--- XIOS/dev/branch_openmp/extern/remap/src/libmapper.cpp	(revision 1545)
+++ XIOS/dev/branch_openmp/extern/remap/src/libmapper.cpp	(revision 1642)
@@ -15,6 +15,4 @@
 #include "cputime.hpp" // cputime
 
-using namespace ep_lib;
-
 using namespace sphereRemap ;
 
@@ -22,5 +20,5 @@
    and deallocated during the second step (computing the weights) */
 Mapper *mapper;
-#pragma omp threadprivate(mapper)
+
 
 /** xxx_bounds_yyy is of length n_vert_per_cell_xxx*n_cell_xxx
@@ -43,9 +41,10 @@
 	assert(n_cell_dst >= 4);
 	assert(1 <= order && order <= 2);
-
-  mapper = new Mapper(MPI_COMM_WORLD);
+  double* src_area=NULL ;
+  double* dst_area=NULL ;
+  mapper = new Mapper(EP_COMM_WORLD);
   mapper->setVerbosity(PROGRESS) ;
-  mapper->setSourceMesh(src_bounds_lon, src_bounds_lat, n_vert_per_cell_src, n_cell_src, src_pole ) ;
-  mapper->setTargetMesh(dst_bounds_lon, dst_bounds_lat, n_vert_per_cell_dst, n_cell_dst, dst_pole ) ;
+  mapper->setSourceMesh(src_bounds_lon, src_bounds_lat, src_area, n_vert_per_cell_src, n_cell_src, src_pole ) ;
+  mapper->setTargetMesh(dst_bounds_lon, dst_bounds_lat, dst_area, n_vert_per_cell_dst, n_cell_dst, dst_pole ) ;
 
 /*
Index: XIOS/dev/branch_openmp/extern/remap/src/mapper.cpp
===================================================================
--- XIOS/dev/branch_openmp/extern/remap/src/mapper.cpp	(revision 1545)
+++ XIOS/dev/branch_openmp/extern/remap/src/mapper.cpp	(revision 1642)
@@ -12,14 +12,6 @@
 
 #include "mapper.hpp"
-using namespace ep_lib;
 
 namespace sphereRemap {
-
-extern CRemapGrid srcGrid;
-#pragma omp threadprivate(srcGrid)
-
-extern CRemapGrid tgtGrid;
-#pragma omp threadprivate(tgtGrid)
-
 
 /* A subdivition of an array into N sub-arays
@@ -35,11 +27,11 @@
 
 
-void Mapper::setSourceMesh(const double* boundsLon, const double* boundsLat, int nVertex, int nbCells, const double* pole, const long int* globalId)
+void Mapper::setSourceMesh(const double* boundsLon, const double* boundsLat, const double* area, int nVertex, int nbCells, const double* pole, const long int* globalId)
 {
   srcGrid.pole = Coord(pole[0], pole[1], pole[2]);
 
   int mpiRank, mpiSize;
-  MPI_Comm_rank(communicator, &mpiRank);
-  MPI_Comm_size(communicator, &mpiSize);
+  ep_lib::MPI_Comm_rank(communicator, &mpiRank);
+  ep_lib::MPI_Comm_size(communicator, &mpiSize);
 
   sourceElements.reserve(nbCells);
@@ -51,5 +43,5 @@
     long int offset ;
     long int nb=nbCells ;
-    MPI_Scan(&nb,&offset,1,MPI_LONG,MPI_SUM,communicator) ;
+    ep_lib::MPI_Scan(&nb,&offset,1,EP_LONG,EP_SUM,communicator) ;
     offset=offset-nb ;
     for(int i=0;i<nbCells;i++) sourceGlobalId[i]=offset+i ;
@@ -67,15 +59,17 @@
     sourceMesh.push_back(Node(elt.x, cptRadius(elt), &sourceElements.back()));
     cptEltGeom(sourceElements[i], Coord(pole[0], pole[1], pole[2]));
-  }
-
-}
-
-void Mapper::setTargetMesh(const double* boundsLon, const double* boundsLat, int nVertex, int nbCells, const double* pole, const long int* globalId)
+    if (area!=NULL) sourceElements[i].given_area=area[i] ;
+    else sourceElements[i].given_area=sourceElements[i].area ;
+  }
+
+}
+
+void Mapper::setTargetMesh(const double* boundsLon, const double* boundsLat, const double* area, int nVertex, int nbCells, const double* pole, const long int* globalId)
 {
   tgtGrid.pole = Coord(pole[0], pole[1], pole[2]);
 
   int mpiRank, mpiSize;
-  MPI_Comm_rank(communicator, &mpiRank);
-  MPI_Comm_size(communicator, &mpiSize);
+  ep_lib::MPI_Comm_rank(communicator, &mpiRank);
+  ep_lib::MPI_Comm_size(communicator, &mpiSize);
 
   targetElements.reserve(nbCells);
@@ -87,5 +81,5 @@
     long int offset ;
     long int nb=nbCells ;
-    MPI_Scan(&nb,&offset,1,MPI_LONG,MPI_SUM,communicator) ;
+    ep_lib::MPI_Scan(&nb,&offset,1,EP_LONG,EP_SUM,communicator) ;
     offset=offset-nb ;
     for(int i=0;i<nbCells;i++) targetGlobalId[i]=offset+i ;
@@ -100,4 +94,6 @@
     targetMesh.push_back(Node(elt.x, cptRadius(elt), &sourceElements.back()));
     cptEltGeom(targetElements[i], Coord(pole[0], pole[1], pole[2]));
+    if (area!=NULL) targetElements[i].given_area=area[i] ;
+    else targetElements[i].given_area=targetElements[i].area ;
   }
 
@@ -121,6 +117,6 @@
   vector<double> timings;
   int mpiSize, mpiRank;
-  MPI_Comm_size(communicator, &mpiSize);
-  MPI_Comm_rank(communicator, &mpiRank);
+  ep_lib::MPI_Comm_size(communicator, &mpiSize);
+  ep_lib::MPI_Comm_rank(communicator, &mpiRank);
 
   this->buildSSTree(sourceMesh, targetMesh);
@@ -132,6 +128,5 @@
 
   tic = cputime();
-  if (interpOrder == 2) 
-  {
+  if (interpOrder == 2) {
     if (mpiRank == 0 && verbose) cout << "Computing grads ..." << endl;
     buildMeshTopology();
@@ -142,5 +137,5 @@
   /* Prepare computation of weights */
   /* compute number of intersections which for the first order case
-     corresponds to the number of edges in the remap matrix */
+           corresponds to the number of edges in the remap matrix */
   int nIntersections = 0;
   for (int j = 0; j < targetElements.size(); j++)
@@ -178,6 +173,6 @@
 {
   int mpiSize, mpiRank;
-  MPI_Comm_size(communicator, &mpiSize);
-  MPI_Comm_rank(communicator, &mpiRank);
+  ep_lib::MPI_Comm_size(communicator, &mpiSize);
+  ep_lib::MPI_Comm_rank(communicator, &mpiRank);
 
   /* create list of intersections (super mesh elements) for each rank */
@@ -187,7 +182,5 @@
     Elt& e = elements[j];
     for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++)
-    {
       elementList[(*it)->id.rank].insert(pair<int, Polyg *>((*it)->id.ind, *it));
-    }
   }
 
@@ -196,4 +189,5 @@
   double **recvValue = new double*[mpiSize];
   double **recvArea = new double*[mpiSize];
+  double **recvGivenArea = new double*[mpiSize];
   Coord **recvGrad = new Coord*[mpiSize];
   GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */
@@ -218,4 +212,5 @@
       recvValue[rank]   = new double[nbSendElement[rank]];
       recvArea[rank]    = new double[nbSendElement[rank]];
+      recvGivenArea[rank] = new double[nbSendElement[rank]];
       if (order == 2)
       {
@@ -240,6 +235,5 @@
   /* communicate sizes of source elements to be sent (index lists and later values and gradients) */
   int *nbRecvElement = new int[mpiSize];
-  MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator);
-  
+  ep_lib::MPI_Alltoall(nbSendElement, 1, EP_INT, nbRecvElement, 1, EP_INT, communicator);
 
   /* communicate indices of source elements on other ranks whoes value and gradient we need (since intersection) */
@@ -249,13 +243,14 @@
   double **sendValue = new double*[mpiSize];
   double **sendArea = new double*[mpiSize];
+  double **sendGivenArea = new double*[mpiSize];
   Coord **sendGrad = new Coord*[mpiSize];
   GloId **sendNeighIds = new GloId*[mpiSize];
-  MPI_Request *sendRequest = new MPI_Request[4*mpiSize];
-  MPI_Request *recvRequest = new MPI_Request[4*mpiSize];
+  ep_lib::MPI_Request *sendRequest = new ep_lib::MPI_Request[5*mpiSize];
+  ep_lib::MPI_Request *recvRequest = new ep_lib::MPI_Request[5*mpiSize];
   for (int rank = 0; rank < mpiSize; rank++)
   {
     if (nbSendElement[rank] > 0)
     {
-      MPI_Issend(sendElement[rank], nbSendElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]);
+      ep_lib::MPI_Issend(sendElement[rank], nbSendElement[rank], EP_INT, rank, 0, communicator, &sendRequest[nbSendRequest]);
       nbSendRequest++;
     }
@@ -266,4 +261,5 @@
       sendValue[rank]   = new double[nbRecvElement[rank]];
       sendArea[rank]   = new double[nbRecvElement[rank]];
+      sendGivenArea[rank] = new double[nbRecvElement[rank]];
       if (order == 2)
       {
@@ -275,12 +271,12 @@
         sendNeighIds[rank] = new GloId[nbRecvElement[rank]];
       }
-      MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+      ep_lib::MPI_Irecv(recvElement[rank], nbRecvElement[rank], EP_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]);
       nbRecvRequest++;
     }
   }
-  MPI_Status *status = new MPI_Status[4*mpiSize];
-
-  MPI_Waitall(nbSendRequest, sendRequest, status);
-  MPI_Waitall(nbRecvRequest, recvRequest, status);
+  ep_lib::MPI_Status *status = new ep_lib::MPI_Status[5*mpiSize];
+  
+  ep_lib::MPI_Waitall(nbSendRequest, sendRequest, status);
+        ep_lib::MPI_Waitall(nbRecvRequest, recvRequest, status);
 
   /* for all indices that have been received from requesting ranks: pack values and gradients, then send */
@@ -296,7 +292,9 @@
         sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val;
         sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area;
+        sendGivenArea[rank][j] = sstree.localElements[recvElement[rank][j]].given_area;
         if (order == 2)
         {
           sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad;
+//          cout<<"grad  "<<jj<<"  "<<recvElement[rank][j]<<"  "<<sendGrad[rank][jj]<<" "<<sstree.localElements[recvElement[rank][j]].grad<<endl ;
           sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id;
           jj++;
@@ -304,61 +302,65 @@
           {
             sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i];
+//            cout<<"grad  "<<jj<<"  "<<sendGrad[rank][jj]<<" "<<sstree.localElements[recvElement[rank][j]].grad<<endl ;
             sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].neighId[i];
-						jj++;
-					}
-				}
-				else
-					sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id;
-			}
-			MPI_Issend(sendValue[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);
-			nbSendRequest++;
-			MPI_Issend(sendArea[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 1, communicator, &sendRequest[nbSendRequest]);
-			nbSendRequest++;
-			if (order == 2)
-			{
-				MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1),
-								MPI_DOUBLE, rank, 2, communicator, &sendRequest[nbSendRequest]);
-				nbSendRequest++;
-				MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), MPI_INT, rank, 3, communicator, &sendRequest[nbSendRequest]);
+            jj++;
+          }
+        }
+        else
+          sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id;
+      }
+      ep_lib::MPI_Issend(sendValue[rank],  nbRecvElement[rank], EP_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);
+      nbSendRequest++;
+      ep_lib::MPI_Issend(sendArea[rank],  nbRecvElement[rank], EP_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);
+      nbSendRequest++;
+      ep_lib::MPI_Issend(sendGivenArea[rank],  nbRecvElement[rank], EP_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);
+      nbSendRequest++;
+      if (order == 2)
+      {
+        ep_lib::MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1), EP_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);
+        nbSendRequest++;
+        ep_lib::MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), EP_INT, rank, 0, communicator, &sendRequest[nbSendRequest]);
 //ym  --> attention taille GloId
-				nbSendRequest++;
-			}
-			else
-			{
-				MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], MPI_INT, rank, 4, communicator, &sendRequest[nbSendRequest]);
+        nbSendRequest++;
+      }
+      else
+      {
+        ep_lib::MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], EP_INT, rank, 0, communicator, &sendRequest[nbSendRequest]);
 //ym  --> attention taille GloId
-				nbSendRequest++;
-			}
-		}
-		if (nbSendElement[rank] > 0)
-		{
-			MPI_Irecv(recvValue[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
-			nbRecvRequest++;
-			MPI_Irecv(recvArea[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 1, communicator, &recvRequest[nbRecvRequest]);
-			nbRecvRequest++;
-			if (order == 2)
-			{
-				MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1),
-						MPI_DOUBLE, rank, 2, communicator, &recvRequest[nbRecvRequest]);
-				nbRecvRequest++;
-				MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), MPI_INT, rank, 3, communicator, &recvRequest[nbRecvRequest]);
-//ym  --> attention taille GloId
-				nbRecvRequest++;
-			}
-			else
-			{
-				MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], MPI_INT, rank, 4, communicator, &recvRequest[nbRecvRequest]);
+        nbSendRequest++;
+      }
+    }
+    if (nbSendElement[rank] > 0)
+    {
+      ep_lib::MPI_Irecv(recvValue[rank],  nbSendElement[rank], EP_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+      nbRecvRequest++;
+      ep_lib::MPI_Irecv(recvArea[rank],  nbSendElement[rank], EP_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+      nbRecvRequest++;
+      ep_lib::MPI_Irecv(recvGivenArea[rank],  nbSendElement[rank], EP_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+      nbRecvRequest++;
+      if (order == 2)
+      {
+        ep_lib::MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1),
+            EP_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+        nbRecvRequest++;
+        ep_lib::MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), EP_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]);
 //ym  --> attention taille GloId
         nbRecvRequest++;
       }
+      else
+      {
+        ep_lib::MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], EP_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+//ym  --> attention taille GloId
+        nbRecvRequest++;
+      }
     }
   }
         
-  MPI_Waitall(nbSendRequest, sendRequest, status);
-  MPI_Waitall(nbRecvRequest, recvRequest, status);
-	
+        ep_lib::MPI_Waitall(nbSendRequest, sendRequest, status);
+  ep_lib::MPI_Waitall(nbRecvRequest, recvRequest, status);
+  
 
   /* now that all values and gradients are available use them to computed interpolated values on target
-    and also to compute weights */
+     and also to compute weights */
   int i = 0;
   for (int j = 0; j < nbElements; j++)
@@ -367,6 +369,6 @@
 
     /* since for the 2nd order case source grid elements can contribute to a destination grid element over several "paths"
-    (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid)
-    accumulate them so that there is only one final weight between two elements */
+       (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid)
+       accumulate them so that there is only one final weight between two elements */
     map<GloId,double> wgt_map;
 
@@ -380,6 +382,8 @@
       double fk = recvValue[rank][n1];
       double srcArea = recvArea[rank][n1];
+      double srcGivenArea = recvGivenArea[rank][n1];
       double w = (*it)->area;
       if (quantity) w/=srcArea ;
+      else w=w*srcGivenArea/srcArea*e.area/e.given_area ;
 
       /* first order: src value times weight (weight = supermesh area), later divide by target area */
@@ -402,5 +406,8 @@
     double renorm=0;
     if (renormalize) 
-      for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second / e.area;
+    {
+      if (quantity) for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second ;
+      else for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second / e.area;
+    }
     else renorm=1. ;
 
@@ -417,6 +424,4 @@
     }
   }
-        
-        //MPI_Barrier(communicator);
 
   /* free all memory allocated in this function */
@@ -428,4 +433,5 @@
       delete[] recvValue[rank];
       delete[] recvArea[rank];
+      delete[] recvGivenArea[rank];
       if (order == 2)
       {
@@ -439,4 +445,5 @@
       delete[] sendValue[rank];
       delete[] sendArea[rank];
+      delete[] sendGivenArea[rank];
       if (order == 2)
         delete[] sendGrad[rank];
@@ -463,14 +470,14 @@
 void Mapper::computeGrads()
 {
-	/* array of pointers to collect local elements and elements received from other cpu */
-	vector<Elt*> globalElements(sstree.nbLocalElements + nbNeighbourElements);
-	int index = 0;
-	for (int i = 0; i < sstree.nbLocalElements; i++, index++)
-		globalElements[index] = &(sstree.localElements[i]);
-	for (int i = 0; i < nbNeighbourElements; i++, index++)
-		globalElements[index] = &neighbourElements[i];
-
-	update_baryc(sstree.localElements, sstree.nbLocalElements);
-	computeGradients(&globalElements[0], sstree.nbLocalElements);
+  /* array of pointers to collect local elements and elements received from other cpu */
+  vector<Elt*> globalElements(sstree.nbLocalElements + nbNeighbourElements);
+  int index = 0;
+  for (int i = 0; i < sstree.nbLocalElements; i++, index++)
+    globalElements[index] = &(sstree.localElements[i]);
+  for (int i = 0; i < nbNeighbourElements; i++, index++)
+    globalElements[index] = &neighbourElements[i];
+
+  update_baryc(sstree.localElements, sstree.nbLocalElements);
+  computeGradients(&globalElements[0], sstree.nbLocalElements);
 }
 
@@ -479,236 +486,237 @@
 void Mapper::buildMeshTopology()
 {
-	int mpiSize, mpiRank;
-	MPI_Comm_size(communicator, &mpiSize);
-	MPI_Comm_rank(communicator, &mpiRank);
-
-	vector<Node> *routingList = new vector<Node>[mpiSize];
-	vector<vector<int> > routes(sstree.localTree.leafs.size());
-
-	sstree.routeIntersections(routes, sstree.localTree.leafs);
-
-	for (int i = 0; i < routes.size(); ++i)
-		for (int k = 0; k < routes[i].size(); ++k)
-			routingList[routes[i][k]].push_back(sstree.localTree.leafs[i]);
-	routingList[mpiRank].clear();
-
-
-	CMPIRouting mpiRoute(communicator);
-	mpiRoute.init(routes);
-	int nRecv = mpiRoute.getTotalSourceElement();
-
-	int *nbSendNode = new int[mpiSize];
-	int *nbRecvNode = new int[mpiSize];
-	int *sendMessageSize = new int[mpiSize];
-	int *recvMessageSize = new int[mpiSize];
-
-	for (int rank = 0; rank < mpiSize; rank++)
-	{
-		nbSendNode[rank] = routingList[rank].size();
-		sendMessageSize[rank] = 0;
-		for (size_t j = 0; j < routingList[rank].size(); j++)
-		{
-			Elt *elt = (Elt *) (routingList[rank][j].data);
-			sendMessageSize[rank] += packedPolygonSize(*elt);
-		}
-	}
-
-	MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator);
-	MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator);
-
-	char **sendBuffer = new char*[mpiSize];
-	char **recvBuffer = new char*[mpiSize];
-	int *pos = new int[mpiSize];
-
-	for (int rank = 0; rank < mpiSize; rank++)
-	{
-		if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sendMessageSize[rank]];
-		if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]];
-	}
-
-	for (int rank = 0; rank < mpiSize; rank++)
-	{
-		pos[rank] = 0;
-		for (size_t j = 0; j < routingList[rank].size(); j++)
-		{
-			Elt *elt = (Elt *) (routingList[rank][j].data);
-			packPolygon(*elt, sendBuffer[rank], pos[rank]);
-		}
-	}
-	delete [] routingList;
-
-
-	int nbSendRequest = 0;
-	int nbRecvRequest = 0;
-	MPI_Request *sendRequest = new MPI_Request[mpiSize];
-	MPI_Request *recvRequest = new MPI_Request[mpiSize];
-	MPI_Status  *status      = new MPI_Status[mpiSize];
-
-	for (int rank = 0; rank < mpiSize; rank++)
-	{
-		if (nbSendNode[rank] > 0)
-		{
-			MPI_Issend(sendBuffer[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
-			nbSendRequest++;
-		}
-		if (nbRecvNode[rank] > 0)
-		{
-			MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
-			nbRecvRequest++;
-		}
-	}
-
-	MPI_Waitall(nbRecvRequest, recvRequest, status);
-	MPI_Waitall(nbSendRequest, sendRequest, status);
-
-	for (int rank = 0; rank < mpiSize; rank++)
-		if (nbSendNode[rank] > 0) delete [] sendBuffer[rank];
-	delete [] sendBuffer;
-
-	char **sendBuffer2 = new char*[mpiSize];
-	char **recvBuffer2 = new char*[mpiSize];
-
-	for (int rank = 0; rank < mpiSize; rank++)
-	{
-		nbSendNode[rank] = 0;
-		sendMessageSize[rank] = 0;
-
-		if (nbRecvNode[rank] > 0)
-		{
-			set<NodePtr> neighbourList;
-			pos[rank] = 0;
-			for (int j = 0; j < nbRecvNode[rank]; j++)
-			{
-				Elt elt;
-				unpackPolygon(elt, recvBuffer[rank], pos[rank]);
-				Node node(elt.x, cptRadius(elt), &elt);
-				findNeighbour(sstree.localTree.root, &node, neighbourList);
-			}
-			nbSendNode[rank] = neighbourList.size();
-			for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++)
-			{
-				Elt *elt = (Elt *) ((*it)->data);
-				sendMessageSize[rank] += packedPolygonSize(*elt);
-			}
-
-			sendBuffer2[rank] = new char[sendMessageSize[rank]];
-			pos[rank] = 0;
-
-			for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++)
-			{
-				Elt *elt = (Elt *) ((*it)->data);
-				packPolygon(*elt, sendBuffer2[rank], pos[rank]);
-			}
-		}
-	}
-	for (int rank = 0; rank < mpiSize; rank++)
-		if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank];
-	delete [] recvBuffer;
-
-
-	MPI_Barrier(communicator);
-	MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator);
-	MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator);
-
-	for (int rank = 0; rank < mpiSize; rank++)
-		if (nbRecvNode[rank] > 0) recvBuffer2[rank] = new char[recvMessageSize[rank]];
-
-	nbSendRequest = 0;
-	nbRecvRequest = 0;
-
-	for (int rank = 0; rank < mpiSize; rank++)
-	{
-		if (nbSendNode[rank] > 0)
-		{
-			MPI_Issend(sendBuffer2[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
-			nbSendRequest++;
-		}
-		if (nbRecvNode[rank] > 0)
-		{
-			MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
-			nbRecvRequest++;
-		}
-	}
-
-	MPI_Waitall(nbRecvRequest, recvRequest, status);
-	MPI_Waitall(nbSendRequest, sendRequest, status);
-
-	int nbNeighbourNodes = 0;
-	for (int rank = 0; rank < mpiSize; rank++)
-		nbNeighbourNodes += nbRecvNode[rank];
-
-	neighbourElements = new Elt[nbNeighbourNodes];
-	nbNeighbourElements = nbNeighbourNodes;
-
-	int index = 0;
-	for (int rank = 0; rank < mpiSize; rank++)
-	{
-		pos[rank] = 0;
-		for (int j = 0; j < nbRecvNode[rank]; j++)
-		{
-			unpackPolygon(neighbourElements[index], recvBuffer2[rank], pos[rank]);
-			neighbourElements[index].id.ind = sstree.localTree.leafs.size() + index;
-			index++;
-		}
-	}
-	for (int rank = 0; rank < mpiSize; rank++)
-	{
-		if (nbRecvNode[rank] > 0) delete [] recvBuffer2[rank];
-		if (nbSendNode[rank] > 0) delete [] sendBuffer2[rank];
-	}
-	delete [] recvBuffer2;
-	delete [] sendBuffer2;
-	delete [] sendMessageSize;
-	delete [] recvMessageSize;
-	delete [] nbSendNode;
-	delete [] nbRecvNode;
-	delete [] sendRequest;
-	delete [] recvRequest;
-	delete [] status;
-	delete [] pos;
-
-	/* re-compute on received elements to avoid having to send this information */
-	neighbourNodes.resize(nbNeighbourNodes);
-	setCirclesAndLinks(neighbourElements, neighbourNodes);
-	cptAllEltsGeom(neighbourElements, nbNeighbourNodes, srcGrid.pole);
-
-	/* the local SS tree must include nodes from other cpus if they are potential
+  int mpiSize, mpiRank;
+  ep_lib::MPI_Comm_size(communicator, &mpiSize);
+  ep_lib::MPI_Comm_rank(communicator, &mpiRank);
+
+  vector<Node> *routingList = new vector<Node>[mpiSize];
+  vector<vector<int> > routes(sstree.localTree.leafs.size());
+
+  sstree.routeIntersections(routes, sstree.localTree.leafs);
+
+  for (int i = 0; i < routes.size(); ++i)
+    for (int k = 0; k < routes[i].size(); ++k)
+      routingList[routes[i][k]].push_back(sstree.localTree.leafs[i]);
+  routingList[mpiRank].clear();
+
+
+  CMPIRouting mpiRoute(communicator);
+  mpiRoute.init(routes);
+  int nRecv = mpiRoute.getTotalSourceElement();
+// cout << mpiRank << " NRECV " << nRecv << "(" << routes.size() << ")"<< endl;
+
+  int *nbSendNode = new int[mpiSize];
+  int *nbRecvNode = new int[mpiSize];
+  int *sendMessageSize = new int[mpiSize];
+  int *recvMessageSize = new int[mpiSize];
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    nbSendNode[rank] = routingList[rank].size();
+    sendMessageSize[rank] = 0;
+    for (size_t j = 0; j < routingList[rank].size(); j++)
+    {
+      Elt *elt = (Elt *) (routingList[rank][j].data);
+      sendMessageSize[rank] += packedPolygonSize(*elt);
+    }
+  }
+
+  ep_lib::MPI_Alltoall(nbSendNode, 1, EP_INT, nbRecvNode, 1, EP_INT, communicator);
+  ep_lib::MPI_Alltoall(sendMessageSize, 1, EP_INT, recvMessageSize, 1, EP_INT, communicator);
+
+  char **sendBuffer = new char*[mpiSize];
+  char **recvBuffer = new char*[mpiSize];
+  int *pos = new int[mpiSize];
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sendMessageSize[rank]];
+    if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]];
+  }
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    pos[rank] = 0;
+    for (size_t j = 0; j < routingList[rank].size(); j++)
+    {
+      Elt *elt = (Elt *) (routingList[rank][j].data);
+      packPolygon(*elt, sendBuffer[rank], pos[rank]);
+    }
+  }
+  delete [] routingList;
+
+
+  int nbSendRequest = 0;
+  int nbRecvRequest = 0;
+  ep_lib::MPI_Request *sendRequest = new ep_lib::MPI_Request[mpiSize];
+  ep_lib::MPI_Request *recvRequest = new ep_lib::MPI_Request[mpiSize];
+  ep_lib::MPI_Status  *status      = new ep_lib::MPI_Status[mpiSize];
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    if (nbSendNode[rank] > 0)
+    {
+      ep_lib::MPI_Issend(sendBuffer[rank], sendMessageSize[rank], EP_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
+      nbSendRequest++;
+    }
+    if (nbRecvNode[rank] > 0)
+    {
+      ep_lib::MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], EP_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+      nbRecvRequest++;
+    }
+  }
+
+  ep_lib::MPI_Waitall(nbRecvRequest, recvRequest, status);
+  ep_lib::MPI_Waitall(nbSendRequest, sendRequest, status);
+
+  for (int rank = 0; rank < mpiSize; rank++)
+    if (nbSendNode[rank] > 0) delete [] sendBuffer[rank];
+  delete [] sendBuffer;
+
+  char **sendBuffer2 = new char*[mpiSize];
+  char **recvBuffer2 = new char*[mpiSize];
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    nbSendNode[rank] = 0;
+    sendMessageSize[rank] = 0;
+
+    if (nbRecvNode[rank] > 0)
+    {
+      set<NodePtr> neighbourList;
+      pos[rank] = 0;
+      for (int j = 0; j < nbRecvNode[rank]; j++)
+      {
+        Elt elt;
+        unpackPolygon(elt, recvBuffer[rank], pos[rank]);
+        Node node(elt.x, cptRadius(elt), &elt);
+        findNeighbour(sstree.localTree.root, &node, neighbourList);
+      }
+      nbSendNode[rank] = neighbourList.size();
+      for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++)
+      {
+        Elt *elt = (Elt *) ((*it)->data);
+        sendMessageSize[rank] += packedPolygonSize(*elt);
+      }
+
+      sendBuffer2[rank] = new char[sendMessageSize[rank]];
+      pos[rank] = 0;
+
+      for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++)
+      {
+        Elt *elt = (Elt *) ((*it)->data);
+        packPolygon(*elt, sendBuffer2[rank], pos[rank]);
+      }
+    }
+  }
+  for (int rank = 0; rank < mpiSize; rank++)
+    if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank];
+  delete [] recvBuffer;
+
+
+  ep_lib::MPI_Barrier(communicator);
+  ep_lib::MPI_Alltoall(nbSendNode, 1, EP_INT, nbRecvNode, 1, EP_INT, communicator);
+  ep_lib::MPI_Alltoall(sendMessageSize, 1, EP_INT, recvMessageSize, 1, EP_INT, communicator);
+
+  for (int rank = 0; rank < mpiSize; rank++)
+    if (nbRecvNode[rank] > 0) recvBuffer2[rank] = new char[recvMessageSize[rank]];
+
+  nbSendRequest = 0;
+  nbRecvRequest = 0;
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    if (nbSendNode[rank] > 0)
+    {
+      ep_lib::MPI_Issend(sendBuffer2[rank], sendMessageSize[rank], EP_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
+      nbSendRequest++;
+    }
+    if (nbRecvNode[rank] > 0)
+    {
+      ep_lib::MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], EP_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+      nbRecvRequest++;
+    }
+  }
+
+  ep_lib::MPI_Waitall(nbRecvRequest, recvRequest, status);
+  ep_lib::MPI_Waitall(nbSendRequest, sendRequest, status);
+
+  int nbNeighbourNodes = 0;
+  for (int rank = 0; rank < mpiSize; rank++)
+    nbNeighbourNodes += nbRecvNode[rank];
+
+  neighbourElements = new Elt[nbNeighbourNodes];
+  nbNeighbourElements = nbNeighbourNodes;
+
+  int index = 0;
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    pos[rank] = 0;
+    for (int j = 0; j < nbRecvNode[rank]; j++)
+    {
+      unpackPolygon(neighbourElements[index], recvBuffer2[rank], pos[rank]);
+      neighbourElements[index].id.ind = sstree.localTree.leafs.size() + index;
+      index++;
+    }
+  }
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    if (nbRecvNode[rank] > 0) delete [] recvBuffer2[rank];
+    if (nbSendNode[rank] > 0) delete [] sendBuffer2[rank];
+  }
+  delete [] recvBuffer2;
+  delete [] sendBuffer2;
+  delete [] sendMessageSize;
+  delete [] recvMessageSize;
+  delete [] nbSendNode;
+  delete [] nbRecvNode;
+  delete [] sendRequest;
+  delete [] recvRequest;
+  delete [] status;
+  delete [] pos;
+
+  /* re-compute on received elements to avoid having to send this information */
+  neighbourNodes.resize(nbNeighbourNodes);
+  setCirclesAndLinks(neighbourElements, neighbourNodes);
+  cptAllEltsGeom(neighbourElements, nbNeighbourNodes, srcGrid.pole);
+
+  /* the local SS tree must include nodes from other cpus if they are potential
            intersector of a local node */
-	sstree.localTree.insertNodes(neighbourNodes);
-
-	/* for every local element,
+  sstree.localTree.insertNodes(neighbourNodes);
+
+  /* for every local element,
            use the SS-tree to find all elements (including neighbourElements)
            who are potential neighbours because their circles intersect,
-	   then check all canditates for common edges to build up connectivity information
-	*/
-	for (int j = 0; j < sstree.localTree.leafs.size(); j++)
-	{
-		Node& node = sstree.localTree.leafs[j];
-
-		/* find all leafs whoes circles that intersect node's circle and save into node->intersectors */
-		node.search(sstree.localTree.root);
-
-		Elt *elt = (Elt *)(node.data);
-
-		for (int i = 0; i < elt->n; i++) elt->neighbour[i] = NOT_FOUND;
-
-		/* for element `elt` loop through all nodes in the SS-tree
+     then check all canditates for common edges to build up connectivity information
+  */
+  for (int j = 0; j < sstree.localTree.leafs.size(); j++)
+  {
+    Node& node = sstree.localTree.leafs[j];
+
+    /* find all leafs whoes circles that intersect node's circle and save into node->intersectors */
+    node.search(sstree.localTree.root);
+
+    Elt *elt = (Elt *)(node.data);
+
+    for (int i = 0; i < elt->n; i++) elt->neighbour[i] = NOT_FOUND;
+
+    /* for element `elt` loop through all nodes in the SS-tree
                    whoes circles intersect with the circle around `elt` (the SS intersectors)
                    and check if they are neighbours in the sense that the two elements share an edge.
                    If they do, save this information for elt */
-		for (list<NodePtr>::iterator it = (node.intersectors).begin(); it != (node.intersectors).end(); ++it)
-		{
-			Elt *elt2 = (Elt *)((*it)->data);
-			set_neighbour(*elt, *elt2);
-		}
+    for (list<NodePtr>::iterator it = (node.intersectors).begin(); it != (node.intersectors).end(); ++it)
+    {
+      Elt *elt2 = (Elt *)((*it)->data);
+      set_neighbour(*elt, *elt2);
+    }
 
 /*
-		for (int i = 0; i < elt->n; i++)
-		{
-			if (elt->neighbour[i] == NOT_FOUND)
-				error_exit("neighbour not found");
-		}
+    for (int i = 0; i < elt->n; i++)
+    {
+      if (elt->neighbour[i] == NOT_FOUND)
+        error_exit("neighbour not found");
+    }
 */
-	}
+  }
 }
 
@@ -716,225 +724,225 @@
 void Mapper::computeIntersection(Elt *elements, int nbElements)
 {
-	int mpiSize, mpiRank;
-	MPI_Comm_size(communicator, &mpiSize);
-	MPI_Comm_rank(communicator, &mpiRank);
-
-	MPI_Barrier(communicator);
-
-	vector<Node> *routingList = new vector<Node>[mpiSize];
-
-	vector<Node> routeNodes;  routeNodes.reserve(nbElements);
-	for (int j = 0; j < nbElements; j++)
-	{
-		elements[j].id.ind = j;
-		elements[j].id.rank = mpiRank;
-		routeNodes.push_back(Node(elements[j].x, cptRadius(elements[j]), &elements[j]));
-	}
-
-	vector<vector<int> > routes(routeNodes.size());
-	sstree.routeIntersections(routes, routeNodes);
-	for (int i = 0; i < routes.size(); ++i)
-		for (int k = 0; k < routes[i].size(); ++k)
-			routingList[routes[i][k]].push_back(routeNodes[i]);
-
-	if (verbose >= 2)
-	{
-		cout << " --> rank  " << mpiRank << " nbElements " << nbElements << " : ";
-		for (int rank = 0; rank < mpiSize; rank++)
-			cout << routingList[rank].size() << "   ";
-		cout << endl;
-	}
-	MPI_Barrier(communicator);
-
-	int *nbSendNode = new int[mpiSize];
-	int *nbRecvNode = new int[mpiSize];
-	int *sentMessageSize = new int[mpiSize];
-	int *recvMessageSize = new int[mpiSize];
-
-	for (int rank = 0; rank < mpiSize; rank++)
-	{
-		nbSendNode[rank] = routingList[rank].size();
-		sentMessageSize[rank] = 0;
-		for (size_t j = 0; j < routingList[rank].size(); j++)
-		{
-			Elt *elt = (Elt *) (routingList[rank][j].data);
-			sentMessageSize[rank] += packedPolygonSize(*elt);
-		}
-	}
-
-	MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator);
-	MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator);
-
-	int total = 0;
-
-	for (int rank = 0; rank < mpiSize; rank++)
-	{
-		total = total + nbRecvNode[rank];
-	}
-
-	if (verbose >= 2) cout << "---> rank " << mpiRank << " : compute intersection : total received nodes  " << total << endl;
-	char **sendBuffer = new char*[mpiSize];
-	char **recvBuffer = new char*[mpiSize];
-	int *pos = new int[mpiSize];
-
-	for (int rank = 0; rank < mpiSize; rank++)
-	{
-		if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sentMessageSize[rank]];
-		if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]];
-	}
-
-	for (int rank = 0; rank < mpiSize; rank++)
-	{
-		pos[rank] = 0;
-		for (size_t j = 0; j < routingList[rank].size(); j++)
-		{
-			Elt* elt = (Elt *) (routingList[rank][j].data);
-			packPolygon(*elt, sendBuffer[rank], pos[rank]);
-		}
-	}
-	delete [] routingList;
-
-	int nbSendRequest = 0;
-	int nbRecvRequest = 0;
-	MPI_Request *sendRequest = new MPI_Request[mpiSize];
-	MPI_Request *recvRequest = new MPI_Request[mpiSize];
-	MPI_Status   *status = new MPI_Status[mpiSize];
-
-	for (int rank = 0; rank < mpiSize; rank++)
-	{
-		if (nbSendNode[rank] > 0)
-		{
-			MPI_Issend(sendBuffer[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
-			nbSendRequest++;
-		}
-		if (nbRecvNode[rank] > 0)
-		{
-			MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
-			nbRecvRequest++;
-		}
-	}
-
-	MPI_Waitall(nbRecvRequest, recvRequest, status);
-	MPI_Waitall(nbSendRequest, sendRequest, status);
-	char **sendBuffer2 = new char*[mpiSize];
-	char **recvBuffer2 = new char*[mpiSize];
-
-	double tic = cputime();
-	for (int rank = 0; rank < mpiSize; rank++)
-	{
-		sentMessageSize[rank] = 0;
-
-		if (nbRecvNode[rank] > 0)
-		{
-			Elt *recvElt = new Elt[nbRecvNode[rank]];
-			pos[rank] = 0;
-			for (int j = 0; j < nbRecvNode[rank]; j++)
-			{
-				unpackPolygon(recvElt[j], recvBuffer[rank], pos[rank]);
-				cptEltGeom(recvElt[j], tgtGrid.pole);
-				Node recvNode(recvElt[j].x, cptRadius(recvElt[j]), &recvElt[j]);
-				recvNode.search(sstree.localTree.root);
-				/* for a node holding an element of the target, loop throught candidates for intersecting source */
-				for (list<NodePtr>::iterator it = (recvNode.intersectors).begin(); it != (recvNode.intersectors).end(); ++it)
-				{
-					Elt *elt2 = (Elt *) ((*it)->data);
-					/* recvElt is target, elt2 is source */
-//					intersect(&recvElt[j], elt2);
-					intersect_ym(&recvElt[j], elt2);
-				}
-
-				if (recvElt[j].is.size() > 0) sentMessageSize[rank] += packIntersectionSize(recvElt[j]);
-
-				// here recvNode goes out of scope
-			}
-
-			if (sentMessageSize[rank] > 0)
-			{
-				sentMessageSize[rank] += sizeof(int);
-				sendBuffer2[rank] = new char[sentMessageSize[rank]];
-				*((int *) sendBuffer2[rank]) = 0;
-				pos[rank] = sizeof(int);
-				for (int j = 0; j < nbRecvNode[rank]; j++)
-				{
-					packIntersection(recvElt[j], sendBuffer2[rank], pos[rank]);
-					//FIXME should be deleted: recvElt[j].delete_intersections(); // intersection areas have been packed to buffer and won't be used any more
-				}
-			}
-			delete [] recvElt;
-
-		}
-	}
-	delete [] pos;
-
-	for (int rank = 0; rank < mpiSize; rank++)
-	{
-		if (nbSendNode[rank] > 0) delete [] sendBuffer[rank];
-		if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank];
-		nbSendNode[rank] = 0;
-	}
-
-	if (verbose >= 2) cout << "Rank " << mpiRank << "  Compute (internal) intersection " << cputime() - tic << " s" << endl;
-	MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator);
-
-	for (int rank = 0; rank < mpiSize; rank++)
-		if (recvMessageSize[rank] > 0)
-			recvBuffer2[rank] = new char[recvMessageSize[rank]];
-
-	nbSendRequest = 0;
-	nbRecvRequest = 0;
-
-	for (int rank = 0; rank < mpiSize; rank++)
-	{
-		if (sentMessageSize[rank] > 0)
-		{
-			MPI_Issend(sendBuffer2[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
-			nbSendRequest++;
-		}
-		if (recvMessageSize[rank] > 0)
-		{
-			MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
-			nbRecvRequest++;
-		}
-	}
-
-	MPI_Waitall(nbRecvRequest, recvRequest, status);
-	MPI_Waitall(nbSendRequest, sendRequest, status);
-
-	delete [] sendRequest;
-	delete [] recvRequest;
-	delete [] status;
-	for (int rank = 0; rank < mpiSize; rank++)
-	{
-		if (nbRecvNode[rank] > 0)
-		{
-			if (sentMessageSize[rank] > 0)
-				delete [] sendBuffer2[rank];
-		}
-
-		if (recvMessageSize[rank] > 0)
-		{
-			unpackIntersection(elements, recvBuffer2[rank]);
-			delete [] recvBuffer2[rank];
-		}
-	}
-	delete [] sendBuffer2;
-	delete [] recvBuffer2;
-	delete [] sendBuffer;
-	delete [] recvBuffer;
-
-	delete [] nbSendNode;
-	delete [] nbRecvNode;
-	delete [] sentMessageSize;
-	delete [] recvMessageSize;
+  int mpiSize, mpiRank;
+  ep_lib::MPI_Comm_size(communicator, &mpiSize);
+  ep_lib::MPI_Comm_rank(communicator, &mpiRank);
+
+  ep_lib::MPI_Barrier(communicator);
+
+  vector<Node> *routingList = new vector<Node>[mpiSize];
+
+  vector<Node> routeNodes;  routeNodes.reserve(nbElements);
+  for (int j = 0; j < nbElements; j++)
+  {
+    elements[j].id.ind = j;
+    elements[j].id.rank = mpiRank;
+    routeNodes.push_back(Node(elements[j].x, cptRadius(elements[j]), &elements[j]));
+  }
+
+  vector<vector<int> > routes(routeNodes.size());
+  sstree.routeIntersections(routes, routeNodes);
+  for (int i = 0; i < routes.size(); ++i)
+    for (int k = 0; k < routes[i].size(); ++k)
+      routingList[routes[i][k]].push_back(routeNodes[i]);
+
+  if (verbose >= 2)
+  {
+    cout << " --> rank  " << mpiRank << " nbElements " << nbElements << " : ";
+    for (int rank = 0; rank < mpiSize; rank++)
+      cout << routingList[rank].size() << "   ";
+    cout << endl;
+  }
+  ep_lib::MPI_Barrier(communicator);
+
+  int *nbSendNode = new int[mpiSize];
+  int *nbRecvNode = new int[mpiSize];
+  int *sentMessageSize = new int[mpiSize];
+  int *recvMessageSize = new int[mpiSize];
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    nbSendNode[rank] = routingList[rank].size();
+    sentMessageSize[rank] = 0;
+    for (size_t j = 0; j < routingList[rank].size(); j++)
+    {
+      Elt *elt = (Elt *) (routingList[rank][j].data);
+      sentMessageSize[rank] += packedPolygonSize(*elt);
+    }
+  }
+
+  ep_lib::MPI_Alltoall(nbSendNode, 1, EP_INT, nbRecvNode, 1, EP_INT, communicator);
+  ep_lib::MPI_Alltoall(sentMessageSize, 1, EP_INT, recvMessageSize, 1, EP_INT, communicator);
+
+  int total = 0;
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    total = total + nbRecvNode[rank];
+  }
+
+  if (verbose >= 2) cout << "---> rank " << mpiRank << " : compute intersection : total received nodes  " << total << endl;
+  char **sendBuffer = new char*[mpiSize];
+  char **recvBuffer = new char*[mpiSize];
+  int *pos = new int[mpiSize];
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sentMessageSize[rank]];
+    if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]];
+  }
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    pos[rank] = 0;
+    for (size_t j = 0; j < routingList[rank].size(); j++)
+    {
+      Elt* elt = (Elt *) (routingList[rank][j].data);
+      packPolygon(*elt, sendBuffer[rank], pos[rank]);
+    }
+  }
+  delete [] routingList;
+
+  int nbSendRequest = 0;
+  int nbRecvRequest = 0;
+  ep_lib::MPI_Request *sendRequest = new ep_lib::MPI_Request[mpiSize];
+  ep_lib::MPI_Request *recvRequest = new ep_lib::MPI_Request[mpiSize];
+  ep_lib::MPI_Status   *status = new ep_lib::MPI_Status[mpiSize];
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    if (nbSendNode[rank] > 0)
+    {
+      ep_lib::MPI_Issend(sendBuffer[rank], sentMessageSize[rank], EP_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
+      nbSendRequest++;
+    }
+    if (nbRecvNode[rank] > 0)
+    {
+      ep_lib::MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], EP_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+      nbRecvRequest++;
+    }
+  }
+
+  ep_lib::MPI_Waitall(nbRecvRequest, recvRequest, status);
+  ep_lib::MPI_Waitall(nbSendRequest, sendRequest, status);
+  char **sendBuffer2 = new char*[mpiSize];
+  char **recvBuffer2 = new char*[mpiSize];
+
+  double tic = cputime();
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    sentMessageSize[rank] = 0;
+
+    if (nbRecvNode[rank] > 0)
+    {
+      Elt *recvElt = new Elt[nbRecvNode[rank]];
+      pos[rank] = 0;
+      for (int j = 0; j < nbRecvNode[rank]; j++)
+      {
+        unpackPolygon(recvElt[j], recvBuffer[rank], pos[rank]);
+        cptEltGeom(recvElt[j], tgtGrid.pole);
+        Node recvNode(recvElt[j].x, cptRadius(recvElt[j]), &recvElt[j]);
+        recvNode.search(sstree.localTree.root);
+        /* for a node holding an element of the target, loop throught candidates for intersecting source */
+        for (list<NodePtr>::iterator it = (recvNode.intersectors).begin(); it != (recvNode.intersectors).end(); ++it)
+        {
+          Elt *elt2 = (Elt *) ((*it)->data);
+          /* recvElt is target, elt2 is source */
+//          intersect(&recvElt[j], elt2);
+          intersect_ym(&recvElt[j], elt2);
+        }
+
+        if (recvElt[j].is.size() > 0) sentMessageSize[rank] += packIntersectionSize(recvElt[j]);
+
+        // here recvNode goes out of scope
+      }
+
+      if (sentMessageSize[rank] > 0)
+      {
+        sentMessageSize[rank] += sizeof(int);
+        sendBuffer2[rank] = new char[sentMessageSize[rank]];
+        *((int *) sendBuffer2[rank]) = 0;
+        pos[rank] = sizeof(int);
+        for (int j = 0; j < nbRecvNode[rank]; j++)
+        {
+          packIntersection(recvElt[j], sendBuffer2[rank], pos[rank]);
+          //FIXME should be deleted: recvElt[j].delete_intersections(); // intersection areas have been packed to buffer and won't be used any more
+        }
+      }
+      delete [] recvElt;
+
+    }
+  }
+  delete [] pos;
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    if (nbSendNode[rank] > 0) delete [] sendBuffer[rank];
+    if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank];
+    nbSendNode[rank] = 0;
+  }
+
+  if (verbose >= 2) cout << "Rank " << mpiRank << "  Compute (internal) intersection " << cputime() - tic << " s" << endl;
+  ep_lib::MPI_Alltoall(sentMessageSize, 1, EP_INT, recvMessageSize, 1, EP_INT, communicator);
+
+  for (int rank = 0; rank < mpiSize; rank++)
+    if (recvMessageSize[rank] > 0)
+      recvBuffer2[rank] = new char[recvMessageSize[rank]];
+
+  nbSendRequest = 0;
+  nbRecvRequest = 0;
+
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    if (sentMessageSize[rank] > 0)
+    {
+      ep_lib::MPI_Issend(sendBuffer2[rank], sentMessageSize[rank], EP_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
+      nbSendRequest++;
+    }
+    if (recvMessageSize[rank] > 0)
+    {
+      ep_lib::MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], EP_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
+      nbRecvRequest++;
+    }
+  }
+
+  ep_lib::MPI_Waitall(nbRecvRequest, recvRequest, status);
+  ep_lib::MPI_Waitall(nbSendRequest, sendRequest, status);
+
+  delete [] sendRequest;
+  delete [] recvRequest;
+  delete [] status;
+  for (int rank = 0; rank < mpiSize; rank++)
+  {
+    if (nbRecvNode[rank] > 0)
+    {
+      if (sentMessageSize[rank] > 0)
+        delete [] sendBuffer2[rank];
+    }
+
+    if (recvMessageSize[rank] > 0)
+    {
+      unpackIntersection(elements, recvBuffer2[rank]);
+      delete [] recvBuffer2[rank];
+    }
+  }
+  delete [] sendBuffer2;
+  delete [] recvBuffer2;
+  delete [] sendBuffer;
+  delete [] recvBuffer;
+
+  delete [] nbSendNode;
+  delete [] nbRecvNode;
+  delete [] sentMessageSize;
+  delete [] recvMessageSize;
 }
 
 Mapper::~Mapper()
 {
-	delete [] remapMatrix;
-	delete [] srcAddress;
-	delete [] srcRank;
-	delete [] dstAddress;
-	if (neighbourElements) delete [] neighbourElements;
-}
-
-}
+  delete [] remapMatrix;
+  delete [] srcAddress;
+  delete [] srcRank;
+  delete [] dstAddress;
+  if (neighbourElements) delete [] neighbourElements;
+}
+
+}
Index: XIOS/dev/branch_openmp/extern/remap/src/mapper.hpp
===================================================================
--- XIOS/dev/branch_openmp/extern/remap/src/mapper.hpp	(revision 1545)
+++ XIOS/dev/branch_openmp/extern/remap/src/mapper.hpp	(revision 1642)
@@ -18,11 +18,10 @@
 {
 public:
-       Mapper(ep_lib::MPI_Comm comm) : communicator(comm), verbose(SILENT), neighbourElements(NULL), sstree(comm) {}
-
+       Mapper(ep_lib::MPI_Comm comm=EP_COMM_WORLD) : communicator(comm), verbose(SILENT), neighbourElements(NULL), sstree(comm) {}
        ~Mapper();
        void setVerbosity(verbosity v) {verbose=v ;}
 
-       void setSourceMesh(const double* boundsLon, const double* boundsLat, int nVertex, int nbCells, const double* pole, const long int* globalId=NULL) ;
-       void setTargetMesh(const double* boundsLon, const double* boundsLat, int nVertex, int nbCells, const double* pole, const long int* globalId=NULL) ;
+       void setSourceMesh(const double* boundsLon, const double* boundsLat, const double* area, int nVertex, int nbCells, const double* pole, const long int* globalId=NULL) ;
+       void setTargetMesh(const double* boundsLon, const double* boundsLat, const double* area, int nVertex, int nbCells, const double* pole, const long int* globalId=NULL) ;
        void setSourceValue(const double* val) ;
        void getTargetValue(double* val) ;
Index: XIOS/dev/branch_openmp/extern/remap/src/meshutil.cpp
===================================================================
--- XIOS/dev/branch_openmp/extern/remap/src/meshutil.cpp	(revision 1545)
+++ XIOS/dev/branch_openmp/extern/remap/src/meshutil.cpp	(revision 1642)
@@ -2,8 +2,65 @@
 #include "elt.hpp"
 #include "polyg.hpp"
+#include "intersection_ym.hpp"
+#include "earcut.hpp"
+#include <vector>
 
 namespace sphereRemap {
 
 using namespace std;
+
+double computePolygoneArea(Elt& a, const Coord &pole)
+{
+  using N = uint32_t;
+  using Point = array<double, 2>;
+  vector<Point> vect_points;
+  vector< vector<Point> > polyline;
+  
+  vector<Coord> dstPolygon ;
+  createGreatCirclePolygon(a, pole, dstPolygon) ;
+
+  int na=dstPolygon.size() ;
+  Coord *a_gno   = new Coord[na];
+  
+  Coord OC=barycentre(a.vertex,a.n) ;
+  Coord Oz=OC ;
+  Coord Ox=crossprod(Coord(0,0,1),Oz) ;
+// choose Ox not too small to avoid rounding error
+  if (norm(Ox)< 0.1) Ox=crossprod(Coord(0,1,0),Oz) ;
+  Ox=Ox*(1./norm(Ox)) ;
+  Coord Oy=crossprod(Oz,Ox) ;
+  double cos_alpha;
+
+  for(int n=0; n<na;n++)
+  {
+    cos_alpha=scalarprod(OC,dstPolygon[n]) ;
+    a_gno[n].x=scalarprod(dstPolygon[n],Ox)/cos_alpha ;
+    a_gno[n].y=scalarprod(dstPolygon[n],Oy)/cos_alpha ;
+    a_gno[n].z=scalarprod(dstPolygon[n],Oz)/cos_alpha ; // must be equal to 1
+
+    vect_points.push_back( array<double, 2>() );
+    vect_points[n][0] = a_gno[n].x;
+    vect_points[n][1] = a_gno[n].y;
+
+  }
+
+  polyline.push_back(vect_points);
+  vector<N> indices_a_gno = mapbox::earcut<N>(polyline);
+  
+  double area_a_gno=0 ;
+  for(int i=0;i<indices_a_gno.size()/3;++i)
+    {
+      Coord x0 = Ox * polyline[0][indices_a_gno[3*i]][0] + Oy* polyline[0][indices_a_gno[3*i]][1] + Oz ;
+      Coord x1 = Ox * polyline[0][indices_a_gno[3*i+1]][0] + Oy* polyline[0][indices_a_gno[3*i+1]][1] + Oz ;
+      Coord x2 = Ox * polyline[0][indices_a_gno[3*i+2]][0] + Oy* polyline[0][indices_a_gno[3*i+2]][1] + Oz ;
+      area_a_gno+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ;
+    }
+
+  vect_points.clear();
+  polyline.clear();
+  indices_a_gno.clear();
+  return area_a_gno ;
+}
+
 
 void cptEltGeom(Elt& elt, const Coord &pole)
@@ -14,5 +71,9 @@
   elt.area = airbar(elt.n, elt.vertex, elt.edge, elt.d, pole, gg);
   elt.x = gg;
-}
+// overwrite area computation 
+
+  elt.area =  computePolygoneArea(elt, pole) ;
+}
+
 
 void cptAllEltsGeom(Elt *elt, int N, const Coord &pole)
@@ -133,5 +194,6 @@
   {
     for (int i = 0; i < elts[j]->n; i++) 
-      if ( elts[j]->neighbour[i]== NOT_FOUND) neighbours[i]=NULL ;
+      if ( elts[j]->neighbour[i]== NOT_FOUND) neighbours[i]=NULL ; // no neighbour
+      else if (elts[elts[j]->neighbour[i]]->is.size() == 0) neighbours[i]=NULL ; // neighbour has none supermesh cell 
       else  neighbours[i] = elts[elts[j]->neighbour[i]];
 
Index: XIOS/dev/branch_openmp/extern/remap/src/mpi_cascade.cpp
===================================================================
--- XIOS/dev/branch_openmp/extern/remap/src/mpi_cascade.cpp	(revision 1545)
+++ XIOS/dev/branch_openmp/extern/remap/src/mpi_cascade.cpp	(revision 1642)
@@ -1,12 +1,11 @@
 #include "mpi_cascade.hpp"
 #include <iostream>
-using namespace ep_lib;
 
 namespace sphereRemap {
 
-CMPICascade::CMPICascade(int nodes_per_level, MPI_Comm comm)
+CMPICascade::CMPICascade(int nodes_per_level, ep_lib::MPI_Comm comm)
 {
 	int remaining_levels;
-	MPI_Comm intraComm;
+	ep_lib::MPI_Comm intraComm;
 	int l = 0; // current level
 	do {
@@ -16,6 +15,6 @@
 		level[l].p_grp_size = level[l].size/level[l].group_size;
 	
-		MPI_Comm_split(comm, level[l].colour(), level[l].key(), &intraComm);
-		MPI_Comm_split(comm, level[l].p_colour(), level[l].p_key(), &(level[l].pg_comm));
+		ep_lib::MPI_Comm_split(comm, level[l].colour(), level[l].key(), &intraComm);
+		ep_lib::MPI_Comm_split(comm, level[l].p_colour(), level[l].p_key(), &(level[l].pg_comm));
 		comm = intraComm;
 		l++;
Index: XIOS/dev/branch_openmp/extern/remap/src/mpi_cascade.hpp
===================================================================
--- XIOS/dev/branch_openmp/extern/remap/src/mpi_cascade.hpp	(revision 1545)
+++ XIOS/dev/branch_openmp/extern/remap/src/mpi_cascade.hpp	(revision 1642)
@@ -12,21 +12,21 @@
 {
 public:
-  CCascadeLevel(ep_lib::MPI_Comm comm) : comm(comm)
-  {
-    ep_lib::MPI_Comm_size(comm, &size);
-    ep_lib::MPI_Comm_rank(comm, &rank);
-  }
-  int colour() const { return rank % group_size; };
-  int key() const { return p_colour() + rank/(p_grp_size*group_size)*p_grp_size; }
+	CCascadeLevel(ep_lib::MPI_Comm comm) : comm(comm)
+	{
+		ep_lib::MPI_Comm_size(comm, &size);
+		ep_lib::MPI_Comm_rank(comm, &rank);
+	}
+	int colour() const { return rank % group_size; };
+	int key() const { return p_colour() + rank/(p_grp_size*group_size)*p_grp_size; }
 
-  // perpendicular group
-  int p_colour() const { return (rank%group_size + rank/group_size) % p_grp_size; }
-  int p_key() const { return colour() + rank/(p_grp_size*group_size)*group_size; }
+	// perpendicular group
+	int p_colour() const { return (rank%group_size + rank/group_size) % p_grp_size; }
+	int p_key() const { return colour() + rank/(p_grp_size*group_size)*group_size; }
 
-  ep_lib::MPI_Comm comm, pg_comm;
-  int rank;
-  int size;
-  int group_size; // group_size and p_grp_size are interchanged?? FIXME
-  int p_grp_size;
+	ep_lib::MPI_Comm comm, pg_comm;
+	int rank;
+	int size;
+	int group_size; // group_size and p_grp_size are interchanged?? FIXME
+	int p_grp_size;
 };
 
@@ -34,8 +34,9 @@
 {
 public:
-  CMPICascade(int nodes_per_level, ep_lib::MPI_Comm comm);
+	// 
+	CMPICascade(int nodes_per_level, ep_lib::MPI_Comm comm);
 
-  int num_levels;
-  std::vector<CCascadeLevel> level;
+	int num_levels;
+	std::vector<CCascadeLevel> level;
 };
 
Index: XIOS/dev/branch_openmp/extern/remap/src/mpi_routing.cpp
===================================================================
--- XIOS/dev/branch_openmp/extern/remap/src/mpi_routing.cpp	(revision 1545)
+++ XIOS/dev/branch_openmp/extern/remap/src/mpi_routing.cpp	(revision 1642)
@@ -5,5 +5,4 @@
 #include "timerRemap.hpp"
 #include <iostream>
-using namespace ep_lib;
 
 namespace sphereRemap {
@@ -11,8 +10,8 @@
 const int verbose = 0;
 
-CMPIRouting::CMPIRouting(MPI_Comm comm) : communicator(comm)
-{
-	MPI_Comm_rank(comm, &mpiRank);
-	MPI_Comm_size(comm, &mpiSize);
+CMPIRouting::CMPIRouting(ep_lib::MPI_Comm comm) : communicator(comm)
+{
+	ep_lib::MPI_Comm_rank(comm, &mpiRank);
+	ep_lib::MPI_Comm_size(comm, &mpiSize);
 }
 
@@ -20,8 +19,8 @@
     but message lengths are *known* to receiver */
 template <typename T>
-void alltoalls_known(const vector<vector<T> >& send, vector<vector<T> >& recv, const vector<int>& ranks, MPI_Comm communicator)
-{
-	vector<MPI_Request> request(ranks.size() * 2);
-	vector<MPI_Status>  status(ranks.size() * 2);
+void alltoalls_known(const vector<vector<T> >& send, vector<vector<T> >& recv, const vector<int>& ranks, ep_lib::MPI_Comm communicator)
+{
+	vector<ep_lib::MPI_Request> request(ranks.size() * 2);
+	vector<ep_lib::MPI_Status>  status(ranks.size() * 2);
 
 	// communicate data
@@ -29,9 +28,9 @@
 	for (int i = 0; i < ranks.size(); i++)
 		if (recv[i].size())
-			MPI_Irecv(&recv[i][0], recv[i].size()*sizeof(T), MPI_CHAR, ranks[i], 0, communicator, &request[nbRequest++]);
+			ep_lib::MPI_Irecv(&recv[i][0], recv[i].size()*sizeof(T), EP_CHAR, ranks[i], 0, communicator, &request[nbRequest++]);
 	for (int i = 0; i < ranks.size(); i++)
 		if (send[i].size())
-			MPI_Isend((void *) &send[i][0], send[i].size()*sizeof(T), MPI_CHAR, ranks[i], 0, communicator, &request[nbRequest++]);
-	MPI_Waitall(nbRequest, &request[0], &status[0]);
+			ep_lib::MPI_Isend((void *) &send[i][0], send[i].size()*sizeof(T), EP_CHAR, ranks[i], 0, communicator, &request[nbRequest++]);
+	ep_lib::MPI_Waitall(nbRequest, &request[0], &status[0]);
 }
 
@@ -39,8 +38,8 @@
     but message lengths are *unknown* to receiver */
 template <typename T>
-void alltoalls_unknown(const vector<vector<T> >& send, vector<vector<T> >& recv, const vector<int>& ranks, MPI_Comm communicator)
-{
-	vector<MPI_Request> request(ranks.size() * 2);
-	vector<MPI_Status>  status(ranks.size() * 2);
+void alltoalls_unknown(const vector<vector<T> >& send, vector<vector<T> >& recv, const vector<int>& ranks, ep_lib::MPI_Comm communicator)
+{
+	vector<ep_lib::MPI_Request> request(ranks.size() * 2);
+	vector<ep_lib::MPI_Status>  status(ranks.size() * 2);
 
 	// communicate sizes
@@ -51,8 +50,8 @@
 		sendSizes[i] = send[i].size();
 	for (int i = 0; i < ranks.size(); i++)
-		MPI_Irecv(&recvSizes[i], 1, MPI_INT, ranks[i], 0, communicator, &request[nbRequest++]);
+		ep_lib::MPI_Irecv(&recvSizes[i], 1, EP_INT, ranks[i], 0, communicator, &request[nbRequest++]);
 	for (int i = 0; i < ranks.size(); i++)
-		MPI_Isend(&sendSizes[i], 1, MPI_INT, ranks[i], 0, communicator, &request[nbRequest++]);
-	MPI_Waitall(nbRequest, &request[0], &status[0]);
+		ep_lib::MPI_Isend(&sendSizes[i], 1, EP_INT, ranks[i], 0, communicator, &request[nbRequest++]);
+	ep_lib::MPI_Waitall(nbRequest, &request[0], &status[0]);
 
 	// allocate
@@ -119,10 +118,10 @@
 	CTimer::get("CMPIRouting::init(reduce_scatter)").reset();
 	CTimer::get("CMPIRouting::init(reduce_scatter)").resume();
-	MPI_Reduce_scatter(toSend, &nbSource, recvCount, MPI_INT, MPI_SUM, communicator);
+	ep_lib::MPI_Reduce_scatter(toSend, &nbSource, recvCount, EP_INT, EP_SUM, communicator);
 	CTimer::get("CMPIRouting::init(reduce_scatter)").suspend();
 	CTimer::get("CMPIRouting::init(reduce_scatter)").print();
 
-	MPI_Alloc_mem(nbTarget *sizeof(int), MPI_INFO_NULL, &targetRank);
-	MPI_Alloc_mem(nbSource *sizeof(int), MPI_INFO_NULL, &sourceRank);
+	ep_lib::MPI_Alloc_mem(nbTarget *sizeof(int), EP_INFO_NULL, &targetRank);
+	ep_lib::MPI_Alloc_mem(nbSource *sizeof(int), EP_INFO_NULL, &sourceRank);
 
 	targetRankToIndex = new int[mpiSize];
@@ -138,10 +137,10 @@
 	}
 
-	MPI_Barrier(communicator);
+	ep_lib::MPI_Barrier(communicator);
 	CTimer::get("CMPIRouting::init(get_source)").reset();
 	CTimer::get("CMPIRouting::init(get_source)").resume();
 
-	MPI_Request *request = new MPI_Request[nbSource + nbTarget];
-	MPI_Status  *status = new MPI_Status[nbSource + nbTarget];
+	ep_lib::MPI_Request *request = new ep_lib::MPI_Request[nbSource + nbTarget];
+	ep_lib::MPI_Status  *status = new ep_lib::MPI_Status[nbSource + nbTarget];
 
 	int indexRequest = 0;
@@ -151,14 +150,16 @@
 	for (int i = 0; i < nbSource; i++)
 	{
-		#ifdef _usingEP
-		MPI_Irecv(&sourceRank[i], 1, MPI_INT, -2, 0, communicator, &request[indexRequest++]);
-		#else
-		MPI_Irecv(&sourceRank[i], 1, MPI_INT, MPI_ANY_SOURCE, 0, communicator, &request[indexRequest++]);
+		#ifdef _usingMPI
+                ep_lib::MPI_Irecv(&sourceRank[i], 1, EP_INT, MPI_ANY_SOURCE, 0, communicator, &request[indexRequest]);
+		#elif _usingEP
+                ep_lib::MPI_Irecv(&sourceRank[i], 1, EP_INT, -2, 0, communicator, &request[indexRequest]);
 		#endif
+                indexRequest++;
 	}
 	MPI_Barrier(communicator);
 	for (int i = 0; i < nbTarget; i++)
 	{
-		MPI_Isend(&mpiRank, 1, MPI_INT, targetRank[i], 0, communicator, &request[indexRequest++]);
+		ep_lib::MPI_Isend(&mpiRank, 1, EP_INT, targetRank[i], 0, communicator, &request[indexRequest]);
+		indexRequest++;
 	}
 	MPI_Waitall(indexRequest, request, status);
@@ -173,15 +174,15 @@
 	for (int i = 0; i < nbSource; i++)
 	{
-		#ifdef _usingEP
-		MPI_Irecv(&sourceRank[i], 1, MPI_INT, -2, 0, communicator, &request[indexRequest]);
-		#else
-		MPI_Irecv(&sourceRank[i], 1, MPI_INT, MPI_ANY_SOURCE, 0, communicator, &request[indexRequest]);
+		#ifdef _usingMPI
+                ep_lib::MPI_Irecv(&sourceRank[i], 1, EP_INT, MPI_ANY_SOURCE, 0, communicator, &request[indexRequest]);
+		#elif _usingEP
+                ep_lib::MPI_Irecv(&sourceRank[i], 1, EP_INT, -2, 0, communicator, &request[indexRequest]);
 		#endif
-		indexRequest++;
-	}
-
-	for (int i = 0; i < nbTarget; i++)
-	{
-		MPI_Isend(&mpiRank, 1, MPI_INT, targetRank[i], 0, communicator, &request[indexRequest]);
+                indexRequest++;
+	}
+
+	for (int i = 0; i < nbTarget; i++)
+	{
+		ep_lib::MPI_Isend(&mpiRank, 1, EP_INT, targetRank[i], 0, communicator, &request[indexRequest]);
 		indexRequest++;
 	}
@@ -208,5 +209,5 @@
 	for (int i = 0; i < nbSource; i++)
 	{
-		MPI_Irecv(&nbSourceElement[i], 1, MPI_INT, sourceRank[i], 0, communicator, &request[indexRequest]);
+		ep_lib::MPI_Irecv(&nbSourceElement[i], 1, EP_INT, sourceRank[i], 0, communicator, &request[indexRequest]);
 		indexRequest++;
 	}
@@ -215,5 +216,5 @@
 	{
 		totalTargetElement += nbTargetElement[i];
-		MPI_Isend(&nbTargetElement[i], 1, MPI_INT, targetRank[i], 0, communicator, &request[indexRequest]);
+		ep_lib::MPI_Isend(&nbTargetElement[i], 1, EP_INT, targetRank[i], 0, communicator, &request[indexRequest]);
 		indexRequest++;
 	}
@@ -283,9 +284,9 @@
 
 
-	MPI_Request* request=new MPI_Request[nbSource+nbTarget];
-	MPI_Status*  status=new MPI_Status[nbSource+nbTarget];
+	ep_lib::MPI_Request* request=new ep_lib::MPI_Request[nbSource+nbTarget];
+	ep_lib::MPI_Status*  status=new ep_lib::MPI_Status[nbSource+nbTarget];
 	int indexRequest=0;
 
-	MPI_Barrier(communicator);
+	ep_lib::MPI_Barrier(communicator);
 	CTimer::get("CMPIRouting::transferToTarget").reset();
 	CTimer::get("CMPIRouting::transferToTarget").resume();
@@ -293,5 +294,5 @@
 	for(int i=0; i<nbSource; i++)
 	{
-		MPI_Irecv(sourceBuffer[i],nbSourceElement[i]*sizeof(T),MPI_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]);
+		ep_lib::MPI_Irecv(sourceBuffer[i],nbSourceElement[i]*sizeof(T),EP_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]);
 		indexRequest++;
 	}
@@ -299,13 +300,13 @@
 	for(int i=0;i<nbTarget; i++)
 	{
-		MPI_Isend(targetBuffer[i],nbTargetElement[i]*sizeof(T), MPI_CHAR, targetRank[i], 0, communicator, &request[indexRequest]);
-		indexRequest++;
-	}
-
-	MPI_Waitall(indexRequest,request,status);
+		ep_lib::MPI_Isend(targetBuffer[i],nbTargetElement[i]*sizeof(T), EP_CHAR, targetRank[i], 0, communicator, &request[indexRequest]);
+		indexRequest++;
+	}
+
+	ep_lib::MPI_Waitall(indexRequest,request,status);
 
 	CTimer::get("CMPIRouting::transferToTarget").suspend();
 	CTimer::get("CMPIRouting::transferToTarget").print();
-	MPI_Barrier(communicator);
+	ep_lib::MPI_Barrier(communicator);
 
 	// unpack the data
@@ -347,9 +348,9 @@
 	}
 
-	MPI_Request *request = new MPI_Request[nbSource + nbTarget];
-	MPI_Status  *status = new MPI_Status[nbSource + nbTarget];
+	ep_lib::MPI_Request *request = new ep_lib::MPI_Request[nbSource + nbTarget];
+	ep_lib::MPI_Status  *status = new ep_lib::MPI_Status[nbSource + nbTarget];
 	int indexRequest = 0;
 
-	MPI_Barrier(communicator);
+	ep_lib::MPI_Barrier(communicator);
 	CTimer::get("CMPIRouting::transferToTarget(messageSize)").reset();
 	CTimer::get("CMPIRouting::transferToTarget(messageSize)").resume();
@@ -357,5 +358,5 @@
 	for(int i=0; i<nbSource; i++)
 	{
-		MPI_Irecv(&sourceMessageSize[i],1,MPI_INT, sourceRank[i], 0, communicator, &request[indexRequest]);
+		ep_lib::MPI_Irecv(&sourceMessageSize[i],1,EP_INT, sourceRank[i], 0, communicator, &request[indexRequest]);
 		indexRequest++;
 	}
@@ -363,11 +364,11 @@
 	for(int i=0; i<nbTarget; i++)
 	{
-		MPI_Isend(&targetMessageSize[i],1, MPI_INT, targetRank[i], 0, communicator, &request[indexRequest]);
-		indexRequest++;
-	}
-
-	MPI_Waitall(indexRequest,request,status);
-
-	MPI_Barrier(communicator);
+		ep_lib::MPI_Isend(&targetMessageSize[i],1, EP_INT, targetRank[i], 0, communicator, &request[indexRequest]);
+		indexRequest++;
+	}
+
+	ep_lib::MPI_Waitall(indexRequest,request,status);
+
+	ep_lib::MPI_Barrier(communicator);
 	CTimer::get("CMPIRouting::transferToTarget(messageSize)").suspend();
 	CTimer::get("CMPIRouting::transferToTarget(messageSize)").print();
@@ -402,5 +403,5 @@
 	for(int i=0; i<nbSource; i++)
 	{
-		MPI_Irecv(sourceBuffer[i],sourceMessageSize[i],MPI_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]);
+		ep_lib::MPI_Irecv(sourceBuffer[i],sourceMessageSize[i],EP_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]);
 		indexRequest++;
 	}
@@ -408,5 +409,5 @@
 	for(int i=0;i<nbTarget; i++)
 	{
-		MPI_Isend(targetBuffer[i],targetMessageSize[i], MPI_CHAR, targetRank[i], 0, communicator, &request[indexRequest]);
+		ep_lib::MPI_Isend(targetBuffer[i],targetMessageSize[i], EP_CHAR, targetRank[i], 0, communicator, &request[indexRequest]);
 		indexRequest++;
 	}
@@ -467,11 +468,11 @@
 	}
 
-	MPI_Request* request=new MPI_Request[nbSource+nbTarget];
-	MPI_Status*  status=new MPI_Status[nbSource+nbTarget];
+	ep_lib::MPI_Request* request=new ep_lib::MPI_Request[nbSource+nbTarget];
+	ep_lib::MPI_Status*  status=new ep_lib::MPI_Status[nbSource+nbTarget];
 	int indexRequest=0;
 
 	for(int i=0; i<nbSource; i++)
 	{
-		MPI_Isend(sourceBuffer[i],nbSourceElement[i]*sizeof(T),MPI_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]);
+		ep_lib::MPI_Isend(sourceBuffer[i],nbSourceElement[i]*sizeof(T),EP_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]);
 		indexRequest++;
 	}
@@ -479,9 +480,9 @@
 	for(int i=0;i<nbTarget; i++)
 	{
-		MPI_Irecv(targetBuffer[i],nbTargetElement[i]*sizeof(T), MPI_CHAR, targetRank[i], 0, communicator, &request[indexRequest]);
-		indexRequest++;
-	}
-
-	MPI_Waitall(indexRequest,request,status);
+		ep_lib::MPI_Irecv(targetBuffer[i],nbTargetElement[i]*sizeof(T), EP_CHAR, targetRank[i], 0, communicator, &request[indexRequest]);
+		indexRequest++;
+	}
+
+	ep_lib::MPI_Waitall(indexRequest,request,status);
 
 	// unpack the data
@@ -523,18 +524,18 @@
 	}
 
-	MPI_Request *request = new MPI_Request[nbSource + nbTarget];
-	MPI_Status  *status = new MPI_Status[nbSource + nbTarget];
+	ep_lib::MPI_Request *request = new ep_lib::MPI_Request[nbSource + nbTarget];
+	ep_lib::MPI_Status  *status = new ep_lib::MPI_Status[nbSource + nbTarget];
 	int indexRequest = 0;
 	for (int i = 0; i < nbSource; i++)
 	{
-		MPI_Isend(&sourceMessageSize[i], 1, MPI_INT, sourceRank[i], 0, communicator, &request[indexRequest]);
-		indexRequest++;
-	}
-	for (int i = 0; i < nbTarget; i++)
-	{
-		MPI_Irecv(&targetMessageSize[i], 1, MPI_INT, targetRank[i], 0, communicator, &request[indexRequest]);
-		indexRequest++;
-	}
-	MPI_Waitall(indexRequest, request, status);
+		ep_lib::MPI_Isend(&sourceMessageSize[i], 1, EP_INT, sourceRank[i], 0, communicator, &request[indexRequest]);
+		indexRequest++;
+	}
+	for (int i = 0; i < nbTarget; i++)
+	{
+		ep_lib::MPI_Irecv(&targetMessageSize[i], 1, EP_INT, targetRank[i], 0, communicator, &request[indexRequest]);
+		indexRequest++;
+	}
+	ep_lib::MPI_Waitall(indexRequest, request, status);
 
 	for (int i = 0; i < nbTarget; i++)
@@ -564,13 +565,13 @@
 	for (int i = 0; i < nbSource; i++)
 	{
-		MPI_Isend(sourceBuffer[i], sourceMessageSize[i], MPI_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]);
-		indexRequest++;
-	}
-	for (int i = 0; i < nbTarget; i++)
-	{
-		MPI_Irecv(targetBuffer[i], targetMessageSize[i], MPI_CHAR, targetRank[i], 0, communicator, &request[indexRequest]);
-		indexRequest++;
-	}
-	MPI_Waitall(indexRequest, request, status);
+		ep_lib::MPI_Isend(sourceBuffer[i], sourceMessageSize[i], EP_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]);
+		indexRequest++;
+	}
+	for (int i = 0; i < nbTarget; i++)
+	{
+		ep_lib::MPI_Irecv(targetBuffer[i], targetMessageSize[i], EP_CHAR, targetRank[i], 0, communicator, &request[indexRequest]);
+		indexRequest++;
+	}
+	ep_lib::MPI_Waitall(indexRequest, request, status);
 
 	// unpack the data
@@ -612,8 +613,8 @@
 
 template void alltoalls_unknown(const std::vector<std::vector<NES> >& send, std::vector<std::vector<NES> >& recv,
-                                const std::vector<int>& ranks, MPI_Comm communicator);
+                                const std::vector<int>& ranks, ep_lib::MPI_Comm communicator);
 
 template void alltoalls_known(const std::vector<std::vector<int> >& send, std::vector<std::vector<int> >& recv,
-                              const std::vector<int>& ranks, MPI_Comm communicator);
-
-}
+                              const std::vector<int>& ranks, ep_lib::MPI_Comm communicator);
+
+}
Index: XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.cpp
===================================================================
--- XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.cpp	(revision 1545)
+++ XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.cpp	(revision 1642)
@@ -12,13 +12,6 @@
 
 #include "parallel_tree.hpp"
-using namespace ep_lib;
 
 namespace sphereRemap {
-
-extern CRemapGrid srcGrid;
-#pragma omp threadprivate(srcGrid)
-
-extern CRemapGrid tgtGrid;
-#pragma omp threadprivate(tgtGrid)
 
 static const int assignLevel = 2;
@@ -121,9 +114,10 @@
 }
 
-CParallelTree::CParallelTree(MPI_Comm comm) : communicator(comm), cascade(MIN_NODE_SZ*MIN_NODE_SZ*2, comm)
-{
-  treeCascade.reserve(cascade.num_levels);
-  for (int lev = 0; lev < cascade.num_levels; lev++)
-    treeCascade.push_back(CSampleTree(cascade.level[lev].group_size, assignLevel));
+//CParallelTree::CParallelTree(MPI_Comm comm) : communicator(comm), cascade(MIN_NODE_SZ*MIN_NODE_SZ, comm)
+CParallelTree::CParallelTree(ep_lib::MPI_Comm comm) : communicator(comm), cascade(MAX_NODE_SZ*MAX_NODE_SZ*2, comm)
+{
+	treeCascade.reserve(cascade.num_levels);
+	for (int lev = 0; lev < cascade.num_levels; lev++)
+		treeCascade.push_back(CSampleTree(cascade.level[lev].group_size, assignLevel));
 }
 
@@ -157,6 +151,5 @@
 
 	int nrecv; // global number of samples  THIS WILL BE THE NUMBER OF LEAFS IN THE SAMPLE TREE
-	MPI_Allreduce(&n, &nrecv, 1, MPI_INT, MPI_SUM, comm.comm); // => size of sample tree does not depend on keepNodes!
-
+	ep_lib::MPI_Allreduce(&n, &nrecv, 1, EP_INT, EP_SUM, comm.comm); // => size of sample tree does not depend on keepNodes!
 	double ratio = blocSize / (1.0 * nrecv);
 	int nsend = ratio * n + 1; // nsend = n_local_samples / n_global_samples * blocksize + 1 = blocksize/comm.size
@@ -164,5 +157,5 @@
 
 	int *counts = new int[comm.size];
-	MPI_Allgather(&nsend, 1, MPI_INT, counts, 1, MPI_INT, comm.comm);
+	ep_lib::MPI_Allgather(&nsend, 1, EP_INT, counts, 1, EP_INT, comm.comm);
 
 	nrecv = 0;
@@ -190,5 +183,5 @@
 	/* each process needs the sample elements from all processes */
 	double *recvBuffer = new double[nrecv*4];
-	MPI_Allgatherv(sendBuffer, 4 * nsend, MPI_DOUBLE, recvBuffer, counts, displs, MPI_DOUBLE, comm.comm);
+	ep_lib::MPI_Allgatherv(sendBuffer, 4 * nsend, EP_DOUBLE, recvBuffer, counts, displs, EP_DOUBLE, comm.comm);
 	delete[] sendBuffer;
 	delete[] counts;
@@ -248,5 +241,5 @@
          << "   node size : "<<node.size()<<"   bloc size : "<<blocSize<<"  total number of leaf : "<<tree.leafs.size()<<endl ;
 /*
-	MPI_Allreduce(&ok, &allok, 1, MPI_INT, MPI_PROD, communicator);
+	MPI_Allreduce(&ok, &allok, 1, EP_INT, MPI_PROD, communicator);
 	if (!allok) {
 		MPI_Finalize();
@@ -254,5 +247,5 @@
 	}
 */
-    MPI_Abort(MPI_COMM_WORLD,-1) ;
+    ep_lib::MPI_Abort(EP_COMM_WORLD,-1) ;
   }
 /*
@@ -272,5 +265,5 @@
 {
 	CMPIRouting MPIRoute(communicator);
-	MPI_Barrier(communicator);
+	ep_lib::MPI_Barrier(communicator);
 	CTimer::get("buildLocalTree(initRoute)").resume();
 	MPIRoute.init(route);
@@ -297,5 +290,5 @@
 
 	int mpiRank;
-	MPI_Comm_rank(communicator, &mpiRank);
+	ep_lib::MPI_Comm_rank(communicator, &mpiRank);
 	localTree.leafs.reserve(nbLocalElements);
 	for (int i = 0; i < nbLocalElements; i++)
@@ -323,7 +316,7 @@
   nb1=node.size() ; nb2=node2.size() ;
   nb=nb1+nb2 ;
-  MPI_Allreduce(&nb, &nbTot, 1, MPI_LONG, MPI_SUM, communicator) ;
+  ep_lib::MPI_Allreduce(&nb, &nbTot, 1, EP_LONG, EP_SUM, communicator) ;
   int commSize ;
-  MPI_Comm_size(communicator,&commSize) ;
+  ep_lib::MPI_Comm_size(communicator,&commSize) ;
   
 	// make multiple of two
@@ -508,9 +501,9 @@
 	// gather circles on this level of the cascade
 	int pg_size;
-	MPI_Comm_size(cascade.level[level].pg_comm, &pg_size);
+	ep_lib::MPI_Comm_size(cascade.level[level].pg_comm, &pg_size);
 	vector<Coord> allRootCentres(pg_size);
 	vector<double> allRootRadia(pg_size);
-	MPI_Allgather(&rootCentre, 3, MPI_DOUBLE, &allRootCentres[0], 3, MPI_DOUBLE, cascade.level[level].pg_comm);
-	MPI_Allgather(&rootRadius, 1, MPI_DOUBLE, &allRootRadia[0],   1, MPI_DOUBLE, cascade.level[level].pg_comm);
+	ep_lib::MPI_Allgather(&rootCentre, 3, EP_DOUBLE, &allRootCentres[0], 3, EP_DOUBLE, cascade.level[level].pg_comm);
+	ep_lib::MPI_Allgather(&rootRadius, 1, EP_DOUBLE, &allRootRadia[0],   1, EP_DOUBLE, cascade.level[level].pg_comm);
 
 	// now allRootsRadia and allRootCentres must be inserted into second levels of us and propagated to root
Index: XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.hpp
===================================================================
--- XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.hpp	(revision 1545)
+++ XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.hpp	(revision 1642)
@@ -12,26 +12,26 @@
 {
 public:
-  CParallelTree(ep_lib::MPI_Comm comm);
-  ~CParallelTree();
+	CParallelTree(ep_lib::MPI_Comm comm);
+	~CParallelTree();
 
-  void build(vector<Node>& node, vector<Node>& node2);
+	void build(vector<Node>& node, vector<Node>& node2);
 
-  void routeNodes(vector<int>& route, vector<Node>& nodes, int level = 0);
-  void routeIntersections(vector<vector<int> >& route, vector<Node>& nodes, int level = 0);
+	void routeNodes(vector<int>& route, vector<Node>& nodes, int level = 0);
+	void routeIntersections(vector<vector<int> >& route, vector<Node>& nodes, int level = 0);
 
-  int nbLocalElements;
-  Elt* localElements;
+	int nbLocalElements;
+	Elt* localElements;
 
-  CTree localTree;
+	CTree localTree;
 
 private:
-  void updateCirclesForRouting(Coord rootCentre, double rootRadius, int level = 0);
-  void buildSampleTreeCascade(vector<Node>& sampleNodes, int level = 0);
-  void buildLocalTree(const vector<Node>& node, const vector<int>& route);
-  void buildRouteTree();
+	void updateCirclesForRouting(Coord rootCentre, double rootRadius, int level = 0);
+	void buildSampleTreeCascade(vector<Node>& sampleNodes, int level = 0);
+	void buildLocalTree(const vector<Node>& node, const vector<int>& route);
+	void buildRouteTree();
 
-  //CSampleTree sampleTree;
-  vector<CSampleTree> treeCascade; // first for sample tree, then for routing tree
-  CMPICascade cascade;
+	//CSampleTree sampleTree;
+	vector<CSampleTree> treeCascade; // first for sample tree, then for routing tree
+	CMPICascade cascade;
   ep_lib::MPI_Comm communicator ;
 
Index: XIOS/dev/branch_openmp/extern/remap/src/polyg.cpp
===================================================================
--- XIOS/dev/branch_openmp/extern/remap/src/polyg.cpp	(revision 1545)
+++ XIOS/dev/branch_openmp/extern/remap/src/polyg.cpp	(revision 1642)
@@ -20,18 +20,18 @@
 void orient(int N, Coord *vertex, Coord *edge, double *d, const Coord &g)
 {
-	Coord ga = vertex[0] - g;
-	Coord gb = vertex[1] - g;
-	Coord vertical = crossprod(ga, gb);
-	if (N > 2 && scalarprod(g, vertical) < 0)  // (GAxGB).G
-	{
-		for (int i = 0; i < N/2; i++)
-			swap(vertex[i], vertex[N-1-i]);
-
-		for (int i = 0; i < (N-1)/2; i++)
-		{
-			swap(edge[N-2-i], edge[i]);
-			swap(d[i], d[N-2-i]);
-		}
-	}
+  Coord ga = vertex[0] - g;
+  Coord gb = vertex[1] - g;
+  Coord vertical = crossprod(ga, gb);
+  if (N > 2 && scalarprod(g, vertical) < 0)  // (GAxGB).G
+  {
+    for (int i = 0; i < N/2; i++)
+      swap(vertex[i], vertex[N-1-i]);
+
+    for (int i = 0; i < (N-1)/2; i++)
+    {
+      swap(edge[N-2-i], edge[i]);
+      swap(d[i], d[N-2-i]);
+    }
+  }
 }
 
@@ -39,21 +39,21 @@
 void normals(Coord *x, int n, Coord *a)
 {
-	for (int i = 0; i<n; i++)
-		a[i] = crossprod(x[(i+1)%n], x[i]);
+  for (int i = 0; i<n; i++)
+    a[i] = crossprod(x[(i+1)%n], x[i]);
 }
 
 Coord barycentre(const Coord *x, int n)
 {
-	if (n == 0) return ORIGIN;
-	Coord bc = ORIGIN;
-	for (int i = 0; i < n; i++)
-		bc = bc + x[i];
-	/* both distances can be equal down to roundoff when norm(bc) < mashineepsilon 
-	   which can occur when weighted with tiny area */
-
-	assert(squaredist(bc, proj(bc)) <= squaredist(bc, proj(bc * (-1.0))));
-	//if (squaredist(bc, proj(bc)) > squaredist(bc, proj(bc * (-1.0)))) return proj(bc * (-1.0));
-
-	return proj(bc);
+  if (n == 0) return ORIGIN;
+  Coord bc = ORIGIN;
+  for (int i = 0; i < n; i++)
+    bc = bc + x[i];
+  /* both distances can be equal down to roundoff when norm(bc) < mashineepsilon 
+     which can occur when weighted with tiny area */
+
+  assert(squaredist(bc, proj(bc)) <= squaredist(bc, proj(bc * (-1.0))));
+  //if (squaredist(bc, proj(bc)) > squaredist(bc, proj(bc * (-1.0)))) return proj(bc * (-1.0));
+
+  return proj(bc);
 }
 
@@ -62,16 +62,16 @@
 static Coord tetrah_side_diff_centre(Coord a, Coord b)
 {
-	Coord n = crossprod(a,b);
+  Coord n = crossprod(a,b);
         double sinc2 = n.x*n.x + n.y*n.y + n.z*n.z;
-	assert(sinc2 < 1.0 + EPS);
+  assert(sinc2 < 1.0 + EPS);
 
         // exact: u = asin(sinc)/sinc - 1; asin(sinc) = geodesic length of arc ab
-	// approx:
+  // approx:
         // double u = sinc2/6. + (3./40.)*sinc2*sinc2;
 
-	// exact  
-	if (sinc2 > 1.0 - EPS) /* if round-off leads to sinc > 1 asin produces NaN */
-		return n * (M_PI_2 - 1);
-	double sinc = sqrt(sinc2);
+  // exact  
+  if (sinc2 > 1.0 - EPS) /* if round-off leads to sinc > 1 asin produces NaN */
+    return n * (M_PI_2 - 1);
+  double sinc = sqrt(sinc2);
         double u = asin(sinc)/sinc - 1;
 
@@ -82,56 +82,65 @@
 Coord gc_normalintegral(const Coord *x, int n)
 {
-	Coord m = barycentre(x, n);
-	Coord bc = crossprod(x[n-1]-m, x[0]-m) + tetrah_side_diff_centre(x[n-1], x[0]);
-	for (int i = 1; i < n; i++)
-		bc = bc + crossprod(x[i-1]-m, x[i]-m) + tetrah_side_diff_centre(x[i-1], x[i]);
-	return bc*0.5;
+  Coord m = barycentre(x, n);
+  Coord bc = crossprod(x[n-1]-m, x[0]-m) + tetrah_side_diff_centre(x[n-1], x[0]);
+  for (int i = 1; i < n; i++)
+    bc = bc + crossprod(x[i-1]-m, x[i]-m) + tetrah_side_diff_centre(x[i-1], x[i]);
+  return bc*0.5;
 }
 
 Coord exact_barycentre(const Coord *x, int n)
 {
-	if (n >= 3)
-	{
-		return  proj(gc_normalintegral(x, n));
-	}
-	else if (n == 0) return ORIGIN;
-	else if (n == 2) return midpoint(x[0], x[1]);
-	else if (n == 1) return x[0];
+  if (n >= 3)
+  {
+    return  proj(gc_normalintegral(x, n));
+  }
+  else if (n == 0) return ORIGIN;
+  else if (n == 2) return midpoint(x[0], x[1]);
+  else if (n == 1) return x[0];
 }
 
 Coord sc_gc_moon_normalintegral(Coord a, Coord b, Coord pole)
 {
-	double hemisphere = (a.z > 0) ? 1: -1;
-
-	double lat = hemisphere * (M_PI_2 - acos(a.z));
-	double lon1 = atan2(a.y, a.x);
-	double lon2 = atan2(b.y, b.x);
-	double lon_diff = lon2 - lon1;
-
-	// wraparound at lon=-pi=pi
-	if (lon_diff < -M_PI) lon_diff += 2.0*M_PI;
-	else if (lon_diff > M_PI) lon_diff -= 2.0*M_PI;
-
-	// integral of the normal over the surface bound by great arcs a-pole and b-pole and small arc a-b 
-	Coord sc_normalintegral = Coord(0.5*(sin(lon2)-sin(lon1))*(M_PI_2 - lat - 0.5*sin(2.0*lat)),
-	                                0.5*(cos(lon1)-cos(lon2))*(M_PI_2 - lat - 0.5*sin(2.0*lat)),
-	                                hemisphere * lon_diff * 0.25 * (cos(2.0*lat) + 1.0));
-	Coord p = Coord(0,0,hemisphere); // TODO assumes north pole is (0,0,1)
-	Coord t[] = {a, b, p};
-	if (hemisphere < 0) swap(t[0], t[1]);
-	return (sc_normalintegral - gc_normalintegral(t, 3)) * hemisphere;
-}
-
-
-double triarea(Coord& A, Coord& B, Coord& C)
-{
-	double a = ds(B, C);
-	double b = ds(C, A);
-	double c = ds(A, B);
-	double s = 0.5 * (a + b + c);
-	double t = tan(0.5*s) * tan(0.5*(s - a)) * tan(0.5*(s - b)) * tan(0.5*(s - c));
-//	assert(t >= 0);
-  if (t<1e-20) return 0. ;
-	return 4 * atan(sqrt(t));
+  double hemisphere = (a.z > 0) ? 1: -1;
+
+  double lat = hemisphere * (M_PI_2 - acos(a.z));
+  double lon1 = atan2(a.y, a.x);
+  double lon2 = atan2(b.y, b.x);
+  double lon_diff = lon2 - lon1;
+
+  // wraparound at lon=-pi=pi
+  if (lon_diff < -M_PI) lon_diff += 2.0*M_PI;
+  else if (lon_diff > M_PI) lon_diff -= 2.0*M_PI;
+
+  // integral of the normal over the surface bound by great arcs a-pole and b-pole and small arc a-b 
+  Coord sc_normalintegral = Coord(0.5*(sin(lon2)-sin(lon1))*(M_PI_2 - lat - 0.5*sin(2.0*lat)),
+                                  0.5*(cos(lon1)-cos(lon2))*(M_PI_2 - lat - 0.5*sin(2.0*lat)),
+                                  hemisphere * lon_diff * 0.25 * (cos(2.0*lat) + 1.0));
+  Coord p = Coord(0,0,hemisphere); // TODO assumes north pole is (0,0,1)
+  Coord t[] = {a, b, p};
+  if (hemisphere < 0) swap(t[0], t[1]);
+  return (sc_normalintegral - gc_normalintegral(t, 3)) * hemisphere;
+}
+
+
+double triarea(const Coord& A, const Coord& B, const Coord& C)
+{
+  double a = ds(B, C);
+  double b = ds(C, A);
+  double c = ds(A, B);
+  double tmp ;
+
+  if (a<b) { tmp=a ; a=b ; b=tmp ; }
+  if (c > a) { tmp=a ; a=c ; c=b, b=tmp;  }
+  else if (c > b) { tmp=c ; c=b ; b=tmp ; }
+  
+  double s = 0.5 * (a + b + c);
+  double t = tan(0.25*(a+(b+c))) * tan(0.25*(c-(a-b))) * tan(0.25*( c + (a-b) )) * tan(0.25*( a + (b - c)));
+  if (t>0) return 4 * atan(sqrt(t));
+  else
+  {
+    std::cout<<"double triarea(const Coord& A, const Coord& B, const Coord& C) : t < 0 !!! t="<<t<<endl ;
+    return 0 ;
+  }
 }
 
@@ -142,11 +151,11 @@
 double alun(double b, double d)
 {
-	double a  = acos(d);
-	assert(b <= 2 * a);
-	double s  = a + 0.5 * b;
-	double t = tan(0.5 * s) * tan(0.5 * (s - a)) * tan(0.5 * (s - a)) * tan(0.5 * (s - b));
-	double r  = sqrt(1 - d*d);
-	double p  = 2 * asin(sin(0.5*b) / r);
-	return p*(1 - d) - 4*atan(sqrt(t));
+  double a  = acos(d);
+  assert(b <= 2 * a);
+  double s  = a + 0.5 * b;
+  double t = tan(0.5 * s) * tan(0.5 * (s - a)) * tan(0.5 * (s - a)) * tan(0.5 * (s - b));
+  double r  = sqrt(1 - d*d);
+  double p  = 2 * asin(sin(0.5*b) / r);
+  return p*(1 - d) - 4*atan(sqrt(t));
 }
 
@@ -160,81 +169,84 @@
 double airbar(int N, const Coord *x, const Coord *c, double *d, const Coord& pole, Coord& gg)
 {
-	if (N < 3)
-	  return 0; /* polygons with less than three vertices have zero area */
-	Coord t[3];
-	t[0] = barycentre(x, N);
-	Coord *g = new Coord[N];
-	double area = 0;
-	Coord gg_exact = gc_normalintegral(x, N);
-	for (int i = 0; i < N; i++)
-	{
-		/* for "spherical circle segment" sum triangular part and "small moon" and => account for small circle */
-		int ii = (i + 1) % N;
-		t[1] = x[i];
-		t[2] = x[ii];
+  if (N < 3)
+    return 0; /* polygons with less then three vertices have zero area */
+  Coord t[3];
+  t[0] = barycentre(x, N);
+  Coord *g = new Coord[N];
+  double area = 0;
+  Coord gg_exact = gc_normalintegral(x, N);
+  for (int i = 0; i < N; i++)
+  {
+    /* for "spherical circle segment" sum triangular part and "small moon" and => account for small circle */
+    int ii = (i + 1) % N;
+    t[1] = x[i];
+    t[2] = x[ii];
                 double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ;
-		assert(sc >= -1e-10); // Error: tri a l'env (wrong orientation)
-		double area_gc = triarea(t[0], t[1], t[2]);
-		double area_sc_gc_moon = 0;
-		if (d[i]) /* handle small circle case */
-		{
-			Coord m = midpoint(t[1], t[2]);
-			double mext = scalarprod(m, c[i]) - d[i];
-			char sgl = (mext > 0) ? -1 : 1;
-			area_sc_gc_moon = sgl * alun(arcdist(t[1], t[2]), fabs(scalarprod(t[1], pole)));
-			gg_exact = gg_exact + sc_gc_moon_normalintegral(t[1], t[2], pole);
-		}
-		area += area_gc + area_sc_gc_moon; /* for "spherical circle segment" sum triangular part (at) and "small moon" and => account for small circle */
-		g[i] = barycentre(t, 3) * (area_gc + area_sc_gc_moon);
-	}
-	gg = barycentre(g, N);
-	gg_exact = proj(gg_exact);
-	delete[] g;
-	gg = gg_exact;
-	return area;
+    assert(sc >= -1e-10); // Error: tri a l'env (wrong orientation)
+    double area_gc = triarea(t[0], t[1], t[2]);
+    double area_sc_gc_moon = 0;
+    if (d[i]) /* handle small circle case */
+    {
+      Coord m = midpoint(t[1], t[2]);
+      double mext = scalarprod(m, c[i]) - d[i];
+      char sgl = (mext > 0) ? -1 : 1;
+      area_sc_gc_moon = sgl * alun(arcdist(t[1], t[2]), fabs(scalarprod(t[1], pole)));
+      gg_exact = gg_exact + sc_gc_moon_normalintegral(t[1], t[2], pole);
+    }
+    area += area_gc + area_sc_gc_moon; /* for "spherical circle segment" sum triangular part (at) and "small moon" and => account for small circle */
+    g[i] = barycentre(t, 3) * (area_gc + area_sc_gc_moon);
+  }
+  gg = barycentre(g, N);
+  gg_exact = proj(gg_exact);
+  delete[] g;
+  gg = gg_exact;
+  return area;
 }
 
 double polygonarea(Coord *vertices, int N)
 {
-	assert(N >= 3);
-
-	/* compute polygon area as sum of triangles */
-	Coord centre = barycentre(vertices, N);
-	double area = 0;
-	for (int i = 0; i < N; i++)
-		area += triarea(centre, vertices[i], vertices[(i+1)%N]);
-	return area;
+  assert(N >= 3);
+
+  /* compute polygon area as sum of triangles */
+  Coord centre = barycentre(vertices, N);
+  double area = 0;
+  for (int i = 0; i < N; i++)
+    area += triarea(centre, vertices[i], vertices[(i+1)%N]);
+  return area;
 }
 
 int packedPolygonSize(const Elt& e)
 {
-	return sizeof(e.id) + sizeof(e.src_id) + sizeof(e.x) + sizeof(e.val) +
-	       sizeof(e.n) + e.n*(sizeof(double)+sizeof(Coord));
+  return sizeof(e.id) + sizeof(e.src_id) + sizeof(e.x) + sizeof(e.val)  + sizeof(e.given_area)+
+         sizeof(e.n) + e.n*(sizeof(double)+sizeof(Coord));
 }
 
 void packPolygon(const Elt& e, char *buffer, int& pos) 
 {
-	*((GloId *) &(buffer[pos])) = e.id;
-	pos += sizeof(e.id);
-	*((GloId *) &(buffer[pos])) = e.src_id;
-	pos += sizeof(e.src_id);
-
-	*((Coord *) &(buffer[pos])) = e.x;
-	pos += sizeof(e.x);
-
-	*((double*) &(buffer[pos])) = e.val;
-	pos += sizeof(e.val);
-
-	*((int *) & (buffer[pos])) = e.n;
-	pos += sizeof(e.n);
-
-	for (int i = 0; i < e.n; i++)
-	{
-		*((double *) & (buffer[pos])) = e.d[i];
-		pos += sizeof(e.d[i]);
-
-		*((Coord *) &(buffer[pos])) = e.vertex[i];
-		pos += sizeof(e.vertex[i]);
-	} 
+  *((GloId *) &(buffer[pos])) = e.id;
+  pos += sizeof(e.id);
+  *((GloId *) &(buffer[pos])) = e.src_id;
+  pos += sizeof(e.src_id);
+
+  *((Coord *) &(buffer[pos])) = e.x;
+  pos += sizeof(e.x);
+
+  *((double*) &(buffer[pos])) = e.val;
+  pos += sizeof(e.val);
+
+  *((double*) &(buffer[pos])) = e.given_area;
+  pos += sizeof(e.val);
+
+  *((int *) & (buffer[pos])) = e.n;
+  pos += sizeof(e.n);
+
+  for (int i = 0; i < e.n; i++)
+  {
+    *((double *) & (buffer[pos])) = e.d[i];
+    pos += sizeof(e.d[i]);
+
+    *((Coord *) &(buffer[pos])) = e.vertex[i];
+    pos += sizeof(e.vertex[i]);
+  } 
 
 }
@@ -242,31 +254,34 @@
 void unpackPolygon(Elt& e, const char *buffer, int& pos) 
 {
-	e.id = *((GloId *) &(buffer[pos]));
-	pos += sizeof(e.id);
-	e.src_id = *((GloId *) &(buffer[pos]));
-	pos += sizeof(e.src_id);
-
-	e.x = *((Coord *) & (buffer[pos]) );
-	pos += sizeof(e.x);
-
-	e.val = *((double *) & (buffer[pos]));
-	pos += sizeof(double);
-
-	e.n = *((int *) & (buffer[pos]));
-	pos += sizeof(int);
-
-	for (int i = 0; i < e.n; i++)
-	{
-		e.d[i] = *((double *) & (buffer[pos]));
-		pos += sizeof(double);
-
-		e.vertex[i] = *((Coord *) & (buffer[pos]));
-		pos += sizeof(Coord);
-	}
+  e.id = *((GloId *) &(buffer[pos]));
+  pos += sizeof(e.id);
+  e.src_id = *((GloId *) &(buffer[pos]));
+  pos += sizeof(e.src_id);
+
+  e.x = *((Coord *) & (buffer[pos]) );
+  pos += sizeof(e.x);
+
+  e.val = *((double *) & (buffer[pos]));
+  pos += sizeof(double);
+
+  e.given_area = *((double *) & (buffer[pos]));
+  pos += sizeof(double);
+
+  e.n = *((int *) & (buffer[pos]));
+  pos += sizeof(int);
+
+  for (int i = 0; i < e.n; i++)
+  {
+    e.d[i] = *((double *) & (buffer[pos]));
+    pos += sizeof(double);
+
+    e.vertex[i] = *((Coord *) & (buffer[pos]));
+    pos += sizeof(Coord);
+  }
 }
 
 int packIntersectionSize(const Elt& elt) 
 {
-	return elt.is.size() * (2*sizeof(int)+ sizeof(GloId) + 5*sizeof(double));
+  return elt.is.size() * (2*sizeof(int)+ sizeof(GloId) + 5*sizeof(double));
 }
 
@@ -274,60 +289,62 @@
 {
   for (list<Polyg *>::const_iterator it = e.is.begin(); it != e.is.end(); ++it)
-	{
-		*((int *) &(buffer[0])) += 1;
-
-		*((int *) &(buffer[pos])) = e.id.ind;
-		pos += sizeof(int);
+  {
+    *((int *) &(buffer[0])) += 1;
+
+    *((int *) &(buffer[pos])) = e.id.ind;
+    pos += sizeof(int);
 
     *((double *) &(buffer[pos])) = e.area;
     pos += sizeof(double);
 
-		*((GloId *) &(buffer[pos])) = (*it)->id;
-		pos += sizeof(GloId);
+
+    *((GloId *) &(buffer[pos])) = (*it)->id;
+    pos += sizeof(GloId);
   
-		*((int *) &(buffer[pos])) = (*it)->n;
-		pos += sizeof(int);
-		*((double *) &(buffer[pos])) = (*it)->area;
-		pos += sizeof(double);
-
-		*((Coord *) &(buffer[pos])) = (*it)->x;
-		pos += sizeof(Coord);
-	}
+    *((int *) &(buffer[pos])) = (*it)->n;
+    pos += sizeof(int);
+    *((double *) &(buffer[pos])) = (*it)->area;
+    pos += sizeof(double);
+
+    *((Coord *) &(buffer[pos])) = (*it)->x;
+    pos += sizeof(Coord);
+  }
 }
 
 void unpackIntersection(Elt* e, const char* buffer) 
 {
-	int ind;
-	int pos = 0;
+  int ind;
+  int pos = 0;
   
-	int n = *((int *) & (buffer[pos]));
-	pos += sizeof(int);
-	for (int i = 0; i < n; i++)
-	{
-		ind = *((int*) &(buffer[pos]));
-		pos+=sizeof(int);
-
-		Elt& elt= e[ind];
+  int n = *((int *) & (buffer[pos]));
+  pos += sizeof(int);
+  for (int i = 0; i < n; i++)
+  {
+    ind = *((int*) &(buffer[pos]));
+    pos+=sizeof(int);
+
+    Elt& elt= e[ind];
 
     elt.area=*((double *) & (buffer[pos]));
-		pos += sizeof(double);
-
-		Polyg *polygon = new Polyg;
-
-		polygon->id =  *((GloId *) & (buffer[pos]));
-		pos += sizeof(GloId);
-
-		polygon->n =  *((int *) & (buffer[pos]));
-		pos += sizeof(int);
-
-		polygon->area =  *((double *) & (buffer[pos]));
-		pos += sizeof(double);
-
-		polygon->x = *((Coord *) & (buffer[pos]));
-		pos += sizeof(Coord);
-
-		elt.is.push_back(polygon);
-	}
-}
-
-}
+    pos += sizeof(double);
+
+
+    Polyg *polygon = new Polyg;
+
+    polygon->id =  *((GloId *) & (buffer[pos]));
+    pos += sizeof(GloId);
+
+    polygon->n =  *((int *) & (buffer[pos]));
+    pos += sizeof(int);
+
+    polygon->area =  *((double *) & (buffer[pos]));
+    pos += sizeof(double);
+
+    polygon->x = *((Coord *) & (buffer[pos]));
+    pos += sizeof(Coord);
+
+    elt.is.push_back(polygon);
+  }
+}
+
+}
Index: XIOS/dev/branch_openmp/extern/remap/src/polyg.hpp
===================================================================
--- XIOS/dev/branch_openmp/extern/remap/src/polyg.hpp	(revision 1545)
+++ XIOS/dev/branch_openmp/extern/remap/src/polyg.hpp	(revision 1642)
@@ -23,4 +23,5 @@
 void unpackIntersection(Elt *e, const char *buffer);
 int packIntersectionSize(const Elt& e);
+double triarea( const Coord& A, const Coord& B, const Coord& C) ;
 
 }
Index: XIOS/dev/branch_openmp/extern/remap/src/timerRemap.cpp
===================================================================
--- XIOS/dev/branch_openmp/extern/remap/src/timerRemap.cpp	(revision 1545)
+++ XIOS/dev/branch_openmp/extern/remap/src/timerRemap.cpp	(revision 1642)
@@ -4,5 +4,4 @@
 #include <map>
 #include <iostream>
-using namespace ep_lib;
 
 namespace sphereRemap {
@@ -10,6 +9,5 @@
 using namespace std;
 
-//map<string,CTimer*> CTimer::allTimer;
-map<string,CTimer*> *CTimer::allTimer_ptr = 0;
+map<string,CTimer*> CTimer::allTimer;
 
 CTimer::CTimer(const string& name_) : name(name_)
@@ -58,9 +56,6 @@
 {
 	map<string,CTimer*>::iterator it;
-	if(allTimer_ptr == 0) allTimer_ptr = new map<string,CTimer*>;
-	//it=allTimer.find(name);
-	it=allTimer_ptr->find(name);
-	//if (it==allTimer.end()) it=allTimer.insert(pair<string,CTimer*>(name,new CTimer(name))).first;
-	if (it==allTimer_ptr->end()) it=allTimer_ptr->insert(pair<string,CTimer*>(name,new CTimer(name))).first;
+	it=allTimer.find(name);
+	if (it==allTimer.end()) it=allTimer.insert(pair<string,CTimer*>(name,new CTimer(name))).first;
 	return *(it->second);
 }
Index: XIOS/dev/branch_openmp/extern/remap/src/timerRemap.hpp
===================================================================
--- XIOS/dev/branch_openmp/extern/remap/src/timerRemap.hpp	(revision 1545)
+++ XIOS/dev/branch_openmp/extern/remap/src/timerRemap.hpp	(revision 1642)
@@ -26,7 +26,5 @@
     double getCumulatedTime(void);
     void print(void);
-    //static map<string,CTimer*> allTimer;
-    static map<string,CTimer*> *allTimer_ptr;
-    #pragma omp threadprivate(allTimer_ptr)
+    static map<string,CTimer*> allTimer;
     static double getTime(void);
     static CTimer& get(string name);
