Some code clean-up and preparation for the computation of topoindex within the HydroData class.

Topoindex will become optional if orography is present.

F90subroutines/routing_interface.f90

@@ -772,7 +772,7 @@ SUBROUTINE correct_inflows(nbpt, nwbas, inflowmax, outflow_grid,&
   INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax), INTENT(inout) :: inflow_grid
 
   ! LOCAL
-  INTEGER(i_std) :: ig, nbas, ib, og, ob, inf, found
+  INTEGER(i_std) :: ig, nbas, ib, og, ob
 
   WRITE(numout,*) "Checking if the HTUs are in the inflows of their outflow"
 

F90subroutines/routing_reg.f90

@@ -2866,7 +2866,6 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, inflowmax, b
     ! Update the information needed in the basin "totakeover"
     ! For the moment only area
     !
-    WRITE(numout, *) "XXXXX Arguments : inflow_number = ", inflow_number(ib,1:10)
     !
     basin_cg(ib, totakeover, 1) = (basin_cg(ib, totakeover, 1)*basin_area(ib, totakeover) &
     & + basin_cg(ib, tokill, 1)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill))

HydroGrid.py

@@ -93,9 +93,14 @@ class HydroGrid :
 class HydroData :
     def __init__(self, nf, box, index) :
         istr, iend, jstr, jend = box[:]
+        #
+        # Flow direction
+        #
         self.trip=gather(nf.variables["trip"][jstr:jend,istr:iend].astype(np.float32), index, 97)
         self.tripdesc=nf.variables["trip"].long_name
         #
+        # ID of basin
+        #
         self.basins=gather(nf.variables["basins"][jstr:jend,istr:iend], index, 999)
         self.basinsdesc=nf.variables["basins"].long_name
         att = getattrcontaining(nf, "basins", "max")
@@ -106,36 +111,33 @@ class HydroData :
             # This variable seems not to be used further
             self.basinsmax = part.domainmax(np.ma.max(ma.masked_where(self.basins < 1.e10, self.basins)))
         #
-        self.topoind=gather(nf.variables["topoind"][jstr:jend,istr:iend].astype(np.float32), index, 10)
-        self.topoinddesc=nf.variables["topoind"].long_name
-        att = getattrcontaining(nf, "topoind", "min")
-        if len(att) > 0 :
-            self.topoindmin=att[0]
-        else :
-            INFO("We need to scan full file to find minimum topoind over domain")
-            self.topoindmin=np.min(np.where(nf.variables["topoind"][:,:] <  1.e15))
-        #
-        #
-        self.topoindh=gather(nf.variables["topoind_h"][jstr:jend,istr:iend].astype(np.float32), index, 10)
-        self.topoindhdesc=nf.variables["topoind_h"].long_name
-        att = getattrcontaining(nf, "topoind_h", "min")
-        if len(att) > 0 :
-            self.topoindhmin=att[0]
-        else :
-            INFO("We need to scan full file to find minimum topoind_h over domain")
-            self.topoindhmin=np.min(np.where(nf.variables["topoind_h"][:,:] <  1.e15))
+        # Distance to ocean
         #
         self.disto=gather(nf.variables["disto"][jstr:jend,istr:iend].astype(np.float32), index, 0)
         self.distodesc=nf.variables["disto"].long_name
         #
+        # Flow accumulation
+        #
         self.fac=gather(nf.variables["fac"][jstr:jend,istr:iend].astype(np.float32), index, 0)
         self.facdesc=nf.variables["fac"].long_name
         #
+        # If Orography is present in the file then we can compute the topographic index.
+        # Else the topographic index needs to be present in the file.
+        #
         if "orog" in nf.variables.keys():
             self.orog = gather(nf.variables["orog"][jstr:jend,istr:iend].astype(np.float32), index, 0)
             self.orogdesc=nf.variables["orog"].long_name
         else:
             self.orog = gather(np.zeros((jend-jstr,iend-istr)).astype(np.float32), index)
+            self.topoind=gather(nf.variables["topoind"][jstr:jend,istr:iend].astype(np.float32), index, 10)
+            self.topoinddesc=nf.variables["topoind"].long_name
+            att = getattrcontaining(nf, "topoind", "min")
+            if len(att) > 0 :
+                self.topoindmin=att[0]
+            else :
+                INFO("We need to scan full file to find minimum topoind over domain")
+                self.topoindmin=np.min(np.where(nf.variables["topoind"][:,:] <  1.e15))
+        #
         #
         if "floodplains" in nf.variables.keys():
             self.floodplains = gather(nf.variables["floodplains"][jstr:jend,istr:iend].astype(np.float32), index, 0)

RoutingPreProc.py

@@ -63,7 +63,10 @@ w = RPP.compweights(wfile, part, modelgrid, hydrogrid)
 #
 #
 nbpt = len(w.index)
-nbvmax = part.domainmax(max(w.hpts))
+nbvmax = part.domainmax(w.maxhpts)
+#
+nbasmax = min(nbasmax,nbvmax)
+#
 INFO("nbpt : {0}".format(nbpt))
 INFO("nbvmax : {0}".format(nbvmax))
 INFO("nbasmax : {0}".format(nbasmax))