Simplification of Interface.py and creating of the GraphHydro.py module to...

Simplification of Interface.py and creating of the GraphHydro.py module to manage the operations on the HTU graph.

GraphHydro.py

+#
+import os
+import sys
+from inspect import currentframe, getframeinfo
+localdir=os.path.dirname(getframeinfo(currentframe()).filename)
+sys.path.append(localdir+'/F90subroutines')
+#
+import numpy as np
+from netCDF4 import Dataset
+#
+import routing_interface
+import NcHelp as NH
+import RPPtools as RPP
+#
+# Add time constants and possibly other prameters to the graph file
+#
+def addparameters(outnf, part, tcst, vtyp, NCFillValue) :
+    if part.rank == 0 :
+        outnf.createDimension('dimpara', 1)
+        stream = outnf.createVariable("StreamTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
+        stream.title = "Time constant for the stream reservoir"
+        stream.units = "day/m"
+        stream[:] = tcst.stream_tcst
+        fast = outnf.createVariable("FastTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
+        fast.title = "Time constant for the fast reservoir"
+        fast.units = "day/m"
+        fast[:] = tcst.fast_tcst
+        slow = outnf.createVariable("SlowTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
+        slow.title = "Time constant for the slow reservoir"
+        slow.units = "day/m"
+        slow[:] = tcst.slow_tcst
+        flood = outnf.createVariable("FloodTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
+        flood.title = "Time constant for the flood reservoir"
+        flood.units = "day/m"
+        flood[:] = tcst.flood_tcst
+        swamp = outnf.createVariable("SwampCst", vtyp, ('dimpara'), fill_value=NCFillValue)
+        swamp.title = "Fraction of the river transport that flows to the swamps"
+        swamp.units = "-"
+        swamp[:] = tcst.swamp_cst
+    return
+#
+#
+#
+class HydroGraph :
+    def __init__(self, nbasmax, hydrosuper, part, modelgrid) :
+        #
+        self.nbasmax = nbasmax
+        self.nbpt = hydrosuper.basin_count.shape[0]
+        nwbas = hydrosuper.basin_topoind.shape[1]
+        nbxmax_in = hydrosuper.inflow_grid.shape[2]
+        #
+        #
+        self.routing_area, self.routing_orog_mean, self.routing_orog_min,self.routing_orog_max, \
+            self.routing_floodp, 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.finish_truncate(nbpt = self.nbpt, inflowmax = nbxmax_in, nbasmax = nbasmax, nwbas = nwbas, \
+                                                                      num_largest = hydrosuper.num_largest, gridarea = modelgrid.area, \
+                                                                      cfrac = modelgrid.contfrac, gridcenters = np.array(modelgrid.centers), \
+                                                                      basin_count = hydrosuper.basin_count, \
+                                                                      basin_notrun = hydrosuper.basin_notrun, basin_area = hydrosuper.basin_area, \
+                                                                      basin_orog_mean = hydrosuper.basin_orog_mean, basin_orog_min = hydrosuper.basin_orog_min,\
+                                                                      basin_orog_max = hydrosuper.basin_orog_max, basin_floodp = hydrosuper.basin_floodp, \
+                                                                      basin_cg = hydrosuper.basin_cg, \
+                                                                      basin_topoind = hydrosuper.basin_topoind, fetch_basin = hydrosuper.fetch_basin, \
+                                                                      basin_id = hydrosuper.basin_id, \
+                                                                      basin_coor = hydrosuper.basin_outcoor, basin_type = hydrosuper.basin_type, \
+                                                                      basin_flowdir = hydrosuper.basin_flowdir, \
+                                                                      outflow_grid = hydrosuper.outflow_grid, outflow_basin = hydrosuper.outflow_basin, \
+                                                                      inflow_number = hydrosuper.inflow_number, inflow_grid = hydrosuper.inflow_grid, \
+                                                                      inflow_basin = hydrosuper.inflow_basin)
+
+        #
+        #
+        #
+        self.num_largest = routing_interface.finalrivclass(part.landcorelist, self.route_togrid, self.route_tobasin, self.routing_fetch, \
+                      hydrosuper.largest_rivarea)
+        #
+        # Inflows
+        self.max_inflow = part.domainmax(np.max(hydrosuper.inflow_number))
+        gingrid = part.l2glandindex( hydrosuper.inflow_grid[:,:,:self.max_inflow])
+        self.route_innum, self.route_ingrid, self.route_inbasin = \
+            routing_interface.finish_inflows(nbpt = self.nbpt, nwbas = nwbas, nbasmax = nbasmax, inf_max = self.max_inflow, \
+                                             basin_count = hydrosuper.basin_count, inflow_number = hydrosuper.inflow_number, \
+                                             inflow_grid = gingrid, inflow_basin = hydrosuper.inflow_basin[:,:,:self.max_inflow])
+
+        return
+    #
+    #
+    #
+    def add_variable(self, outnf, procgrid, NCFillValue, part, coord, name, title, units, data, vtyp, orig_type = "float"):
+        var = procgrid.landscatter(data.astype(vtyp), order='F')
+
+        if orig_type == "float":
+            var[np.isnan(var)] = NCFillValue
+        elif orig_type == "int":
+            var[np.isnan(var)] = RPP.IntFillValue
+            var[var==RPP.IntFillValue] = NCFillValue
+
+        if part.rank == 0:
+            ncvar = outnf.createVariable(name, vtyp, coord, fill_value=NCFillValue)
+            ncvar.title = title
+            ncvar.units = units
+        else :
+            ncvar = np.zeros((1,1,1))
+        ncvar[:] = part.gather(var)
+    #
+    #
+    #
+    def dumpnetcdf(self, filename, globalgrid, procgrid, part, tcst) :
+        #
+        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('z', self.nbasmax)
+            outnf.createDimension('bnd', nbcorners)
+            outnf.createDimension('inflow', self.max_inflow)
+        else :
+            outnf = None
+        #
+        NH.addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorners, cornerind)
+        NH.addenvironment(outnf, procgrid, part, vtyp, NCFillValue, self.nbpt)
+        addparameters(outnf, part, tcst, vtyp, NCFillValue)
+        #
+        # land grid index -> to facilitate the analyses of the routing
+        #
+        nbpt_loc = np.zeros((self.nbpt,1)).astype(np.int32)
+        nbpt_loc[:,0] = np.arange(1, self.nbpt+1)
+        nbpt_glo = part.l2glandindex(nbpt_loc)
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), "nbpt_glo", "Land point indices on global grid", "-", nbpt_glo[:,0], vtyp)
+        #
+        ################
+        #
+        # Conversion of grid indices for route_togrid
+        grgrid = part.l2glandindex(self.route_togrid[:,:])
+        #
+        #
+        # The field route_togrid is with indices on the local grid. That needs to be converted to the global grid.
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "routetogrid", "Grid into which the basin flows", "-", grgrid, vtyp, "int")
+        # route to basin
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "routetobasin", "Basin in to which the water goes", "-", self.route_tobasin[:,:], vtyp, "int")
+        # basin id
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basinid", "ID of basin", "-", self.global_basinid[:,:], vtyp, "int")
+        #
+        # basin area
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basin_area", "area of basin", "m^2", self.routing_area[:,:], vtyp, "float")
+
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basin_orog_mean", "Mean orography", "m", self.routing_orog_mean[:,:], vtyp, "float")
+
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basin_orog_min", "Min orography", "m", self.routing_orog_min[:,:], vtyp, "float")
+
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basin_orog_max", "Max orography", "m", self.routing_orog_max[:,:], vtyp, "float")
+
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basin_floodp", "Fraction of floodplains", "-", self.routing_floodp[:,:], vtyp, "float")
+
+        # route number into basin
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), "routenbintobas", "Number of basin into current one", "-", self.route_nbintobas[:], vtyp, "int")
+        #
+        # original number into basin
+        self.add_variable(outnf, procgrid, NCFillValue, part, ( 'y','x'), "originnbintobas", "Number of sub-grid basin into current one before truncation", "-", self.origin_nbintobas[:], vtyp, "int")
+        #
+        # latitude of outlet
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "outletlat", "Latitude of Outlet", "degrees north", self.route_outlet[:,:,0], vtyp, "float")
+        # longitude of outlet
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "outletlon", "Longitude of Outlet", "degrees east", self.route_outlet[:,:,1], vtyp, "float")
+        # type of outlet
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "outlettype", "Type of outlet", "code", self.route_type[:,:], vtyp, "float")
+        #
+        # topographic index
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "topoindex", "Topographic index of the retention time", "m", self.topo_resid[:,:], vtyp, "float")
+        #
+
+        # Inflow number
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "route_innum", "Number of inflow", "-", self.route_innum[:,:], vtyp, "int")
+        #
+        # Inflow grid
+        #gingrid = part.l2glandindex(self.inflow_grid[:,:,:inflow_size])
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow', 'z','y','x'), "route_ingrid", "Grid from which the water flows", "-", self.route_ingrid[:,:,:], vtyp, "int")
+        #
+        # Inflow basin
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow', 'z','y','x'), "route_inbasin", "Basin from which the water flows", "-", self.route_inbasin[:,:,:], vtyp, "int")
+
+        #
+        # Save centre of gravity of HTU
+        #
+        # Check if it works
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "CG_lon", "Longitude of centre of gravity of HTU", "degrees east", self.routing_cg[:,:,1], vtyp, "float")
+
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "CG_lat", "Latitude of centre of gravity of HTU", "degrees north", self.routing_cg[:,:,0], vtyp, "float")
+        #
+        # Save the fetch of each basin
+        #
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "fetch", "Fetch contributing to each HTU", "m^2", self.routing_fetch[:,:], vtyp, "float")
+        #
+        # Close the file
+        if part.rank == 0:
+            outnf.close()
+        #
+        return

