Correcting some errors

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()):

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: