From e8810ecec5da20a9de032cda0e1a0677e755d22b Mon Sep 17 00:00:00 2001
From: POLCHER Jan <jan.polcher@lmd.jussieu.fr>
Date: Fri, 10 Jan 2020 00:23:55 +0100
Subject: [PATCH] Added the function to dump the graph of HydroSuper, that is
 before the truncation.

---
 F90subroutines/routing_interface.f90 |   3 +-
 Interface.py                         | 347 +++++++++++++++++----------
 RoutingPreProc.py                    |   1 +
 3 files changed, 225 insertions(+), 126 deletions(-)

diff --git a/F90subroutines/routing_interface.f90 b/F90subroutines/routing_interface.f90
index e3d6ba7..ebbacd6 100644
--- a/F90subroutines/routing_interface.f90
+++ b/F90subroutines/routing_interface.f90
@@ -154,7 +154,7 @@ SUBROUTINE findbasins(nbpt, nbvmax_in, nbxmax_in, nbi, nbj, trip_bx, basin_bx, f
   !
   INTEGER :: ib
   !
-  diaglalo(1,:) = (/ 38.25, -7.25 /)
+  diaglalo(1,:) = (/ 39.6791, 2.6687 /)
   !
   IF ( nbvmax_in .NE. nbvmax .OR. nbxmax_in .NE. nbxmax ) THEN
      WRITE(*,*) "nbvmax or nbxmax have changed !!"
@@ -508,6 +508,7 @@ SUBROUTINE finalfetch(nbpt, nbasmax, nbcore, corepts, routing_area, basin_count,
   !
   USE defprec
   USE ioipsl
+  USE constantes_var
   !
   !! INPUT VARIABLES
   INTEGER(i_std), INTENT(in)   :: nbpt !! Domain size (unitless)
diff --git a/Interface.py b/Interface.py
index ebda363..60df07c 100644
--- a/Interface.py
+++ b/Interface.py
@@ -71,6 +71,106 @@ def closeinterface(comm) :
 #
 #
 #
+def addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorners, cornerind) :
+    #
+    # Longitude
+    longitude = part.gather(procgrid.lon_full)
+    if part.rank == 0 :
+        lon=outnf.createVariable("lon", vtyp, ('y','x'), fill_value=NCFillValue)
+        lon.units="grid box centre degrees east"
+        lon.title="Longitude"
+        lon.axis="X"
+        lon[:,:] = longitude[:,:]
+    #
+    # Latitude
+    latitude = part.gather(procgrid.lat_full)
+    if part.rank == 0 :
+        lat=outnf.createVariable("lat", vtyp, ('y','x'), fill_value=NCFillValue)
+        lat.units="grid box centre degrees north"
+        lat.standard_name="grid latitude"
+        lat.title="Latitude"
+        lat.axis="Y"
+        lat[:] = latitude[:,:]
+    #
+    # Bounds of grid box
+    #
+    llonpoly=np.zeros((nbcorners,procgrid.nbland))
+    llatpoly=np.zeros((nbcorners,procgrid.nbland))
+    for i in range(procgrid.nbland) :
+        llonpoly[:,i] = [procgrid.polyll[i][ic][0] for ic in cornerind]
+        llatpoly[:,i] = [procgrid.polyll[i][ic][1] for ic in cornerind]
+        lon_bnd = procgrid.landscatter(llonpoly)
+        lat_bnd = procgrid.landscatter(llatpoly)
+    if part.rank == 0 :
+        lonbnd=outnf.createVariable("lon_bnd", vtyp, ('bnd','y','x'), fill_value=NCFillValue)
+        lonbnd.units="grid box corners degrees east"
+        lonbnd.title="Longitude of Corners"
+        latbnd=outnf.createVariable("lat_bnd", vtyp, ('bnd','y','x'), fill_value=NCFillValue)
+        latbnd.units="grid box corners degrees north"
+        latbnd.title="Latitude of Corners"
+    else :
+        lonbnd= np.zeros((1,1,1))
+        latbnd= np.zeros((1,1,1))
+    lonbnd[:,:,:] = part.gather(lon_bnd[:,:,:])
+    latbnd[:,:,:] = part.gather(lat_bnd[:,:,:])
+    #
+    # Land sea mask
+    #
+    if part.rank == 0 :
+        land=outnf.createVariable("land", vtyp, ('y','x'), fill_value=NCFillValue)
+        land.units="Land Sea mask"
+        land.standard_name="landsea mask"
+        land.title="Land"
+        land[:,:] = globalgrid.land[:,:]
+    # Area
+    areas = procgrid.landscatter(np.array(procgrid.area, dtype=np.float64))
+    areas[np.isnan(areas)] = NCFillValue
+    if part.rank == 0 :
+        area=outnf.createVariable("area", vtyp, ('y','x'), fill_value=NCFillValue)
+        area.units="m^2"
+        area.standard_name="grid area"
+        area.title="Area"
+    else :
+        area = np.zeros((1,1))
+    area[:,:] = part.gather(areas[:,:])
+    #
+    return
+#
+# Add environment to netCDF file
+#
+def addenvironment(outnf, procgrid, part, vtyp, NCFillValue, nbpt) :
+    #
+    nbpt_proc = np.arange(1,nbpt+1, dtype=vtyp)
+    proc = np.full(nbpt, part.rank, dtype=vtyp)
+    # Environment
+    # nbpt_proc
+    subpt = procgrid.landscatter(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))
+    subptgrid[:,:] = part.gather(subpt)
+    #
+    # rank
+    procnum = procgrid.landscatter(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))
+    procn[:,:] = part.gather(procnum)
+    #
+    return
+#
+#
+#
 def finalfetch(part, routing_area, basin_count, route_togrid, route_tobasin, fetch_in) :
     #
     fetch_out = np.zeros(routing_area.shape, dtype=np.float32, order='F')
