Commit 2da2d684 authored by POLCHER Jan's avatar POLCHER Jan 🚴🏾
Browse files

Time constantes are now in s/km in the routing-graph file.

parent a2f12108
......@@ -37,19 +37,19 @@ def addparameters(outnf, version, part, tcst, vtyp, NCFillValue) :
outnf.createDimension('dimpara', 1)
stream = outnf.createVariable("StreamTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
stream.title = "Time constant for the stream reservoir"
stream.units = "10+3.day/km"
stream.units = "s/km"
stream[:] = tcst.stream_tcst
fast = outnf.createVariable("FastTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
fast.title = "Time constant for the fast reservoir"
fast.units = "10+3.day/km"
fast.units = "s/km"
fast[:] = tcst.fast_tcst
slow = outnf.createVariable("SlowTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
slow.title = "Time constant for the slow reservoir"
slow.units = "10+3.day/km"
slow.units = "s/km"
slow[:] = tcst.slow_tcst
flood = outnf.createVariable("FloodTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
flood.title = "Time constant for the flood reservoir"
flood.units = "10+3.day/km"
flood.units = "s/km"
flood[:] = tcst.flood_tcst
swamp = outnf.createVariable("SwampCst", vtyp, ('dimpara'), fill_value=NCFillValue)
swamp.title = "Fraction of the river transport that flows to the swamps"
......@@ -172,7 +172,7 @@ class HydroGraph :
if part.rank == 0:
gvar[gvar >= NCFillValue] = np.nan
hist, bin_edges = np.histogram(gvar*tcst.stream_tcst/1000*RPP.OneDay, bins=nbbins, range=(0,RPP.OneDay/12))
hist, bin_edges = np.histogram(gvar*tcst.stream_tcst, bins=nbbins, range=(0,RPP.OneDay/12))
ncbins = outnf.createVariable("timebin", vtyp, "nbbins")
ncbins.title = "Time bins"
ncbins.units = "s"
......
......@@ -133,12 +133,13 @@ def calctopoindex(nf, istr, iend, jstr, jend, index, hydrogrid) :
dl = (rx+ry)*0.5*0.5
topoind[j,i]=np.sqrt(dl**3./(dz*1.0e6))
#
minitopo = np.nanmin(topoind)
if minitopo <= np.finfo(float).eps :
mintopo = np.nanmin(topoind)
maxtopo = np.nanmax(topoind)
if mintopo <= np.finfo(float).eps :
ERROR("Topoind close to zero encoutered.")
sys.exit()
#
return gather(topoind, index, default = 10), "Topographic index which serves for the residence time", minitopo
return gather(topoind, index, default = 10), "Topographic index which serves for the residence time", mintopo, maxtopo
#
class HydroGrid :
def __init__(self, lolacorners, wfile) :
......@@ -253,8 +254,9 @@ class HydroData :
self.topoinddesc = "desc"
mintopo = 1
else:
self.topoind, self.topoinddesc, mintopo = calctopoindex(nf, istr, iend, jstr, jend, index, hydrogrid)
self.topoind, self.topoinddesc, mintopo, maxtopo = calctopoindex(nf, istr, iend, jstr, jend, index, hydrogrid)
self.topoindmin = part.domainmin(mintopo)
self.topoindmax = part.domainmax(maxtopo)
elif "topoind" in nf.variables.keys() and "orog" in nf.variables.keys() :
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, default = 10)
......@@ -265,9 +267,16 @@ class HydroData :
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))
att = getattrcontaining(nf, "topoind", "max")
if len(att) > 0 :
self.topoindmax=att[0]
else :
INFO("We need to scan full file to find minimum topoind over domain")
self.topoindmax=np.max(np.where(nf.variables["topoind"][:,:] < 1.e15))
else :
ERROR("Not enough information in the HydroFile in order to compute topographic residence time.")
sys.exit()
INFO("Range of topoindex in hydrological file :"+str(self.topoindmin)+" < "+str(self.topoindmax))
#
#
if FloodplainsFile is not None:
......@@ -283,26 +292,27 @@ class HydroData :
class HydroParameter :
# Class to set, compute and manage the time constants of the routing scheme.
def __init__(self, hydrogrid) :
# Now values are in s/km ... the conversion is kept for memory.
if hydrogrid.hd >= 0.5 :
# Case of the Fekete/Vörösmarty hydrological data set
self.stream_tcst = 0.24
self.fast_tcst = 3.0
self.slow_tcst = 25.0
self.flood_tcst = 4.0
self.stream_tcst = 0.24/1000*RPP.OneDay
self.fast_tcst = 3.0/1000*RPP.OneDay
self.slow_tcst = 25.0/1000*RPP.OneDay
self.flood_tcst = 4.0/1000*RPP.OneDay
self.swamp_cst = 0.2
elif hydrogrid.hd >= 0.016 :
# Case for MERIT
self.stream_tcst = 0.0014
self.fast_tcst = 0.07
self.slow_tcst = 0.98
self.flood_tcst = 4.0
self.stream_tcst = 0.0014/1000*RPP.OneDay
self.fast_tcst = 0.07/1000*RPP.OneDay
self.slow_tcst = 0.98/1000*RPP.OneDay
self.flood_tcst = 4.0/1000*RPP.OneDay
self.swamp_cst = 0.2
elif hydrogrid.hd >= 0.008 :
# Case for HydroSEHDS
self.stream_tcst = 0.01
self.fast_tcst = 0.5
self.slow_tcst = 7.0
self.flood_tcst = 4.0
self.stream_tcst = 0.01/1000*RPP.OneDay
self.fast_tcst = 0.5/1000*RPP.OneDay
self.slow_tcst = 7.0/1000*RPP.OneDay
self.flood_tcst = 4.0/1000*RPP.OneDay
self.swamp_cst = 0.2
else :
print("For the resolution ",hydrogrid.hd," deg. we do not yet have a parameter set.")
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment