diff --git a/DocumentationInterface b/DocumentationInterface
index 28ae7d02cf129055ba7f7f4efdac2117da73c0d9..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 100644
--- a/DocumentationInterface
+++ b/DocumentationInterface
@@ -1,268 +0,0 @@
-initatmgrid(rank_bn,nbcore,polygons_in,centre_in,area_in,contfrac_in,neighbours_in,[nbpt,nbseg])
-
-Wrapper for ``initatmgrid``.
-
-Parameters
-----------
-rank_bn : input int
-nbcore : input int
-polygons_in : input rank-3 array('f') with bounds (nbpt,2 * nbseg + 1,2)
-centre_in : input rank-2 array('f') with bounds (nbpt,2)
-area_in : input rank-1 array('f') with bounds (nbpt)
-contfrac_in : input rank-1 array('f') with bounds (nbpt)
-neighbours_in : input rank-2 array('i') with bounds (nbpt,2 * nbseg)
-
-Other Parameters
-----------------
-nbpt : input int, optional
-    Default: shape(polygons_in,0)
-nbseg : input int, optional
-    Default: (shape(polygons_in,1)-1)/(2)
-====================================================================
-nbi,nbj,area_bx,trip_bx,basin_bx,topoind_bx,fac_bx,hierarchy_bx,lon_bx,lat_bx,lshead_bx,nwbas = gethydrogrid(nbxmax_in,sub_pts,sub_index,sub_area,max_basins,min_topoind,lon_rel,lat_rel,trip,basins,topoindex,fac,hierarchy,[nbpt,nbvmax_in])
-
-Wrapper for ``gethydrogrid``.
-
-Parameters
-----------
-nbxmax_in : input int
-sub_pts : input rank-1 array('i') with bounds (nbpt)
-sub_index : input rank-3 array('i') with bounds (nbpt,nbvmax_in,2)
-sub_area : input rank-2 array('f') with bounds (nbpt,nbvmax_in)
-max_basins : in/output rank-0 array(float,'f')
-min_topoind : input float
-lon_rel : input rank-2 array('f') with bounds (nbpt,nbvmax_in)
-lat_rel : input rank-2 array('f') with bounds (nbpt,nbvmax_in)
-trip : in/output rank-2 array('f') with bounds (nbpt,nbvmax_in)
-basins : in/output rank-2 array('f') with bounds (nbpt,nbvmax_in)
-topoindex : in/output rank-2 array('f') with bounds (nbpt,nbvmax_in)
-fac : in/output rank-2 array('f') with bounds (nbpt,nbvmax_in)
-hierarchy : in/output rank-2 array('f') with bounds (nbpt,nbvmax_in)
-
-Other Parameters
-----------------
-nbpt : input int, optional
-    Default: len(sub_pts)
-nbvmax_in : input int, optional
-    Default: shape(sub_index,1)
-
-Returns
--------
-nbi : rank-1 array('i') with bounds (nbpt)
-nbj : rank-1 array('i') with bounds (nbpt)
-area_bx : rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
-trip_bx : rank-3 array('i') with bounds (nbpt,nbxmax_in,nbxmax_in)
-basin_bx : rank-3 array('i') with bounds (nbpt,nbxmax_in,nbxmax_in)
-topoind_bx : rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
-fac_bx : rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
-hierarchy_bx : rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
-lon_bx : rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
-lat_bx : rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
-lshead_bx : rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
-nwbas : int
-====================================================================
-nb_basin,basin_inbxid,basin_outlet,basin_outtp,basin_sz,basin_bxout,basin_bbout,basin_pts,basin_lshead,coast_pts = findbasins(nbvmax_in,nbi,nbj,trip_bx,basin_bx,fac_bx,hierarchy_bx,topoind_bx,lshead_bx,lontmp,lattmp,[nbpt,nbxmax_in])
-
-Wrapper for ``findbasins``.
-
-Parameters
-----------
-nbvmax_in : input int
-nbi : input rank-1 array('i') with bounds (nbpt)
-nbj : input rank-1 array('i') with bounds (nbpt)
-trip_bx : in/output rank-3 array('i') with bounds (nbpt,nbxmax_in,nbxmax_in)
-basin_bx : in/output rank-3 array('i') with bounds (nbpt,nbxmax_in,nbxmax_in)
-fac_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
-hierarchy_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
-topoind_bx : in/output rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
-lshead_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
-lontmp : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
-lattmp : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
-
-Other Parameters
-----------------
-nbpt : input int, optional
-    Default: len(nbi)
-nbxmax_in : input int, optional
-    Default: shape(trip_bx,1)
-
-Returns
--------
-nb_basin : rank-1 array('i') with bounds (nbpt)
-basin_inbxid : rank-2 array('i') with bounds (nbpt,nbvmax_in)
-basin_outlet : rank-3 array('f') with bounds (nbpt,nbvmax_in,2)
-basin_outtp : rank-2 array('f') with bounds (nbpt,nbvmax_in)
-basin_sz : rank-2 array('i') with bounds (nbpt,nbvmax_in)
-basin_bxout : rank-2 array('i') with bounds (nbpt,nbvmax_in)
-basin_bbout : rank-2 array('i') with bounds (nbpt,nbvmax_in)
-basin_pts : rank-4 array('i') with bounds (nbpt,nbvmax_in,nbvmax_in,2)
-basin_lshead : rank-2 array('f') with bounds (nbpt,nbvmax_in)
-coast_pts : rank-2 array('i') with bounds (nbpt,nbvmax_in)
-====================================================================
-basin_count,basin_notrun,basin_area,basin_cg,basin_hierarchy,basin_fac,basin_topoind,basin_id,basin_coor,basin_type,basin_flowdir,basin_lshead,outflow_grid,outflow_basin,nbcoastal,coastal_basin = globalize(area_bx,lon_bx,lat_bx,trip_bx,hierarchy_bx,fac_bx,topoind_bx,min_topoind,nb_basin,basin_inbxid,basin_outlet,basin_outtp,basin_sz,basin_pts,basin_bxout,basin_bbout,lshead,coast_pts,nwbas,[nbpt,nbvmax_in,nbxmax_in])
-
-Wrapper for ``globalize``.
-
-Parameters
-----------
-area_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
-lon_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
-lat_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
-trip_bx : input rank-3 array('i') with bounds (nbpt,nbxmax_in,nbxmax_in)
-hierarchy_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
-fac_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
-topoind_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
-min_topoind : input float
-nb_basin : input rank-1 array('i') with bounds (nbpt)
-basin_inbxid : input rank-2 array('i') with bounds (nbpt,nbvmax_in)
-basin_outlet : input rank-3 array('f') with bounds (nbpt,nbvmax_in,2)
-basin_outtp : input rank-2 array('f') with bounds (nbpt,nbvmax_in)
-basin_sz : input rank-2 array('i') with bounds (nbpt,nbvmax_in)
-basin_pts : input rank-4 array('i') with bounds (nbpt,nbvmax_in,nbvmax_in,2)
-basin_bxout : input rank-2 array('i') with bounds (nbpt,nbvmax_in)
-basin_bbout : input rank-2 array('i') with bounds (nbpt,nbvmax_in)
-lshead : input rank-2 array('f') with bounds (nbpt,nbvmax_in)
-coast_pts : input rank-2 array('i') with bounds (nbpt,nbvmax_in)
-nwbas : input int
-
-Other Parameters
-----------------
-nbpt : input int, optional
-    Default: shape(area_bx,0)
-nbvmax_in : input int, optional
-    Default: shape(basin_inbxid,1)
-nbxmax_in : input int, optional
-    Default: shape(area_bx,1)
-
-Returns
--------
-basin_count : rank-1 array('i') with bounds (nbpt)
-basin_notrun : rank-1 array('i') with bounds (nbpt)
-basin_area : rank-2 array('f') with bounds (nbpt,nwbas)
-basin_cg : rank-3 array('f') with bounds (nbpt,nwbas,2)
-basin_hierarchy : rank-2 array('f') with bounds (nbpt,nwbas)
-basin_fac : rank-2 array('f') with bounds (nbpt,nwbas)
-basin_topoind : rank-2 array('f') with bounds (nbpt,nwbas)
-basin_id : rank-2 array('i') with bounds (nbpt,nwbas)
-basin_coor : rank-3 array('f') with bounds (nbpt,nwbas,2)
-basin_type : rank-2 array('f') with bounds (nbpt,nwbas)
-basin_flowdir : rank-2 array('i') with bounds (nbpt,nwbas)
-basin_lshead : rank-2 array('f') with bounds (nbpt,nwbas)
-outflow_grid : rank-2 array('i') with bounds (nbpt,nwbas)
-outflow_basin : rank-2 array('i') with bounds (nbpt,nwbas)
-nbcoastal : rank-1 array('i') with bounds (nbpt)
-coastal_basin : rank-2 array('i') with bounds (nbpt,nwbas)
-====================================================================
-inflow_number,inflow_grid,inflow_basin = linkup(nbxmax_in,basin_count,basin_area,basin_id,basin_flowdir,basin_lshead,basin_hierarchy,basin_fac,outflow_grid,outflow_basin,nbcoastal,coastal_basin,invented_basins,[nbpt,nwbas])
-
-Wrapper for ``linkup``.
-
-Parameters
-----------
-nbxmax_in : input int
-basin_count : input rank-1 array('i') with bounds (nbpt)
-basin_area : input rank-2 array('f') with bounds (nbpt,nwbas)
-basin_id : input rank-2 array('i') with bounds (nbpt,nwbas)
-basin_flowdir : input rank-2 array('i') with bounds (nbpt,nwbas)
-basin_lshead : input rank-2 array('f') with bounds (nbpt,nwbas)
-basin_hierarchy : in/output rank-2 array('f') with bounds (nbpt,nwbas)
-basin_fac : input rank-2 array('f') with bounds (nbpt,nwbas)
-outflow_grid : in/output rank-2 array('i') with bounds (nbpt,nwbas)
-outflow_basin : in/output rank-2 array('i') with bounds (nbpt,nwbas)
-nbcoastal : in/output rank-1 array('i') with bounds (nbpt)
-coastal_basin : in/output rank-2 array('i') with bounds (nbpt,nwbas)
-invented_basins : input float
-
-Other Parameters
-----------------
-nbpt : input int, optional
-    Default: len(basin_count)
-nwbas : input int, optional
-    Default: shape(basin_area,1)
-
-Returns
--------
-inflow_number : rank-2 array('i') with bounds (nbpt,8 * nbxmax_in)
-inflow_grid : rank-3 array('i') with bounds (nbpt,8 * nbxmax_in,8 * nbxmax_in)
-inflow_basin : rank-3 array('i') with bounds (nbpt,8 * nbxmax_in,8 * nbxmax_in)
-====================================================================
-outflow_uparea = fetch(corepts,basin_count,basin_area,basin_id,basin_hierarchy,basin_fac,outflow_grid,outflow_basin,partial_sum,fetch_basin,[nbpt,nwbas,nbcore])
-
-Wrapper for ``fetch``.
-
-Parameters
-----------
-corepts : input rank-1 array('i') with bounds (nbcore)
-basin_count : input rank-1 array('i') with bounds (nbpt)
-basin_area : input rank-2 array('f') with bounds (nbpt,nwbas)
-basin_id : input rank-2 array('i') with bounds (nbpt,nwbas)
-basin_hierarchy : input rank-2 array('f') with bounds (nbpt,nwbas)
-basin_fac : input rank-2 array('f') with bounds (nbpt,nwbas)
-outflow_grid : input rank-2 array('i') with bounds (nbpt,nwbas)
-outflow_basin : input rank-2 array('i') with bounds (nbpt,nwbas)
-partial_sum : input rank-2 array('f') with bounds (nbpt,nwbas)
-fetch_basin : in/output rank-2 array('f') with bounds (nbpt,nwbas)
-
-Other Parameters
-----------------
-nbpt : input int, optional
-    Default: len(basin_count)
-nwbas : input int, optional
-    Default: shape(basin_area,1)
-nbcore : input int, optional
-    Default: len(corepts)
-
-Returns
--------
-outflow_uparea : rank-1 array('f') with bounds (nbpt*nwbas)
-====================================================================
-routing_area,routing_cg,topo_resid,route_nbbasin,route_togrid,route_tobasin,route_nbintobas,global_basinid,route_outlet,route_type,origin_nbintobas,routing_fetch = truncate(nbasmax,num_largest,corepts,basin_count,basin_notrun,basin_area,basin_cg,basin_topoind,fetch_basin,basin_id,basin_coor,basin_type,basin_flowdir,outflow_grid,outflow_basin,inflow_number,inflow_grid,inflow_basin,[nbpt,nbxmax_in,nwbas,nbcore])
-
-Wrapper for ``truncate``.
-
-Parameters
-----------
-nbasmax : input int
-num_largest : input int
-corepts : input rank-1 array('i') with bounds (nbcore)
-basin_count : in/output rank-1 array('i') with bounds (nbpt)
-basin_notrun : in/output rank-1 array('i') with bounds (nbpt)
-basin_area : in/output rank-2 array('f') with bounds (nbpt,nwbas)
-basin_cg : in/output rank-3 array('f') with bounds (nbpt,nwbas,2)
-basin_topoind : in/output rank-2 array('f') with bounds (nbpt,nwbas)
-fetch_basin : in/output rank-2 array('f') with bounds (nbpt,nwbas)
-basin_id : in/output rank-2 array('i') with bounds (nbpt,nwbas)
-basin_coor : in/output rank-3 array('f') with bounds (nbpt,nwbas,2)
-basin_type : in/output rank-2 array('f') with bounds (nbpt,nwbas)
-basin_flowdir : in/output rank-2 array('i') with bounds (nbpt,nwbas)
-outflow_grid : in/output rank-2 array('i') with bounds (nbpt,nwbas)
-outflow_basin : in/output rank-2 array('i') with bounds (nbpt,nwbas)
-inflow_number : in/output rank-2 array('i') with bounds (nbpt,8 * nbxmax_in)
-inflow_grid : in/output rank-3 array('i') with bounds (nbpt,8 * nbxmax_in,8 * nbxmax_in)
-inflow_basin : in/output rank-3 array('i') with bounds (nbpt,8 * nbxmax_in,8 * nbxmax_in)
-
-Other Parameters
-----------------
-nbpt : input int, optional
-    Default: len(basin_count)
-nbxmax_in : input int, optional
-    Default: (shape(inflow_number,1))/(8)
-nwbas : input int, optional
-    Default: shape(basin_area,1)
-nbcore : input int, optional
-    Default: len(corepts)
-
-Returns
--------
-routing_area : rank-2 array('f') with bounds (nbpt,nbasmax)
-routing_cg : rank-3 array('f') with bounds (nbpt,nbasmax,2)
-topo_resid : rank-2 array('f') with bounds (nbpt,nbasmax)
-route_nbbasin : rank-1 array('i') with bounds (nbpt)
-route_togrid : rank-2 array('i') with bounds (nbpt,nbasmax)
-route_tobasin : rank-2 array('i') with bounds (nbpt,nbasmax)
-route_nbintobas : rank-1 array('i') with bounds (nbpt)
-global_basinid : rank-2 array('i') with bounds (nbpt,nbasmax)
-route_outlet : rank-3 array('f') with bounds (nbpt,nbasmax,2)
-route_type : rank-2 array('f') with bounds (nbpt,nbasmax)
-origin_nbintobas : rank-1 array('i') with bounds (nbpt)
-routing_fetch : rank-2 array('f') with bounds (nbpt,nbasmax)
diff --git a/Interface.py b/Interface.py
index 2e8639133dab68d340f61448c2243b55f117e32b..ebda3639c299db3840918b3d4b344ba9d2172a77 100644
--- a/Interface.py
+++ b/Interface.py
@@ -269,6 +269,9 @@ class HydroGraph :
         self.num_largest = routing_interface.finalrivclass(part.landcorelist, self.route_togrid, self.route_tobasin, self.routing_fetch, \
                                                            hydrosuper.largest_rivarea)
         #
