Add the distribution of the time step limitation for all HTUs. Serves as a...

Add the distribution of the time step limitation for all HTUs. Serves as a diagnostic for the most approriate time step.

GraphHydro.py

@@ -140,7 +140,7 @@ class HydroGraph :
                         else :
                             INFO(typename+" number : "+str(lid)+" No appropriate HTU with an upstream "+\
                                  "error less than "+str(MaxUpstrErr)+" % ")
-                        #
+    #
     #
     #
     def add_variable(self, outnf, procgrid, NCFillValue, part, coord, name, title, units, data, vtyp, orig_type = "float", arrayorder = 'F'):
@@ -160,6 +160,28 @@ class HydroGraph :
         else :
             ncvar = np.zeros((1,)*len(coord))
         ncvar[:] = part.gather(var, default=NCFillValue)
+        return
+    #
+    #
+    #
+    def add_tstepdistrib(self, outnf, procgrid, NCFillValue, part, data, vtyp, nbbins, tcst, arrayorder = 'F') :
+
+        
+        var = procgrid.landscatter(data.astype(vtyp), order=arrayorder)
+        gvar = part.gather(var, default=NCFillValue)
+        
+        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))
+            ncbins = outnf.createVariable("timebin", vtyp, "nbbins")
+            ncbins.title = "Time bins"
+            ncbins.units = "s"
+            ncbins[:] = (bin_edges[0:-1]+bin_edges[1:])/2.0
+            ncvar = outnf.createVariable("stabtstep", vtyp, "nbbins", fill_value=NCFillValue)
+            ncvar.title = "Stable timestep distribution"
+            ncvar.units = "count"
+            ncvar[:] = hist[:]
+        return
     #
     #
     #
@@ -168,6 +190,7 @@ class HydroGraph :
         vtyp=np.float64
         cornerind=[0,2,4,6]
         nbcorners = len(cornerind)
+        nbbins = 800
         #
         nlmax, nblocated = locations.maxlocations(self.nbasmax, part)
         #
@@ -181,6 +204,7 @@ class HydroGraph :
             outnf.createDimension('bnd', nbcorners)
             outnf.createDimension('inflow', self.max_inflow)
             outnf.createDimension('stnperhtu', nlmax)
+            outnf.createDimension('nbbins', nbbins)
             outnf.createDimension('locations', nblocated)
         else :
             outnf = None
@@ -194,7 +218,8 @@ class HydroGraph :
         nbpt_loc = np.zeros((self.nbpt,1)).astype(np.int32)
         nbpt_loc[:,0] = np.arange(1, self.nbpt+1)
         nbpt_glo = part.l2glandindex(nbpt_loc)
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), "nbpt_glo", "Land point indices on global grid", "-", nbpt_glo[:,0], vtyp)
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), \
+                          "nbpt_glo", "Land point indices on global grid", "-", nbpt_glo[:,0], vtyp)
         #
         ################
         #
@@ -203,68 +228,91 @@ class HydroGraph :
         #
         #
         # The field route_togrid is with indices on the local grid. That needs to be converted to the global grid.
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "routetogrid", "Grid into which the basin flows", "-", grgrid, vtyp, orig_type="int")
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
+                          "routetogrid", "Grid into which the basin flows", "-", grgrid, vtyp, orig_type="int")
         # route to basin
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "routetobasin", "Basin in to which the water goes", "-", self.route_tobasin[:,:], vtyp, orig_type="int")
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
+                          "routetobasin", "Basin in to which the water goes", "-", self.route_tobasin[:,:], vtyp, orig_type="int")
         # basin id
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "basinid", "ID of basin", "-", self.global_basinid[:,:], vtyp, orig_type="int")
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
+                          "basinid", "ID of basin", "-", self.global_basinid[:,:], vtyp, orig_type="int")
         #
         # basin area
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "basin_area", "area of basin", "m^2", self.routing_area[:,:], vtyp)
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
+                          "basin_area", "area of basin", "m^2", self.routing_area[:,:], vtyp)
 
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "basin_orog_mean", "Mean orography", "m", self.routing_orog_mean[:,:], vtyp)
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
+                          "basin_orog_mean", "Mean orography", "m", self.routing_orog_mean[:,:], vtyp)
 
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "basin_orog_min", "Min orography", "m", self.routing_orog_min[:,:], vtyp)
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
+                          "basin_orog_min", "Min orography", "m", self.routing_orog_min[:,:], vtyp)
 
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "basin_orog_max", "Max orography", "m", self.routing_orog_max[:,:], vtyp)
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
+                          "basin_orog_max", "Max orography", "m", self.routing_orog_max[:,:], vtyp)
 
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "basin_floodp", "Fraction of floodplains", "-", self.routing_floodp[:,:], vtyp)
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
+                          "basin_floodp", "Fraction of floodplains", "-", self.routing_floodp[:,:], vtyp)
 
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "floodcri", "Height at which all basin is flooded", "mm", self.floodcri[:,:], vtyp)
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
+                          "floodcri", "Height at which all basin is flooded", "mm", self.floodcri[:,:], vtyp)
 
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "basin_beta_fp", "Shape of the floodplins HTUs (concave / convex)", "-", self.routing_beta[:,:], vtyp)
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
+                          "basin_beta_fp", "Shape of the floodplins HTUs (concave / convex)", "-", self.routing_beta[:,:], vtyp)
 
         # route number into basin
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), "routenbintobas", "Number of basin into current one", "-", self.route_nbintobas[:], vtyp, orig_type="int")
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), \
+                          "routenbintobas", "Number of basin into current one", "-", self.route_nbintobas[:], vtyp, orig_type="int")
         #
         # original number into basin
