From 55225f03f3a355c71709e389aa3829168203649c Mon Sep 17 00:00:00 2001
From: POLCHER Jan <jan.polcher@lmd.jussieu.fr>
Date: Fri, 5 Jun 2020 15:24:40 +0200
Subject: [PATCH] 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 | 2 +-
F90subroutines/routing_reg.f90 | 1 -
HydroGrid.py | 38 +++++++++++++++-------------
RoutingPreProc.py | 5 +++-
4 files changed, 25 insertions(+), 21 deletions(-)
diff --git a/F90subroutines/routing_interface.f90 b/F90subroutines/routing_interface.f90
index a163c64..3d49719 100644
--- a/F90subroutines/routing_interface.f90
+++ b/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"
diff --git a/F90subroutines/routing_reg.f90 b/F90subroutines/routing_reg.f90
index e569b5b..2e3dada 100644
--- a/F90subroutines/routing_reg.f90
+++ b/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))
diff --git a/HydroGrid.py b/HydroGrid.py
index da9dc9c..0ebb1b0 100644
--- a/HydroGrid.py
+++ b/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)
diff --git a/RoutingPreProc.py b/RoutingPreProc.py
index b3062b8..d9bd597 100644
--- a/RoutingPreProc.py
+++ b/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))
--
GitLab