Interface.py

#
#
import numpy as np
import os
import pickle
from netCDF4 import Dataset
from mpi4py import MPI
import gc
#
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 NcHelp as NH
#
import getargs
config = getargs.SetupConfig()
gendoc = config.get("OverAll", "Documentation", fallback='false')
nbxmax = config.getint("OverAll", "nbxmax", fallback=63)
largest_pos = config.getint("OverAll", "ROUTING_RIVERS", fallback=200)
#
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.finish_truncate.__doc__)
    docwrapper.write("====================================================================\n")
    docwrapper.write(routing_interface.killbas.__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
#
#
#
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.int32, 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][:,ip]
        #
        part.landsendtohalo(np.array(sub_area), order='F')
        #
        ijdim=[]
        for ib in range(nbpt) :
            ijdim.append(max(np.max(sub_index[ib,:sub_pts[ib],0])-np.min(sub_index[ib,:sub_pts[ib],0])+1,np.max(sub_index[ib,:sub_pts[ib],1])-np.min(sub_index[ib,:sub_pts[ib],1])+1))
        ijdimmax = max(ijdim)
        #
        print("GETHYDROGRID : nbpt = {0}".format(nbpt))
        print("GETHYDROGRID : nbvmax = {0}".format(nbvmax))
        print("GETHYDROGRID : ijdimmax = {0}".format(ijdimmax))
        #
        del sub_area_in; del sub_lon_in; del sub_lat_in; del sub_index_in
        #
        #
        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')
        orog_tmp = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
        floodp_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][:])
            orog_tmp[ib,0:sub_pts[ib]] = np.asarray(hydrodata.orog[ib][:])
            floodp_tmp[ib,0:sub_pts[ib]] = np.asarray(hydrodata.floodplains[ib][:])
        #
        del hydrodata.trip; del hydrodata.basins; del hydrodata.topoind
        del hydrodata.fac; del hydrodata.disto; del hydrodata.orog;
        del hydrodata.floodplains
        #
        trip_tmp[np.isnan(trip_tmp)] = undef_int
        basins_tmp[np.isnan(trip_tmp)] = undef_int
        #
        # Go to the call of the FORTRAN interface
        #
        self.nbi, self.nbj, self.area_bx, self.trip_bx, self.basin_bx, self.topoind_bx, self.fac_bx, self.hierarchy_bx, \
            self.orog_bx, self.floodp_bx, \
            self.lon_bx, self.lat_bx, self.lshead_bx = \
                    routing_interface.gethydrogrid(ijdimmax, sub_pts, sub_index, sub_area, \
                    hydrodata.basinsmax, hydrodata.topoindmin, sub_lon, sub_lat, trip_tmp, basins_tmp, topoind_tmp, fac_tmp,\
                        hierarchy_tmp, orog_tmp, floodp_tmp)
        #
        del trip_tmp; del basins_tmp; del topoind_tmp
        del fac_tmp; del hierarchy_tmp; del orog_tmp
        del floodp_tmp
        #
        self.nwbas = nbvmax
        # 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, nbasmax, part) :
        #
        # Keep largest possible number of HTUs
        #
        self.nbasmax = nbasmax
        self.nbhtuext = nbvmax
        self.nbpt = hydrooverlap.nbi.shape[0]
        #
        # nb_htu can be adjusted with self.nwbas
        # nb_htu can be lowered with a larger maxpercent (routing_reg.f90)
        nb_htu = nbvmax
        nbv = 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(nbpt = self.nbpt, nb_htu = self.nbhtuext, nbv = nbv, nbi = hydrooverlap.nbi, \
                                                 nbj = hydrooverlap.nbj, trip_bx = hydrooverlap.trip_bx, \
                                                 basin_bx = hydrooverlap.basin_bx, fac_bx = hydrooverlap.fac_bx, \
                                                 hierarchy_bx = hydrooverlap.hierarchy_bx, \
                                                 topoind_bx = hydrooverlap.topoind_bx, lshead_bx = hydrooverlap.lshead_bx, \
                                                 lontmp = hydrooverlap.lon_bx, lattmp = hydrooverlap.lat_bx)
        #
        # Adjust nwbas to the maximum found over the domain
        #
        self.nwbas = part.domainmax(np.max(nb_basin))
        # Set the number of inflows per basin. For the moment twice the maximum number of basins.
        self.inflowmax = max(10, self.nwbas*2)
        print("Maximum number of basin created : {0}".format(self.nwbas))
        ijdim=[]
        for i in range(self.nbpt) :
            ijdim.append(max(hydrooverlap.area_bx[i,:,:].shape))
        self.ijdimmax = max(ijdim)
        #
        # 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_orog_mean, self.basin_orog_min, self.basin_orog_max, self.basin_floodp, \
            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(nbpt = self.nbpt, nb_htu = self.nbhtuext, nbv = nbv, ijdimmax = self.ijdimmax, \
                                                area_bx = hydrooverlap.area_bx, lon_bx = lon_bx_tmp, lat_bx = lat_bx_tmp, trip_bx = hydrooverlap.trip_bx, \
                                                hierarchy_bx = hydrooverlap.hierarchy_bx, orog_bx = hydrooverlap.orog_bx, floodp_bx =  hydrooverlap.floodp_bx,\
                                                fac_bx = hydrooverlap.fac_bx, topoind_bx = hydrooverlap.topoind_bx, min_topoind = hydrodata.topoindmin, \
                                                nb_basin = nb_basin, basin_inbxid = basin_inbxid, basin_outlet = basin_outlet, basin_outtp = basin_outtp, \
                                                basin_sz = self.basin_sz, basin_pts = self.basin_pts, basin_bxout = basin_bxout, \
                                                basin_bbout = basin_bbout, lshead = basin_lshead, coast_pts = coast_pts, nwbas = self.nwbas)

        # Memory management
        del basin_bbout; del basin_lshead; del coast_pts; del basin_bxout; del self.basin_pts;
        del basin_outtp; del basin_outlet; del basin_inbxid; del nb_basin
        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(self.ijdimmax, self.inflowmax, 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))
        self.nbxmax_in = self.inflow_number.shape[1]
        #
        #
        #
        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
        #
        fhalolist = [i+1 for i in range(self.nbpt) if i not in part.landcorelist]
        #
        while iter_count < part.size*3 and maxdiff_sorted > prec :
            fetch_basin[:,:] = 0.0
            outflow_uparea = routing_interface.fetch(fhalolist, 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.
            l = min(sorted_outareas.shape[0],largest_pos)
            if part.size == 1 :
                maxdiff_sorted = 0.0
            else :
                maxdiff_sorted = np.max(np.abs(sorted_outareas[0:l]-old_sorted[0:l]))
                old_sorted[:l] = 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[l-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 check_fetch(self):

        routing_interface.checkfetch(nbpt = self.nbpt, nwbas = self.nwbas, fetch_basin = self.fetch_basin, outflow_grid = self.outflow_grid, \
                                     outflow_basin = self.outflow_basin, basin_count = self.basin_count)

        return

    def check_routing(self):

        routing_interface.checkrouting(nbpt = self.nbpt, nwbas = self.nwbas, outflow_grid = self.outflow_grid, outflow_basin = self.outflow_basin, \
                                       basin_count = self.basin_count)

        return

    def correct_outflows(self, part):
        # Global index of the proc domain
        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)
        # Halo points
        fhalo = np.array([pt+1 for pt in range(self.nbpt) if pt not in part.landcorelist], order = "F")
        #
        # Outflow grid in global index and send to halo
        hg = np.copy(self.outflow_grid)
        hg = part.l2glandindex(self.outflow_grid)
        part.landsendtohalo(hg, order='F')
        # Convert to local index
        outflows = np.unique(hg)
        outflows_out = [a for a in outflows if (a not in nbpt_glo and a>0)]
        # Work in a copy to avoid error
        hg1 = np.copy(hg)
        for a in outflows_out:
          hg1[hg == a] = 0
        for loc, glo in zip(nbpt_loc,nbpt_glo):
          hg1[hg == glo[0]] = loc[0]
        hg = hg1
        del hg1
        # Send Outflow basin to the halo and adapt it
        hb = np.copy(self.outflow_basin)
        part.landsendtohalo(hb, order='F')
        hb[hg <= 0] = 999999999
        for ig in range(self.nbpt):
            hb[ig,self.basin_count[ig]:] = 0
        #
        # Correct the routing graph in the halo
        routing_interface.correct_outflows(nbpt = self.nbpt, nwbas = self.nwbas, nbhalo = fhalo.shape[0], \
                    outflow_grid = self.outflow_grid, outflow_basin = self.outflow_basin, \
                    basin_count = self.basin_count, hg = hg, hb = hb, halopts = fhalo)
        #
        # Correct the inflows
        nbxmax_tmp = self.inflow_grid.shape[2]
        routing_interface.correct_inflows(nbpt = self.nbpt, nwbas = self.nwbas, inflowmax = nbxmax_tmp, outflow_grid = self.outflow_grid, outflow_basin = self.outflow_basin, basin_count = self.basin_count, inflow_number = self.inflow_number, inflow_grid = self.inflow_grid, inflow_basin = self.inflow_basin)

        return


    def killbas(self, tokill, totakeover, numops):
        ops = tokill.shape[1]
        #
        nbxmax_tmp = self.inflow_grid.shape[2]
        #
        routing_interface.killbas(nbpt = self.nbpt, inflowmax = nbxmax_tmp, \
                nbasmax = self.nbasmax, nwbas = self.nwbas, ops = ops, tokill = tokill,\
                totakeover = totakeover, numops = numops, basin_count = self.basin_count,\
                basin_area = self.basin_area, basin_orog_mean = self.basin_orog_mean, \
                basin_orog_min = self.basin_orog_min, basin_orog_max = self.basin_orog_max, basin_floodp = self.basin_floodp, \
                basin_cg = self.basin_cg, basin_topoind = self.basin_topoind, fetch_basin = self.fetch_basin,\
                basin_id = self.basin_id, basin_coor = self.basin_outcoor, basin_type = self.basin_type,\
                basin_flowdir = self.basin_flowdir, outflow_grid = self.outflow_grid, outflow_basin = self.outflow_basin, \
                inflow_number = self.inflow_number, inflow_grid = self.inflow_grid, inflow_basin = self.inflow_basin)

    #
    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.basin_id.shape[1])
            outnf.createDimension('htu', self.inflow_number.shape[1])
            outnf.createDimension('in',inflow_size )
            outnf.createDimension('bnd', nbcorners)
        else :
            outnf = None
        #
        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 !
        #
        # 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_orog
        self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "basin_orog", "Basin orography", "-", self.basin_orog, vtyp)
        #
        # self.basin_floodp
        self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "basin_floodp", "Basin floodplains", "-", self.basin_floodp, 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
#
#
#