Interface.py 28.3 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
#
#
import numpy as np
import os
import pickle
from netCDF4 import Dataset
import RPPtools as RPP
from mpi4py import MPI
#
import sys
sys.path.append(os.getcwd()+'/F90subroutines')
if MPI.COMM_WORLD.Get_rank() == 0 :
    err=os.system("cd F90subroutines; 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
24
config = configparser.ConfigParser({'Documentation':'false', 'nbxmax':'63', 'ROUTING_RIVERS':'50'})
25
config.read("run.def")
26 27 28
gendoc = config.get("OverAll", "Documentation")
nbxmax = config.getint("OverAll", "nbxmax")
largest_pos = config.getint("OverAll", "ROUTING_RIVERS")
29 30
#
undef_int = 999999999.9
31 32
# Order of magnitude for the area precision in m^2.
prec = 100.0
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57
#
# 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.
#
58
def initatmgrid(rank, nbcore, nbpt, modelgrid) :
59 60 61
    print("INITATMGRID corners", np.array(modelgrid.polyll).shape)
    print("INITATMGRID area", np.array(modelgrid.area).shape)
    print("INITATMGRID neighbours", np.array(modelgrid.neighbours).shape)
62 63 64 65 66 67 68 69
    routing_interface.initatmgrid(rank, nbcore, modelgrid.polyll, modelgrid.coordll, modelgrid.area, modelgrid.contfrac, modelgrid.neighbours)
    return
#
#
#
def closeinterface(comm) :
    comm.Barrier()
    routing_interface.closeinterface()
70 71 72
    return
#
#
73
#
74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173
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
#
#
#
174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214
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
#
#
#
215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265
class HydroOverlap :
#
    def __init__(self, nbpt, nbvmax, sub_pts, sub_index_in, sub_area_in, sub_lon_in, sub_lat_in,  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]]
        #
        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
        #
266 267 268 269
        return
#
#
#
270 271
class HydroSuper :
    def __init__(self, nbvmax, hydrodata, hydrooverlap) :
272 273 274 275
        #
        # Keep largest possible number of HTUs
        #
        self.nbhtuext = nbvmax
276 277 278 279 280 281 282 283 284 285 286
        #
        # 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
        #
287 288 289 290 291 292
        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, \
293
            self.basin_lshead, self.outflow_grid, self.outflow_basin, self.nbcoastal, self.coastal_basin = \
294 295 296
                    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, \
297
                                                basin_bbout, basin_lshead, coast_pts, hydrooverlap.nwbas)
298
        self.nbpt = self.basin_count.shape[0]
299 300
        return
    #
301 302 303 304 305 306 307 308 309
    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))
310 311
        return
    #
POLCHER Jan's avatar
POLCHER Jan committed
312
    def fetch(self, part) :
313
        #
314
        fetch_basin = np.zeros(self.basin_area.shape, dtype=np.float32, order='F')
315
        #
316
        self.basin_area = routing_interface.areanorm(self.basin_count, self.basin_area, self.outflow_grid)
317
        partial_sum = np.zeros(self.basin_area.shape, dtype=np.float32, order='F')
318
        #
319
        old_sorted = np.zeros(largest_pos, dtype=np.float32, order='F')
320 321 322 323 324
        #
        maxdiff_sorted = prec*prec
        iter_count = 0
        #
        while iter_count < part.size*3 and maxdiff_sorted > prec :
325
            fetch_basin[:,:] = 0.0
326
            outflow_uparea = routing_interface.fetch(part.landcorelist, self.basin_count, self.basin_area, self.basin_id, self.basin_hierarchy, \
327
                                                         self.basin_fac, self.outflow_grid, self.outflow_basin, partial_sum, fetch_basin)
328 329
            partial_sum = np.copy(fetch_basin)
            part.landsendtohalo(partial_sum, order='F')
330 331
            partial_sum = part.zerocore(partial_sum, order='F')
            #
332 333 334 335 336 337 338 339 340 341 342 343
            # 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
344

345 346
        self.fetch_basin = np.copy(fetch_basin)
        #
