From dc6837f52c8f56c79f018e7ced6f194c8a48b0c4 Mon Sep 17 00:00:00 2001 From: Anthony <anthony.schrapffer@polytechnique.fr> Date: Mon, 18 May 2020 23:07:32 +0200 Subject: [PATCH] Correcting some errors --- Diags/Routing_Test_and_Correction.py | 2 -- Interface.py | 34 +++++++++++++++++----------- 2 files changed, 21 insertions(+), 15 deletions(-) diff --git a/Diags/Routing_Test_and_Correction.py b/Diags/Routing_Test_and_Correction.py index 298d925..29c527e 100644 --- a/Diags/Routing_Test_and_Correction.py +++ b/Diags/Routing_Test_and_Correction.py @@ -151,7 +151,6 @@ class routing_test: if np.prod(ovar.shape) > 300*300*100: for it in range(0,ovar.shape[0], 5): - print(" ", iit, ',', eit) iit = it eit = it + 5 if varn in list(self.D_corr.keys()): @@ -159,7 +158,6 @@ class routing_test: else: newvar[iit:eit,] = ovar[iit:eit,] if eit != ovar.shape[0]: - print(" ", iit, ',', eit) iit = eit eit = ovar.shape[0] if varn in list(self.D_corr.keys()): diff --git a/Interface.py b/Interface.py index f8f0e67..3e1daaa 100644 --- a/Interface.py +++ b/Interface.py @@ -211,7 +211,7 @@ def finalfetch(part, routing_area, basin_count, route_togrid, route_tobasin, fet # fetch_error = np.sum(np.abs(fetch_out[part.landcorelist,:]-fetch_in[part.landcorelist,:]), axis=1)\ - /np.sum(routing_area[part.landcorelist,:], axis=1) + / np.ma.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) @@ -240,6 +240,18 @@ class HydroOverlap : # 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') @@ -259,22 +271,15 @@ class HydroOverlap : 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 # - # Compute nbxmax - # - ijdim=[] - for ib in range(nbpt) : - ijdim.append(max(np.max(sub_index[ib,:,0])-np.min(sub_index[ib,:,0])+1,np.max(sub_index[ib,:,1])-np.min(sub_index[ib,:,1])+1)) - ijdimmax = max(ijdim) - # # 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.orog_bx, self.floodp_bx, \ self.lon_bx, self.lat_bx, self.lshead_bx = \ @@ -282,7 +287,9 @@ class HydroOverlap : hydrodata.basinsmax, hydrodata.topoindmin, sub_lon, sub_lat, trip_tmp, basins_tmp, topoind_tmp, fac_tmp,\ hierarchy_tmp, orog_tmp, floodp_tmp) # - # Plot some diagnostics for the hydrology grid within the atmospheric meshes. + 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. @@ -619,6 +626,7 @@ class HydroGraph : if orig_type == "float": var[np.isnan(var)] = NCFillValue elif orig_type == "int": + var[np.isnan(var)] = RPP.IntFillValue var[var>=np.abs(RPP.IntFillValue)] = NCFillValue if part.rank == 0: -- GitLab