+        nbpt = hydrosuper.basin_count.shape[0]
+        self.nbpt_proc = np.arange(1,nbpt+1)
+        self.proc = np.full(nbpt, part.rank)
         return
     #
     def dumpnetcdf(self, filename, globalgrid, procgrid, part) :
@@ -349,6 +352,32 @@ class HydroGraph :
         else :
             area = np.zeros((1,1))
         area[:,:] = part.gather(areas[:,:])
+        #
+        # Environment
+        # nbpt_proc
+        subpt = procgrid.landscatter(self.nbpt_proc[:], order='F')
+        subpt = subpt.astype(vtyp, copy=False)
+        subpt[np.isnan(subpt)] = NCFillValue
+        if part.rank == 0 :
+            subptgrid = outnf.createVariable("nbpt_proc", vtyp, ('y','x'), fill_value=NCFillValue)
+            subptgrid.title = "gridpoint reference inside each proc"
+            subptgrid.units = "-"
+        else :
+            subptgrid = np.zeros((1,1,1))
+        subptgrid[:,:] = part.gather(subpt)
+        #
+        # rank
+        procnum = procgrid.landscatter(self.proc[:], order='F')
+        procnum = procnum.astype(vtyp, copy=False)
+        procnum[np.isnan(procnum)] = NCFillValue
+        if part.rank == 0 :
+            procn = outnf.createVariable("proc_num", vtyp, ('y','x'), fill_value=NCFillValue)
+            procn.title = "rank"
+            procn.units = "-"
+        else :
+            procn = np.zeros((1,1,1))
+        procn[:,:] = part.gather(procnum)
+
         #
         # Variables
         # Once gathered on root_proc we transform them into float64. Careful, Integer variables do not have NaN ! 
