#
#
import numpy as np
import os
import pickle
from netCDF4 import Dataset
import RPPtools as RPP
from mpi4py import MPI
#
import sys
from inspect import currentframe, getframeinfo
#
localdir=os.path.dirname(getframeinfo(currentframe()).filename)
sys.path.append(localdir+'/F90subroutines')
F90=localdir+'/F90subroutines'
if MPI.COMM_WORLD.Get_rank() == 0 :
    err=os.system("cd "+F90+"; make all")
    if err != 0 :
        print("Compilation error in the FORTRAN interface")
        sys.exit()
else :
    print("Not compiling on other cores")
MPI.COMM_WORLD.Barrier()
#
import routing_interface
#
import configparser
config = configparser.ConfigParser({'Documentation':'false', 'nbxmax':'63', 'ROUTING_RIVERS':'50'})
config.read("run.def")
gendoc = config.get("OverAll", "Documentation")
nbxmax = config.getint("OverAll", "nbxmax")
largest_pos = config.getint("OverAll", "ROUTING_RIVERS")
#
undef_int = 999999999.9
# Order of magnitude for the area precision in m^2.
prec = 100.0
#
# Print the documentation for the FORTRAN interface
#
if gendoc.lower() == "true" : 
    docwrapper = open('DocumentationInterface', 'w')
    docwrapper.write(routing_interface.initatmgrid.__doc__)
    docwrapper.write("====================================================================\n")
    docwrapper.write(routing_interface.gethydrogrid.__doc__)
    docwrapper.write("====================================================================\n")
    docwrapper.write(routing_interface.findbasins.__doc__)
    docwrapper.write("====================================================================\n")
    docwrapper.write(routing_interface.globalize.__doc__)
    docwrapper.write("====================================================================\n")
    docwrapper.write(routing_interface.linkup.__doc__)
    docwrapper.write("====================================================================\n")
    docwrapper.write(routing_interface.fetch.__doc__)
    docwrapper.write("====================================================================\n")
    docwrapper.write(routing_interface.truncate.__doc__)
    docwrapper.close
#
# Functions to access the interfaces
#
#
# initatmgrid : Initialises the grid.f90 module and passes the description of the atmospheric grid.
#
def initatmgrid(rank, nbcore, nbpt, modelgrid) :
    print("INITATMGRID corners", np.array(modelgrid.polyll).shape)
    print("INITATMGRID area", np.array(modelgrid.area).shape)
    print("INITATMGRID neighbours", np.array(modelgrid.neighbours).shape)
    routing_interface.initatmgrid(rank, nbcore, modelgrid.polyll, modelgrid.coordll, modelgrid.area, modelgrid.contfrac, modelgrid.neighbours)
    return
#
#
#
def closeinterface(comm) :
    comm.Barrier()
    routing_interface.closeinterface()
    return
#
#
#
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')
    partial_sum = np.zeros(routing_area.shape, dtype=np.float32, order='F')
    old_sorted = np.zeros(largest_pos, dtype=np.float32, order='F')
    #
    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(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 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.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) :
        #
        # Reshape stuff so that it fits into arrays
        #
        sub_index = np.zeros((nbpt,nbvmax,2), dtype=np.int8, order='F')
        sub_area = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
        sub_lon = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
        sub_lat = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
        for ib in range(nbpt) :
            sub_area[ib,0:sub_pts[ib]] = sub_area_in[ib][:]
            sub_lon[ib,0:sub_pts[ib]] = sub_lon_in[ib][:]
            sub_lat[ib,0:sub_pts[ib]] = sub_lat_in[ib][:]
            for ip in range(sub_pts[ib]) :
                sub_index[ib,ip,:] = [sub_index_in[ib][0][ip],sub_index_in[ib][1][ip]]
        #
        part.landsendtohalo(sub_area, order='F')
        #
        trip_tmp = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
        basins_tmp = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
        topoind_tmp = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
        fac_tmp = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
        hierarchy_tmp = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
        #
        trip_tmp[:,:] = np.nan
        basins_tmp[:,:] = np.nan
        for ib in range(nbpt) :
            trip_tmp[ib,0:sub_pts[ib]] = np.asarray(hydrodata.trip[ib][:])
            basins_tmp[ib,0:sub_pts[ib]] = np.asarray(hydrodata.basins[ib][:])
            topoind_tmp[ib,0:sub_pts[ib]] = np.asarray(hydrodata.topoind[ib][:])
            fac_tmp[ib,0:sub_pts[ib]] = np.asarray(hydrodata.fac[ib][:])
            hierarchy_tmp[ib,0:sub_pts[ib]] = np.asarray(hydrodata.disto[ib][:])
        #
        trip_tmp[np.isnan(trip_tmp)] = undef_int
        basins_tmp[np.isnan(trip_tmp)] = undef_int
        #
        # Go to the call of the FORTRAN interface
        #
        print("GETHYDROGRID : nbpt = ", nbpt, nbvmax)
        print("GETHYDROGRID : nbvmax = ", nbvmax)
        print("GETHYDROGRID : nbxmax = ", nbxmax)
        self.nbi, self.nbj, self.area_bx, self.trip_bx, self.basin_bx, self.topoind_bx, self.fac_bx, self.hierarchy_bx, \
            self.lon_bx, self.lat_bx, self.lshead_bx, self.nwbas = \
                    routing_interface.gethydrogrid(nbxmax, sub_pts, sub_index, sub_area, \
                    hydrodata.basinsmax, hydrodata.topoindmin, sub_lon, sub_lat, trip_tmp, basins_tmp, topoind_tmp, fac_tmp, hierarchy_tmp)
        #
        # Plot some diagnostics for the hydrology grid within the atmospheric meshes.
        #
        # Clean-up these arrays so that they are easy to use in Python.
        self.lon_bx[self.lon_bx > 360.]=np.nan
        self.lat_bx[self.lat_bx > 90.]=np.nan
        #
        return