347
        # Upstream area of the smalest river we call largest rivers. 
348
        #
349
        self.largest_rivarea = sorted_outareas[largest_pos-1]
350 351 352
        #
        #
        #
353 354 355
        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)
356
        return
357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456
    #
    #
    #
    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
457 458 459
#
#
#
460
class HydroGraph :
POLCHER Jan's avatar
POLCHER Jan committed
461
    def __init__(self, nbasmax, hydrosuper, part) :
462
        #
463
        self.nbasmax = nbasmax
464 465
        self.nbpt = hydrosuper.basin_count.shape[0]
        #
466 467 468 469
        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, \
470 471 472 473
                                                               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)
474 475 476 477 478 479
        #
        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)
        #
480 481
        return
    #
482 483 484 485 486
    def dumpnetcdf(self, filename, globalgrid, procgrid, part) :
        #
        NCFillValue=1.0e20
        vtyp=np.float64
        cornerind=[0,2,4,6]
487
        nbcorners = len(cornerind)
488 489 490 491 492 493 494 495
        #
        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)
496
            outnf.createDimension('bnd', nbcorners)
497
        else :
498
            outnf = None
499
        #
500 501
        addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorners, cornerind)
        addenvironment(outnf, procgrid, part, vtyp, NCFillValue, self.nbpt)
502
        #
503 504
        # The field route_togrid is with indices on the local grid. That needs to be converted to the global grid.
        #
505 506 507 508 509 510 511 512 513
        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 :
POLCHER Jan's avatar
POLCHER Jan committed
514
            print(part.rank, itarget, " Before route_togrid = ", self.route_togrid[itarget,:])
515
        grgrid = part.l2glandindex(self.route_togrid[:,:])
516
        if itarget >+ 0 :
POLCHER Jan's avatar
POLCHER Jan committed
517
            print(part.rank, itarget, " After route_togrid = ", self.route_togrid[itarget,:])
518
        rgrid = procgrid.landscatter(grgrid.astype(vtyp), order='F')
519 520 521 522 523 524 525 526 527
        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)
        #
528
        rtobasin = procgrid.landscatter(self.route_tobasin[:,:].astype(vtyp), order='F')
529
        rtobasin = rtobasin.astype(vtyp, copy=False)
530 531 532 533 534 535 536 537 538
        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)
        #
539
        rid = procgrid.landscatter(self.global_basinid[:,:].astype(vtyp), order='F')
540 541 542 543 544 545 546 547 548
        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)
        #
549
        rintobas = procgrid.landscatter(self.route_nbintobas[:].astype(vtyp))
550 551 552 553 554 555 556 557 558
        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)
        #
559
        onbintobas = procgrid.landscatter(self.origin_nbintobas[:].astype(vtyp))
560 561 562 563 564 565 566 567 568
        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)
        #
569
        olat = procgrid.landscatter(self.route_outlet[:,:,0].astype(vtyp), order='F')
570 571 572 573 574 575 576 577 578
        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)
        #
579
        olon = procgrid.landscatter(self.route_outlet[:,:,1].astype(vtyp), order='F')
580 581 582 583 584 585 586 587 588
        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)
        #
589
        otype = procgrid.landscatter(self.route_type[:,:].astype(vtyp), order='F')
590 591 592 593 594 595 596 597 598
        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)
        #
599
        tind = procgrid.landscatter(self.topo_resid[:,:].astype(vtyp), order='F')
600 601 602 603 604 605 606 607 608
        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)
        #
609 610
        # Save centre of gravity of HTU
        #
611
        cg = procgrid.landscatter(self.routing_cg[:,:,:].astype(vtyp), order='F')
612 613 614 615 616 617 618 619 620 621 622 623 624 625
        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,:,:,:])
        #
626 627
        # Save the fetch of each basin
        #
628
        fe =  procgrid.landscatter(self.routing_fetch[:,:].astype(vtyp), order='F')
629 630 631 632 633 634 635 636 637
        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)
        #
638 639
        if part.rank == 0 :
            outnf.close()
640 641
        #
        return
642 643 644 645