-        self.add_variable(outnf, procgrid, NCFillValue, part, ( 'y','x'), "originnbintobas", "Number of sub-grid basin into current one before truncation", "-", \
+        self.add_variable(outnf, procgrid, NCFillValue, part, ( 'y','x'), \
+                          "originnbintobas", "Number of sub-grid basin into current one before truncation", "-", \
                           self.origin_nbintobas[:], vtyp, orig_type="int")
         #
         # latitude of outlet
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "outletlat", "Latitude of Outlet", "degrees north", self.route_outlet[:,:,0], vtyp)
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
+                          "outletlat", "Latitude of Outlet", "degrees north", self.route_outlet[:,:,0], vtyp)
         # longitude of outlet
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "outletlon", "Longitude of Outlet", "degrees east", self.route_outlet[:,:,1], vtyp)
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
+                          "outletlon", "Longitude of Outlet", "degrees east", self.route_outlet[:,:,1], vtyp)
         # type of outlet
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "outlettype", "Type of outlet", "code", self.route_type[:,:], vtyp)
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
+                          "outlettype", "Type of outlet", "code", self.route_type[:,:], vtyp)
         #
         # topographic index
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "topoindex", "Topographic index of the retention time", "km", self.topo_resid[:,:], vtyp)
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
+                          "topoindex", "Topographic index of the retention time", "km", self.topo_resid[:,:], vtyp)
+        self.add_tstepdistrib(outnf, procgrid, NCFillValue, part, self.topo_resid[:,:], vtyp, nbbins, tcst)
         #
 
         # Inflow number
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "route_innum", "Number of inflow", "-", self.route_innum[:,:], vtyp, orig_type="int")
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
+                          "route_innum", "Number of inflow", "-", self.route_innum[:,:], vtyp, orig_type="int")
         #
         # Inflow grid
         #gingrid = part.l2glandindex(self.inflow_grid[:,:,:inflow_size])
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow','z','y','x'), "route_ingrid", "Grid from which the water flows", "-", \
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow','z','y','x'), \
+                          "route_ingrid", "Grid from which the water flows", "-", \
                           self.route_ingrid[:,:,:], vtyp, orig_type="int")
         #
         # Inflow basin
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow','z','y','x'), "route_inbasin", "Basin from which the water flows", "-", \
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow','z','y','x'), \
+                          "route_inbasin", "Basin from which the water flows", "-", \
                           self.route_inbasin[:,:,:], vtyp, orig_type="int")
 
         #
         # Save centre of gravity of HTU
         #
         # Check if it works
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "CG_lon", "Longitude of centre of gravity of HTU", "degrees east", self.routing_cg[:,:,1], vtyp)
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
+                          "CG_lon", "Longitude of centre of gravity of HTU", "degrees east", self.routing_cg[:,:,1], vtyp)
 
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "CG_lat", "Latitude of centre of gravity of HTU", "degrees north", self.routing_cg[:,:,0], vtyp)
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
+                          "CG_lat", "Latitude of centre of gravity of HTU", "degrees north", self.routing_cg[:,:,0], vtyp)
         #
         # Save the fetch of each basin
         #
-        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "fetch", "Fetch contributing to each HTU", "m^2", self.routing_fetch[:,:], vtyp)
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
+                          "fetch", "Fetch contributing to each HTU", "m^2", self.routing_fetch[:,:], vtyp)
         #
         # Build and save a map of all the structures located on the graph
         #

RPPtools.py

@@ -7,6 +7,7 @@ FillValue=np.nan
 IntFillValue=-999999
 vtyp=np.float64
 vinttyp=np.int32
+OneDay=86400
 #
 from netCDF4 import Dataset
 import pickle