diff --git a/RoutingPreProc.py b/RoutingPreProc.py
index bdf36b56aa5d84a11e407f26af8bcb6396128977..eb56fe1dc0752664a84c366d63f69dd86015c1dd 100644
--- a/RoutingPreProc.py
+++ b/RoutingPreProc.py
@@ -128,12 +128,14 @@ print("nbasmax : ", nbasmax)
 #
 # Extract hydo data from file
 #
+INFO("hydrodata")
 hydrodata = HG.HydroData(hydrogrid.ncfile, hydrogrid.box, sub_index)
 
+INFO("nitiatmgrid")
 IF.initatmgrid(rank, nbcore, nbpt, modelgrid)
 
 
-
+INFO("hoverlap")
 hoverlap = IF.HydroOverlap(nbpt, nbvmax, sub_pts, sub_index, sub_area, sub_lon, sub_lat, modelgrid, hydrodata)
 
 #
@@ -146,7 +148,7 @@ del sub_lon
 del sub_lat
 gc.collect()
 
-if nbcore == 1 :
+if rank ==0 :
     if lonint[0] != lonint[1] and latint[0] != latint[1] :
         DP.MAPPLOT("MapHydroGrid", lonint, latint, hoverlap, hoverlap.hierarchy_bx, modelgrid.polyll, title="Distance to ocean")
 
@@ -154,13 +156,13 @@ hsuper = IF.HydroSuper(nbvmax, hydrodata, hoverlap)
 #
 # Plot the hydrological supermesh
 #
-if nbcore == 1 :
+if rank == 0 :
     if lonint[0] != lonint[1] and latint[0] != latint[1] :
         DP.SUPERMESHPLOT("MapSuperGrid_Beforelinkup", lonint, latint, hoverlap, hsuper, modelgrid.polyll)
 
 print("=================== LINKUP ====================")
 hsuper.linkup(hydrodata)
-if nbcore == 1 :
+if rank ==0 :
     if lonint[0] != lonint[1] and latint[0] != latint[1] :
         DP.SUPERMESHPLOT("MapSuperGrid_Afterlinkup", lonint, latint, hoverlap, hsuper, modelgrid.polyll)