@@ -169,6 +269,10 @@ class HydroOverlap :
 #
 class HydroSuper :
     def __init__(self, nbvmax, hydrodata, hydrooverlap) :
+        #
+        # Keep largest possible number of HTUs
+        #
+        self.nbhtuext = nbvmax
         #
         # Call findbasins
         #
@@ -191,6 +295,7 @@ class HydroSuper :
                                                 hydrooverlap.hierarchy_bx, hydrooverlap.fac_bx, hydrooverlap.topoind_bx, hydrodata.topoindmin, \
                                                 nb_basin, basin_inbxid, basin_outlet, basin_outtp, self.basin_sz, self.basin_pts, basin_bxout, \
                                                 basin_bbout, basin_lshead, coast_pts, hydrooverlap.nwbas)
+        self.nbpt = self.basin_count.shape[0]
         return
     #
     def linkup(self, hydrodata) :
@@ -249,12 +354,115 @@ class HydroSuper :
                                                                self.fetch_basin, self.largest_rivarea)
         print("Rank :"+str(part.rank)+" Area of smallest large rivers : ", self.largest_rivarea, " Nb of Large rivers on proc : ",self.num_largest)
         return
+    #
+    #
+    #
+    def dumpnetcdf(self, filename, globalgrid, procgrid, part) :
+        #
+        NCFillValue=1.0e20
+        vtyp=np.float64
+        cornerind=[0,2,4,6]
+        nbcorners = len(cornerind)
+        #
+        if part.rank == 0 :
+            outnf=Dataset(filename, 'w', format='NETCDF4_CLASSIC')
+            # Dimensions
+            outnf.createDimension('x', globalgrid.ni)
+            outnf.createDimension('y', globalgrid.nj)
+            outnf.createDimension('land', len(procgrid.area))
+            outnf.createDimension('htuext', self.nbhtuext)
+            outnf.createDimension('bnd', nbcorners)
+        else :
+            outnf = None
+        #
+        addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorners, cornerind)
+        addenvironment(outnf, procgrid, part, vtyp, NCFillValue, self.nbpt)
+        # 
+        # Variables
+        # Once gathered on root_proc we transform them into float64. Careful, Integer variables do not have NaN !
+        #
+        # self.basin_id
+        bid = procgrid.landscatter(self.basin_id[:,:].astype(vtyp), order='F')
+        bid[np.isnan(bid)] = NCFillValue
+        if part.rank == 0 :
+            basinid = outnf.createVariable("HTUid", vtyp, ('htuext','y','x'), fill_value=NCFillValue)
+            basinid.title = "ID for each HTU"
+            basinid.units = "-"
+        else :
+            basinid = np.zeros((1,1,1))
+        basinid[:,:,:] = part.gather(bid)
+        #
+        #self.basin_count
+        bcnt = procgrid.landscatter(self.basin_count[:].astype(vtyp), order='F')
+        bcnt[np.isnan(bcnt)] = NCFillValue
+        if part.rank == 0 :
+            basincnt = outnf.createVariable("HTUcount", vtyp, ('y','x'), fill_value=NCFillValue)
+            basincnt.title = "HTU count"
+            basincnt.units = "-"
+        else :
+            basincnt = np.zeros((1,1))
+        basincnt[:,:] = part.gather(bcnt)
+        #
+        #self.basin_area
+        ba = procgrid.landscatter(self.basin_area[:,:].astype(vtyp), order='F')
+        ba[np.isnan(ba)] = NCFillValue
+        if part.rank == 0 :
+            basinarea = outnf.createVariable("HTUarea", vtyp, ('htuext','y','x'), fill_value=NCFillValue)
+            basinarea.title = "HTU area"
+            basinarea.units = "m^2"
+        else :
+            basinarea = np.zeros((1,1,1))
+        basinarea[:,:,:] = part.gather(ba)
+        #
+        #self.outflow_grid
+        og = procgrid.landscatter(self.outflow_grid[:,:].astype(vtyp), order='F')
+        og[np.isnan(og)] = NCFillValue
+        if part.rank == 0 :
+            outgrid = outnf.createVariable("HTUoutgrid", vtyp, ('htuext','y','x'), fill_value=NCFillValue)
+            outgrid.title = "HTU outflow grid"
+            outgrid.units = "-"
+        else :
+            outgrid = np.zeros((1,1,1))
+        outgrid[:,:,:] = part.gather(og)
+        #
+        #self.outflow_basin
+        ob = procgrid.landscatter(self.outflow_basin[:,:].astype(vtyp), order='F')
+        ob[np.isnan(ob)] = NCFillValue
+        if part.rank == 0 :
+            outbas = outnf.createVariable("HTUoutbasin", vtyp, ('htuext','y','x'), fill_value=NCFillValue)
+            outbas.title = "Outflow HTU of grid"
+            outbas.units = "-"
+        else :
+            outbas = np.zeros((1,1,1))
+        outbas[:,:,:] = part.gather(ob)
+        #
+        # Save the fetch of each basin
+        #
+        fe = procgrid.landscatter(self.fetch_basin[:,:].astype(vtyp), order='F')
+        fe[np.isnan(fe)] = NCFillValue
+        if part.rank == 0 :
+            fetch = outnf.createVariable("fetch", vtyp, ('htuext','y','x'), fill_value=NCFillValue)
+            fetch.title = "Fetch contributing to each HTU"
+            fetch.units = "m^2"
+        else :
+            fetch = np.zeros((1,1,1))
+        fetch[:,:,:] = part.gather(fe)
+        #
+        # Close file
+        #
+        if part.rank == 0 :
+            outnf.close()
+        #
+        return
 #
 #
 #
 class HydroGraph :
     def __init__(self, nbasmax, hydrosuper, part) :
