# # 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.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 # # # 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(np.array(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 = \ 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. # 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) : # # Keep largest possible number of HTUs # self.nbasmax = nbasmax self.nbhtuext = nbvmax self.nwbas = hydrooverlap.nwbas # # 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)) 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 # 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 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 killbas(self, tokill, totakeover, numops): ops = tokill.shape[1] routing_interface.killbas(nbpt = self.nbpt, nbxmax_in = self.nbxmax_in, 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_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.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, modelgrid) : # self.nbasmax = nbasmax self.nbpt = hydrosuper.basin_count.shape[0] nwbas = hydrosuper.basin_topoind.shape[1] nbxmax_in = hydrosuper.inflow_grid.shape[1] # 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.finish_truncate(nbpt = self.nbpt, nbxmax_in = nbxmax_in, nbasmax = nbasmax, nwbas = nwbas, num_largest = hydrosuper.num_largest, gridarea = modelgrid.area, cfrac = modelgrid.contfrac, basin_count = hydrosuper.basin_count, \ basin_notrun = hydrosuper.basin_notrun, basin_area = hydrosuper.basin_area, 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.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, nbxmax_in = nbxmax_in, 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[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) : # 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) 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) # # 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", "Grid point Global", "-", nbpt_glo[:,0], vtyp) # ################ # # TEST: l2glandindex 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,:]) # Conversion grgrid = part.l2glandindex(self.route_togrid[:,:]) if itarget >+ 0 : print(part.rank, itarget, " After route_togrid = ", self.route_togrid[itarget,:]) ################ # # 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, ('htu', 'y','x'), "routetogrid", "Grid into which the basin flows", "-", grgrid, vtyp, "int") # route to basin self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "routetobasin", "Basin in to which the water goes", "-", self.route_tobasin[:,:], vtyp, "int") # basin id self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "basinid", "ID of basin", "-", self.global_basinid[:,:], vtyp, "int") # # basin area self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "basin_area", "area of basin", "m^2", self.routing_area[:,:], 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, ('htu','y','x'), "outletlat", "Latitude of Outlet", "degrees north", self.route_outlet[:,:,0], vtyp, "float") # longitude of outlet self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "outletlon", "Longitude of Outlet", "degrees east", self.route_outlet[:,:,1], vtyp, "float") # type of outlet self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "outlettype", "Type of outlet", "code", self.route_type[:,:], vtyp, "float") # # topographic index self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "topoindex", "Topographic index of the retention time", "m", self.topo_resid[:,:], vtyp, "float") # # Inflow number self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','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', 'htu','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', 'htu','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, ('htu','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, ('htu','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, ('htu','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