# # 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) maxpercent = config.getint("OverAll", "maxpercent", fallback=2) # 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, maxpercent, 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, comm) : # # 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.basin_beta_fp, 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 # # Correct all variables over the halo before the Linkup # except : outflow_grid (corrected after linkup), basin_cg and basin_outcoor (not needed and not valid with part.landsendtohalo) for variable in [self.basin_count, self.basin_notrun, self.basin_area, 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_type,self.basin_flowdir,\ self.basin_lshead, self.outflow_basin,self.nbcoastal,self.coastal_basin]: part.landsendtohalo(variable, order='F') comm.Barrier() # # 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 ind in range(self.nbpt): hg1[hg == nbpt_glo[ind,0]] = nbpt_loc[ind,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, basin_beta_fp = self.basin_beta_fp, 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 get_floodcri(self): nbxmax_tmp = self.inflow_grid.shape[2] self.floodcri = routing_interface.get_floodcri(nbpt = self.nbpt,\ nwbas=self.nwbas, inflowmax = nbxmax_tmp, inflow_number = self.inflow_number,\ inflow_basin = self.inflow_basin, inflow_grid = self.inflow_grid, basin_count = self.basin_count,\ basin_floodp = self.basin_floodp, basin_orog_min = self.basin_orog_min) # 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 # # #