+        #
         self.nbasmax = nbasmax
+        self.nbpt = hydrosuper.basin_count.shape[0]
+        #
         self.routing_area, self.routing_cg, self.topo_resid, self.route_nbbasin, self.route_togrid, self.route_tobasin, self.route_nbintobas, \
             self.global_basinid, self.route_outlet, self.route_type, self.origin_nbintobas, self.routing_fetch = \
                                     routing_interface.truncate(nbasmax, hydrosuper.num_largest, part.landcorelist, hydrosuper.basin_count, \
@@ -269,9 +477,6 @@ 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) :
@@ -289,109 +494,11 @@ class HydroGraph :
             outnf.createDimension('land', len(procgrid.area))
             outnf.createDimension('htu', self.nbasmax)
             outnf.createDimension('bnd', nbcorners)
-        #
-        # Coordinates
-        #
-        # Longitude
-        longitude = part.gather(procgrid.lon_full)
-        if part.rank == 0 :
-            lon=outnf.createVariable("lon", vtyp, ('y','x'), fill_value=NCFillValue)
-            lon.units="grid box centre degrees east"
-            lon.title="Longitude"
-            lon.axis="X"
-            lon[:,:] = longitude[:,:]
-        #
-        # Latitude
-        latitude = part.gather(procgrid.lat_full)
-        if part.rank == 0 :
-            lat=outnf.createVariable("lat", vtyp, ('y','x'), fill_value=NCFillValue)
-            lat.units="grid box centre degrees north"
-            lat.standard_name="grid latitude"
-            lat.title="Latitude"
-            lat.axis="Y"
-            lat[:] = latitude[:,:]
-        #
-        # Bounds of grid box
-        #
-        llonpoly=np.zeros((nbcorners,procgrid.nbland))
-        llatpoly=np.zeros((nbcorners,procgrid.nbland))
-        for i in range(procgrid.nbland) :
-            llonpoly[:,i] = [procgrid.polyll[i][ic][0] for ic in cornerind]
-            llatpoly[:,i] = [procgrid.polyll[i][ic][1] for ic in cornerind]
-        lon_bnd = procgrid.landscatter(llonpoly)
-        lat_bnd = procgrid.landscatter(llatpoly)
-        if part.rank == 0 :
-            lonbnd=outnf.createVariable("lon_bnd", vtyp, ('bnd','y','x'), fill_value=NCFillValue)
-            lonbnd.units="grid box corners degrees east"
-            lonbnd.title="Longitude of Corners"
-            latbnd=outnf.createVariable("lat_bnd", vtyp, ('bnd','y','x'), fill_value=NCFillValue)
-            latbnd.units="grid box corners degrees north"
-            latbnd.title="Latitude of Corners"
-        else :
-            lonbnd= np.zeros((1,1,1))
-            latbnd= np.zeros((1,1,1))
-        lonbnd[:,:,:] = part.gather(lon_bnd[:,:,:])
-        latbnd[:,:,:] = part.gather(lat_bnd[:,:,:])
-        #
-        # Land sea mask
-        #
-        if part.rank == 0 :
-                land=outnf.createVariable("land", vtyp, ('y','x'), fill_value=NCFillValue)
-                land.units="Land Sea mask"
-                land.standard_name="landsea mask"
-                land.title="Land"
-                land[:,:] = globalgrid.land[:,:]
-        # Area
-        areas = procgrid.landscatter(np.array(procgrid.area, dtype=np.float64))
-        areas[np.isnan(areas)] = NCFillValue
-        if part.rank == 0 :
-            area=outnf.createVariable("area", vtyp, ('y','x'), fill_value=NCFillValue)
-            area.units="m^2"
-            area.standard_name="grid area"
-            area.title="Area"
-        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 ! 
+            outnf = None
         #