#
#
#
class HydroSuper :
    def __init__(self, nbvmax, hydrodata, hydrooverlap) :
        #
        # Keep largest possible number of HTUs
        #
        self.nbhtuext = nbvmax
        #
        # Call findbasins
        #
        nb_basin, basin_inbxid, basin_outlet, basin_outtp, self.basin_sz, basin_bxout, basin_bbout, self.basin_pts, basin_lshead, coast_pts = \
                    routing_interface.findbasins(nbvmax, hydrooverlap.nbi, hydrooverlap.nbj, hydrooverlap.trip_bx, \
                                                 hydrooverlap.basin_bx, hydrooverlap.fac_bx, hydrooverlap.hierarchy_bx, \
                                                 hydrooverlap.topoind_bx, hydrooverlap.lshead_bx, \
                                                 hydrooverlap.lon_bx, hydrooverlap.lat_bx)
        #
        # Call Globalize
        #
        lon_bx_tmp = hydrooverlap.lon_bx
        lon_bx_tmp[np.isnan(lon_bx_tmp)] = undef_int
        lat_bx_tmp = hydrooverlap.lat_bx
        lat_bx_tmp[np.isnan(lat_bx_tmp)] = undef_int
        self.basin_count, self.basin_notrun, self.basin_area, self.basin_cg, self.basin_hierarchy, self.basin_fac, self.basin_topoind, \
            self.basin_id, self.basin_outcoor, self.basin_type, self.basin_flowdir, \
            self.basin_lshead, self.outflow_grid, self.outflow_basin, self.nbcoastal, self.coastal_basin = \
                    routing_interface.globalize(hydrooverlap.area_bx, lon_bx_tmp, lat_bx_tmp, hydrooverlap.trip_bx, \
                                                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) :
        #
        # Call the linkup routine in routing_reg.
        #
        print("Invented basins =", hydrodata.basinsmax)
        self.inflow_number,self.inflow_grid,self.inflow_basin = routing_interface.linkup(nbxmax, self.basin_count, self.basin_area, self.basin_id, \
                                                                       self.basin_flowdir, self.basin_lshead, self.basin_hierarchy, \
                                                                       self.basin_fac, self.outflow_grid, self.outflow_basin, \
                                                                       self.nbcoastal, self.coastal_basin, float(hydrodata.basinsmax))
        return
    #
    def fetch(self, part) :
        #
        fetch_basin = np.zeros(self.basin_area.shape, dtype=np.float32, order='F')
        #
        self.basin_area = routing_interface.areanorm(self.basin_count, self.basin_area, self.outflow_grid)
        partial_sum = np.zeros(self.basin_area.shape, dtype=np.float32, order='F')
        #
        old_sorted = np.zeros(largest_pos, dtype=np.float32, order='F')
        #
        maxdiff_sorted = prec*prec
        iter_count = 0
        #
        while iter_count < part.size*3 and maxdiff_sorted > prec :
            fetch_basin[:,:] = 0.0
            outflow_uparea = routing_interface.fetch(part.landcorelist, self.basin_count, self.basin_area, self.basin_id, self.basin_hierarchy, \
                                                         self.basin_fac, self.outflow_grid, self.outflow_basin, partial_sum, fetch_basin)
            partial_sum = np.copy(fetch_basin)
            part.landsendtohalo(partial_sum, order='F')
            partial_sum = part.zerocore(partial_sum, order='F')
            #
            # Find area the largest basins need at least to have.
            #
            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 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

        self.fetch_basin = np.copy(fetch_basin)
        #
        # Upstream area of the smalest river we call largest rivers. 
        #
        self.largest_rivarea = sorted_outareas[largest_pos-1]
        #
        #
        #
        self.num_largest = routing_interface.rivclassification(part.landcorelist, self.basin_count, self.outflow_grid, self.outflow_basin, \
                                                               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 add_variable(self,outnf, procgrid, NCFillValue, part, coord, name, title, units, data, vtyp):
        var = procgrid.landscatter(data.astype(vtyp), order='F')
        var[np.isnan(var)] = 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) :
        #
        NCFillValue=1.0e20
        vtyp=np.float64
        inflow_size = 100
        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('htu', self.inflow_number.shape[1])
            outnf.createDimension('in',inflow_size )
            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 !
        #
        # nbpt_glo
        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", "Grid point Global", "-", nbpt_glo[:,0], vtyp)
        #
        # contfrac
        contfrac = np.array(procgrid.contfrac)
        self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), "contfrac", "Land fraction", "-", np.array(procgrid.contfrac), vtyp)
        #
        # basin_id
        self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "basin_id", "ID for each HTU", "-", self.basin_id, vtyp)
        #
        #self.basin_count
        self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), "basin_count", "HTU count", "-", self.basin_count, vtyp)
        #
        # self.basin_notrun
        self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), "basin_notrun", "Not run", "-", self.basin_notrun, vtyp)
        #
        # self.basin_area
        self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "basin_area", "Basin area", "-", self.basin_area, vtyp)
        #
        # self.basin_cg
        self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "CG_lon", "CG lon", "-", self.basin_cg[:,:,1], vtyp)
        self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "CG_lat", "CG lat", "-", self.basin_cg[:,:,0], vtyp)
        #
        # self.topoind
        self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "basin_topoind", "Topoindex", "-", self.basin_topoind, vtyp)
        # 
        # outcoor
        self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "outcoor_lon", "outcoor lon", "-", self.basin_outcoor[:,:,1], vtyp)
        self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "outcoor_lat", "outcoor lat", "-", self.basin_outcoor[:,:,0], vtyp)
        # 
        # type
        self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "basin_type", "type", "-", self.basin_type, vtyp)
        #
        # flowdir
        self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "basin_flowdir", "flowdir", "-", self.basin_flowdir, vtyp)
        #
        # 
        #self.outflow_grid
        grgrid = part.l2glandindex(self.outflow_grid)
        grgrid[self.outflow_grid == 0 ] = -2 # in case it flows out of the domain, the 0 should not remain
        grgrid[self.outflow_grid == -1 ] = -1
        grgrid[self.outflow_grid == -2 ] = -2
        grgrid[self.outflow_grid == -3 ] = -3
        self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "HTUoutgrid", "HTU outflow grid", "-", grgrid, vtyp)
        #
        #self.outflow_basin
        self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "HTUoutbasin", "Outflow HTU of grid", "-", self.outflow_basin, vtyp)
        #
        # self.inflow_number
        self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "HTUinnum", "Inflow number", "-", self.inflow_number, vtyp)
        #
        # Inflow Grid -> convert to global
        gingrid = part.l2glandindex(self.inflow_grid[:,:,:inflow_size])
        self.add_variable(outnf, procgrid, NCFillValue, part, ('in','htu','y','x'), "HTUingrid", "Inflow grid", "-", gingrid, vtyp)
        #
        # Inflow Basin
        self.add_variable(outnf, procgrid, NCFillValue, part, ('in','htu','y','x'), "HTUinbas", "Inflow basin", "-", self.inflow_basin[:,:,:inflow_size], vtyp)
        #
        # Save the fetch of each basin
        #
        self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "fetch_basin", "Fetch contributing to each HTU", "m^2", self.fetch_basin, vtyp)
        #
        # 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, \
                                                               hydrosuper.basin_notrun, hydrosuper.basin_area, hydrosuper.basin_cg, \
                                                               hydrosuper.basin_topoind, hydrosuper.fetch_basin, hydrosuper.basin_id, \
                                                               hydrosuper.basin_outcoor, hydrosuper.basin_type, hydrosuper.basin_flowdir, \
                                                               hydrosuper.outflow_grid, hydrosuper.outflow_basin, \
                                                               hydrosuper.inflow_number,hydrosuper.inflow_grid,hydrosuper.inflow_basin)
        #
        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)
        #
        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('htu', self.nbasmax)
            outnf.createDimension('bnd', nbcorners)
        else :
            outnf = None
        #
        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.
        #
        itarget=-1
        for il in range(procgrid.nbland) :
            lo = procgrid.lon_full[procgrid.indP[il][0],procgrid.indP[il][1]]
            la = procgrid.lat_full[procgrid.indP[il][0],procgrid.indP[il][1]]
            d=np.sqrt((lo-3.13)**2+(la-39.70)**2)
            if d < 0.05 :
                itarget = il
                
        if itarget >+ 0 :
            print(part.rank, itarget, " Before route_togrid = ", self.route_togrid[itarget,:])
        grgrid = part.l2glandindex(self.route_togrid[:,:])
        if itarget >+ 0 :
            print(part.rank, itarget, " After route_togrid = ", self.route_togrid[itarget,:])
        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)
            routetogrid.title = "Grid into which the basin flows"
            routetogrid.units = "-"
        else :
            routetogrid = np.zeros((1,1,1))    
        routetogrid[:,:,:] = part.gather(rgrid)
        #
        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 :
            routetobasin = outnf.createVariable("routetobasin", vtyp, ('htu','y','x'), fill_value=NCFillValue)
            routetobasin.title = "Basin in to which the water goes"
            routetobasin.units = "-"
        else :
            routetobasin = np.zeros((1,1,1))
        routetobasin[:,:,:] = part.gather(rtobasin)
        #
        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)
            basinid.title = "ID of basin"
            basinid.units = "-"
        else :
            basinid = np.zeros((1,1,1))
        basinid[:,:,:] = part.gather(rid)
        #
        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)
            routenbintobas.title = "Number of basin into current one"
            routenbintobas.units = "-"
        else :
            routenbintobas = np.zeros((1,1))
        routenbintobas[:,:] = part.gather(rintobas)
        #
        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)
            originnbintobas.title = "Number of sub-grid basin into current one before truncation"
            originnbintobas.units = "-"
        else :
            originnbintobas = np.zeros((1,1))
        originnbintobas[:,:] = part.gather(onbintobas)
        #
        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)
            outletlat.title = "Latitude of Outlet"
            outletlat.title = "degrees north"
        else :
            outletlat = np.zeros((1,1,1))
        outletlat[:,:,:] = part.gather(olat)
        #
        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)
            outletlon.title = "Longitude of outlet"
            outletlon.units = "degrees east"
        else :
            outletlon = np.zeros((1,1,1))
        outletlon[:,:,:] = part.gather(olon)
        #
        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)
            outlettype.title = "Type of outlet"
            outlettype.units = "code"
        else :
            outlettype = np.zeros((1,1,1))
        outlettype[:,:,:] = part.gather(otype)
        #
        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)
            topoindex.title = "Topographic index of the retention time"
            topoindex.units = "m"
        else :
            topoindex = np.zeros((1,1,1))
        topoindex[:,:,:] = part.gather(tind)
        #
        # Save centre of gravity of HTU
        #
        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)
            CG_lon.title = "Longitude of centre of gravity of HTU"
            CG_lon.units = "degrees east"
            CG_lat = outnf.createVariable("CG_lat", vtyp, ('htu','y','x'), fill_value=NCFillValue)
            CG_lat.title = "Latitude of centre of gravity of HTU"
            CG_lat.units = "degrees north"
        else :
            CG_lon = np.zeros((1,1,1))
            CG_lat = np.zeros((1,1,1))
        CG_lon[:,:,:] = part.gather(cg[1,:,:,:])
        CG_lat[:,:,:] = part.gather(cg[0,:,:,:])
        #
        # Save the fetch of each basin
        #
        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)
            fetch.title = "Fetch contributing to each HTU"
            fetch.units = "m^2"
        else :
            fetch = np.zeros((1,1,1))
        fetch[:,:,:] = part.gather(fe)
        #
        if part.rank == 0 :
            outnf.close()
        #
        return