Interface.py

@@ -4,7 +4,6 @@ import numpy as np
 import os
 import pickle
 from netCDF4 import Dataset
-import RPPtools as RPP
 from mpi4py import MPI
 import gc
 #
@@ -24,6 +23,7 @@ else :
 MPI.COMM_WORLD.Barrier()
 #
 import routing_interface
+import NcHelp as NH
 #
 import getargs
 config = getargs.SetupConfig()
@@ -78,183 +78,6 @@ 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[:,:] = globalgrid.lonmat[:,:]
-    #
-    # 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[:,:] = globalgrid.latmat[:,:]
-    #
-    # 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
-#
-# Add time constants and possibly other prameters to the graph file
-#
-def addparameters(outnf, part, tcst, vtyp, NCFillValue) :
-    if part.rank == 0 :
-        outnf.createDimension('dimpara', 1)
-        stream = outnf.createVariable("StreamTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
-        stream.title = "Time constant for the stream reservoir"
-        stream.units = "day/m"
-        stream[:] = tcst.stream_tcst
-        fast = outnf.createVariable("FastTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
-        fast.title = "Time constant for the fast reservoir"
-        fast.units = "day/m"
-        fast[:] = tcst.fast_tcst
-        slow = outnf.createVariable("SlowTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
-        slow.title = "Time constant for the slow reservoir"
-        slow.units = "day/m"
-        slow[:] = tcst.slow_tcst
-        flood = outnf.createVariable("FloodTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
-        flood.title = "Time constant for the flood reservoir"
-        flood.units = "day/m"
-        flood[:] = tcst.flood_tcst
-        swamp = outnf.createVariable("SwampCst", vtyp, ('dimpara'), fill_value=NCFillValue)
-        swamp.title = "Fraction of the river transport that flows to the swamps"
-        swamp.units = "-"
-        swamp[:] = tcst.swamp_cst
-    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')
-    partial_sum = np.zeros(routing_area.shape, dtype=np.float32, order='F')
-    old_sorted = np.zeros(largest_pos, dtype=np.float32, order='F')
-    #
-    nbpt = routing_area.shape[0]
-    fhalolist = [i+1 for i in range(nbpt) if i not in part.landcorelist]
-    in_core = np.array([1 if i in part.landcorelist else 0 for i in range(nbpt)])
-    #
-    maxdiff_sorted = prec*prec
-    iter_count = 0
-    #
-    while iter_count < part.size*3 and maxdiff_sorted > prec :
-        fetch_out[:,:] = 0.0
-        outflow_uparea = routing_interface.finalfetch(in_core,fhalolist, part.landcorelist, routing_area, basin_count, route_togrid, \
-                                                      route_tobasin, partial_sum, fetch_out)
-        partial_sum = np.copy(fetch_out)
-        part.landsendtohalo(partial_sum, order='F')
-        partial_sum = part.zerocore(partial_sum, order='F')
-        #
-        # Find area the largest basins we need to have right.
-        #
-        xtmp = np.hstack(part.comm.allgather(outflow_uparea[np.where(outflow_uparea > 0.0)]))
-        # Precision in m^2 of the upstream areas when sorting.
-        sorted_outareas = (np.unique(np.rint(np.array(xtmp)/prec))*prec)[::-1]
-        if sorted_outareas.shape[0]<largest_pos:
-           s = sorted_outareas[:]
-           sorted_outareas = np.zeros(largest_pos, dtype=np.float32, order='F')
-           sorted_outareas[:s.shape[0]] = s[:]
-        # If mono-proc no need to iterate as fetch produces the full result.
-        if part.size == 1 :
-            maxdiff_sorted = 0.0
-        else :
-            maxdiff_sorted = np.max(np.abs(sorted_outareas[0:largest_pos]-old_sorted))
-            old_sorted[:] = sorted_outareas[0:largest_pos]
-        iter_count += 1
-
-    #
-    fetch_error = np.sum(np.abs(fetch_out[part.landcorelist,:]-fetch_in[part.landcorelist,:]), axis=1)\
-                                                    / np.ma.sum(routing_area[part.landcorelist,:], axis=1)
-    if np.max(fetch_error) > prec :
-        print("Rank :"+str(part.rank)+" Too large fetch error (fraction of greid area) : ", fetch_error)
-
-    print("Total fetch error in fraction of grid box : ", np.sum(fetch_error))
-    #
-    return fetch_out
-#
-#
-#
 class HydroOverlap :
 #
     def __init__(self, nbpt, nbvmax, sub_pts, sub_index_in, sub_area_in, sub_lon_in, sub_lat_in, part, modelgrid, hydrodata) :
@@ -564,8 +387,8 @@ class HydroSuper :
         else :
             outnf = None
         #
-        addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorners, cornerind)
-        addenvironment(outnf, procgrid, part, vtyp, NCFillValue, self.nbpt)
+        NH.addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorners, cornerind)
+        NH.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 !
@@ -650,168 +473,3 @@ class HydroSuper :
 #
 #
 #
-class HydroGraph :
-    def __init__(self, nbasmax, hydrosuper, part, modelgrid) :
-        #
-        self.nbasmax = nbasmax
-        self.nbpt = hydrosuper.basin_count.shape[0]
-        nwbas = hydrosuper.basin_topoind.shape[1]
-        nbxmax_in = hydrosuper.inflow_grid.shape[2]
-        #
-        #
-        self.routing_area, self.routing_orog_mean, self.routing_orog_min,self.routing_orog_max, \
-            self.routing_floodp, 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.finish_truncate(nbpt = self.nbpt, inflowmax = nbxmax_in, nbasmax = nbasmax, nwbas = nwbas, \
-                                                                      num_largest = hydrosuper.num_largest, gridarea = modelgrid.area, \
-                                                                      cfrac = modelgrid.contfrac, gridcenters = np.array(modelgrid.centers), \
-                                                                      basin_count = hydrosuper.basin_count, \
-                                                                      basin_notrun = hydrosuper.basin_notrun, basin_area = hydrosuper.basin_area, \
-                                                                      basin_orog_mean = hydrosuper.basin_orog_mean, basin_orog_min = hydrosuper.basin_orog_min,\
-                                                                      basin_orog_max = hydrosuper.basin_orog_max, basin_floodp = hydrosuper.basin_floodp, \
-                                                                      basin_cg = hydrosuper.basin_cg, \
-                                                                      basin_topoind = hydrosuper.basin_topoind, fetch_basin = hydrosuper.fetch_basin, \
-                                                                      basin_id = hydrosuper.basin_id, \
-                                                                      basin_coor = hydrosuper.basin_outcoor, basin_type = hydrosuper.basin_type, \
-                                                                      basin_flowdir = hydrosuper.basin_flowdir, \
-                                                                      outflow_grid = hydrosuper.outflow_grid, outflow_basin = hydrosuper.outflow_basin, \
-                                                                      inflow_number = hydrosuper.inflow_number, inflow_grid = hydrosuper.inflow_grid, \
-                                                                      inflow_basin = hydrosuper.inflow_basin)
-
-        #
-        # This step is no more necessary
-        #self.routing_fetch = finalfetch(part, self.routing_area, self.route_nbbasin, self.route_togrid, self.route_tobasin, self.routing_fetch)
-        #
-        self.num_largest = routing_interface.finalrivclass(part.landcorelist, self.route_togrid, self.route_tobasin, self.routing_fetch, \
-                      hydrosuper.largest_rivarea)
-        #
-        # Inflows
-        self.max_inflow = part.domainmax(np.max(hydrosuper.inflow_number))
-        gingrid = part.l2glandindex( hydrosuper.inflow_grid[:,:,:self.max_inflow])
-        self.route_innum, self.route_ingrid, self.route_inbasin = \
-            routing_interface.finish_inflows(nbpt = self.nbpt, nwbas = nwbas, nbasmax = nbasmax, inf_max = self.max_inflow, \
-                                             basin_count = hydrosuper.basin_count, inflow_number = hydrosuper.inflow_number, \
-                                             inflow_grid = gingrid, inflow_basin = hydrosuper.inflow_basin[:,:,:self.max_inflow])
-
-        return
-    #
-    #
-    #
-    def add_variable(self, outnf, procgrid, NCFillValue, part, coord, name, title, units, data, vtyp, orig_type = "float"):
-        var = procgrid.landscatter(data.astype(vtyp), order='F')
-
-        if orig_type == "float":
-            var[np.isnan(var)] = NCFillValue
-        elif orig_type == "int":
-            var[np.isnan(var)] = RPP.IntFillValue
-            var[var==RPP.IntFillValue] = NCFillValue
-
-        if part.rank == 0:
-            ncvar = outnf.createVariable(name, vtyp, coord, fill_value=NCFillValue)
-            ncvar.title = title
-            ncvar.units = units
-        else :
-            ncvar = np.zeros((1,1,1))
-        ncvar[:] = part.gather(var)
-    #
-    #
-    #
-    def dumpnetcdf(self, filename, globalgrid, procgrid, part, tcst) :
-        #
-        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('z', self.nbasmax)
-            outnf.createDimension('bnd', nbcorners)
-            outnf.createDimension('inflow', self.max_inflow)
-        else :
-            outnf = None
-        #
-        addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorners, cornerind)
-        addenvironment(outnf, procgrid, part, vtyp, NCFillValue, self.nbpt)
-        addparameters(outnf, part, tcst, vtyp, NCFillValue)
-        #
-        # land grid index -> to facilitate the analyses of the routing
-        #
-        nbpt_loc = np.zeros((self.nbpt,1)).astype(np.int32)
-        nbpt_loc[:,0] = np.arange(1, self.nbpt+1)
-        nbpt_glo = part.l2glandindex(nbpt_loc)
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), "nbpt_glo", "Land point indices on global grid", "-", nbpt_glo[:,0], vtyp)
-        #
-        ################
-        #
-        # Conversion of grid indices for route_togrid
-        grgrid = part.l2glandindex(self.route_togrid[:,:])
-        #
-        #
-        # The field route_togrid is with indices on the local grid. That needs to be converted to the global grid.
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "routetogrid", "Grid into which the basin flows", "-", grgrid, vtyp, "int")
-        # route to basin
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "routetobasin", "Basin in to which the water goes", "-", self.route_tobasin[:,:], vtyp, "int")
-        # basin id
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basinid", "ID of basin", "-", self.global_basinid[:,:], vtyp, "int")
-        #
-        # basin area
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basin_area", "area of basin", "m^2", self.routing_area[:,:], vtyp, "float")
-
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basin_orog_mean", "Mean orography", "m", self.routing_orog_mean[:,:], vtyp, "float")
-
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basin_orog_min", "Min orography", "m", self.routing_orog_min[:,:], vtyp, "float")
-
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basin_orog_max", "Max orography", "m", self.routing_orog_max[:,:], vtyp, "float")
-
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basin_floodp", "Fraction of floodplains", "-", self.routing_floodp[:,:], vtyp, "float")
-
-        # route number into basin
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), "routenbintobas", "Number of basin into current one", "-", self.route_nbintobas[:], vtyp, "int")
-        #
-        # original number into basin
-        self.add_variable(outnf, procgrid, NCFillValue, part, ( 'y','x'), "originnbintobas", "Number of sub-grid basin into current one before truncation", "-", self.origin_nbintobas[:], vtyp, "int")
-        #
-        # latitude of outlet
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "outletlat", "Latitude of Outlet", "degrees north", self.route_outlet[:,:,0], vtyp, "float")
-        # longitude of outlet
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "outletlon", "Longitude of Outlet", "degrees east", self.route_outlet[:,:,1], vtyp, "float")
-        # type of outlet
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "outlettype", "Type of outlet", "code", self.route_type[:,:], vtyp, "float")
-        #
-        # topographic index
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "topoindex", "Topographic index of the retention time", "m", self.topo_resid[:,:], vtyp, "float")
-        #
-
-        # Inflow number
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "route_innum", "Number of inflow", "-", self.route_innum[:,:], vtyp, "int")
-        #
-        # Inflow grid
-        #gingrid = part.l2glandindex(self.inflow_grid[:,:,:inflow_size])
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow', 'z','y','x'), "route_ingrid", "Grid from which the water flows", "-", self.route_ingrid[:,:,:], vtyp, "int")
-        #
-        # Inflow basin
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow', 'z','y','x'), "route_inbasin", "Basin from which the water flows", "-", self.route_inbasin[:,:,:], vtyp, "int")
-
-        #
-        # Save centre of gravity of HTU
-        #
-        # Check if it works
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "CG_lon", "Longitude of centre of gravity of HTU", "degrees east", self.routing_cg[:,:,1], vtyp, "float")
-
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "CG_lat", "Latitude of centre of gravity of HTU", "degrees north", self.routing_cg[:,:,0], vtyp, "float")
-        #
-        # Save the fetch of each basin
-        #
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "fetch", "Fetch contributing to each HTU", "m^2", self.routing_fetch[:,:], vtyp, "float")
-        #