-        rarea = procgrid.landscatter(self.routing_area[:,:], order='F')
-        rarea = rarea.astype(vtyp, copy=False)
-        rarea[np.isnan(rarea)] = NCFillValue
-        if part.rank == 0 :
-            routingarea = outnf.createVariable("routingarea", vtyp, ('htu','y','x'), fill_value=NCFillValue)
-            routingarea.title = "Surface of basin"
-            routingarea.units = "m^2"
-        else :
-            routingarea = np.zeros((1,1,1))
-        routingarea[:,:,:] = part.gather(rarea)
+        addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorners, cornerind)
+        addenvironment(outnf, procgrid, part, vtyp, NCFillValue, self.nbpt)
         #
         # The field route_togrid is with indices on the local grid. That needs to be converted to the global grid.
         #
@@ -408,8 +515,7 @@ class HydroGraph :
         grgrid = part.l2glandindex(self.route_togrid[:,:])
         if itarget >+ 0 :
             print(part.rank, itarget, " After route_togrid = ", self.route_togrid[itarget,:])
-        rgrid = procgrid.landscatter(grgrid, order='F')
-        rgrid = rgrid.astype(vtyp, copy=False)
+        rgrid = procgrid.landscatter(grgrid.astype(vtyp), order='F')
         rgrid[rgrid >= RPP.IntFillValue] = NCFillValue
         if part.rank == 0 :
             routetogrid = outnf.createVariable("routetogrid", vtyp, ('htu','y','x'), fill_value=NCFillValue)
@@ -419,7 +525,7 @@ class HydroGraph :
             routetogrid = np.zeros((1,1,1))    
         routetogrid[:,:,:] = part.gather(rgrid)
         #
-        rtobasin = procgrid.landscatter(self.route_tobasin[:,:], order='F')
+        rtobasin = procgrid.landscatter(self.route_tobasin[:,:].astype(vtyp), order='F')
         rtobasin = rtobasin.astype(vtyp, copy=False)
         rtobasin[rtobasin >= RPP.IntFillValue] = NCFillValue
         if part.rank == 0 :
@@ -430,8 +536,7 @@ class HydroGraph :
             routetobasin = np.zeros((1,1,1))
         routetobasin[:,:,:] = part.gather(rtobasin)
         #
-        rid = procgrid.landscatter(self.global_basinid[:,:], order='F')
-        rid = rid.astype(vtyp, copy=False)
+        rid = procgrid.landscatter(self.global_basinid[:,:].astype(vtyp), order='F')
         rid[rid >= RPP.IntFillValue] = NCFillValue
         if part.rank == 0 :                           
             basinid = outnf.createVariable("basinid", vtyp, ('htu','y','x'), fill_value=NCFillValue)
@@ -441,8 +546,7 @@ class HydroGraph :
             basinid = np.zeros((1,1,1))
         basinid[:,:,:] = part.gather(rid)
         #
-        rintobas = procgrid.landscatter(self.route_nbintobas[:])
-        rintobas = rintobas.astype(vtyp, copy=False)
+        rintobas = procgrid.landscatter(self.route_nbintobas[:].astype(vtyp))
         rintobas[rintobas >= RPP.IntFillValue] = NCFillValue
         if part.rank == 0 : 
             routenbintobas = outnf.createVariable("routenbintobas", vtyp, ('y','x'), fill_value=NCFillValue)
@@ -452,8 +556,7 @@ class HydroGraph :
             routenbintobas = np.zeros((1,1))
         routenbintobas[:,:] = part.gather(rintobas)
         #
-        onbintobas = procgrid.landscatter(self.origin_nbintobas[:])
-        onbintobas = onbintobas.astype(vtyp, copy=False)
+        onbintobas = procgrid.landscatter(self.origin_nbintobas[:].astype(vtyp))
         onbintobas[onbintobas >= RPP.IntFillValue] = NCFillValue
         if part.rank == 0 :
             originnbintobas = outnf.createVariable("originnbintobas", vtyp, ('y','x'), fill_value=NCFillValue)
@@ -463,8 +566,7 @@ class HydroGraph :
             originnbintobas = np.zeros((1,1))
         originnbintobas[:,:] = part.gather(onbintobas)
         #
-        olat = procgrid.landscatter(self.route_outlet[:,:,0], order='F')
-        olat = olat.astype(vtyp, copy=False)
+        olat = procgrid.landscatter(self.route_outlet[:,:,0].astype(vtyp), order='F')
         olat[np.isnan(olat)] = NCFillValue
         if part.rank == 0 :
             outletlat = outnf.createVariable("outletlat", vtyp, ('htu','y','x'), fill_value=NCFillValue)
@@ -474,8 +576,7 @@ class HydroGraph :
             outletlat = np.zeros((1,1,1))
         outletlat[:,:,:] = part.gather(olat)
         #
-        olon = procgrid.landscatter(self.route_outlet[:,:,1], order='F')
-        olon = olon.astype(vtyp, copy=False)
+        olon = procgrid.landscatter(self.route_outlet[:,:,1].astype(vtyp), order='F')
         olon[np.isnan(olon)] = NCFillValue
         if part.rank == 0 :
             outletlon = outnf.createVariable("outletlon", vtyp, ('htu','y','x'), fill_value=NCFillValue)
@@ -485,8 +586,7 @@ class HydroGraph :
             outletlon = np.zeros((1,1,1))
         outletlon[:,:,:] = part.gather(olon)
         #
-        otype = procgrid.landscatter(self.route_type[:,:], order='F')
-        otype = otype.astype(vtyp, copy=False)
+        otype = procgrid.landscatter(self.route_type[:,:].astype(vtyp), order='F')
         otype[np.isnan(otype)] = NCFillValue
         if part.rank == 0 :
             outlettype = outnf.createVariable("outlettype", vtyp, ('htu','y','x'), fill_value=NCFillValue)
@@ -496,8 +596,7 @@ class HydroGraph :
             outlettype = np.zeros((1,1,1))
         outlettype[:,:,:] = part.gather(otype)
         #
-        tind = procgrid.landscatter(self.topo_resid[:,:], order='F')
-        tind = tind.astype(vtyp, copy=False)
+        tind = procgrid.landscatter(self.topo_resid[:,:].astype(vtyp), order='F')
         tind[np.isnan(tind)] = NCFillValue
         if part.rank == 0 :
             topoindex = outnf.createVariable("topoindex", vtyp, ('htu','y','x'), fill_value=NCFillValue)
@@ -509,8 +608,7 @@ class HydroGraph :
         #
         # Save centre of gravity of HTU
         #
-        cg = procgrid.landscatter(self.routing_cg[:,:,:], order='F')
-        cg = cg.astype(vtyp, copy=False)
+        cg = procgrid.landscatter(self.routing_cg[:,:,:].astype(vtyp), order='F')
         cg[np.isnan(cg)] = NCFillValue
         if part.rank == 0 :
             CG_lon = outnf.createVariable("CG_lon", vtyp, ('htu','y','x'), fill_value=NCFillValue)
@@ -527,8 +625,7 @@ class HydroGraph :
         #
         # Save the fetch of each basin
         #
-        fe =  procgrid.landscatter(self.routing_fetch[:,:], order='F')
-        fe = fe.astype(vtyp, copy=False)
+        fe =  procgrid.landscatter(self.routing_fetch[:,:].astype(vtyp), order='F')
         fe[np.isnan(fe)] = NCFillValue
         if part.rank == 0 :
             fetch = outnf.createVariable("fetch", vtyp, ('htu','y','x'), fill_value=NCFillValue)
diff --git a/RoutingPreProc.py b/RoutingPreProc.py
index eb56fe1..699318b 100644
--- a/RoutingPreProc.py
+++ b/RoutingPreProc.py
@@ -175,6 +175,7 @@ gc.collect()
 
 print("=================== Compute fetch ====================")
 hsuper.fetch(part)
+hsuper.dumpnetcdf(OutGraphFile.replace(".nc","_HydroSuper.nc"), gg, modelgrid, part)
 
 hgraph = IF.HydroGraph(nbasmax, hsuper, part)
 
-- 
GitLab