Version which writes the fetch of each HTU to the netCDF file

DocumentationInterface

@@ -186,14 +186,14 @@ inflow_number : rank-2 array('i') with bounds (nbpt,8 * nbxmax_in)
 inflow_grid : rank-3 array('i') with bounds (nbpt,8 * nbxmax_in,8 * nbxmax_in)
 inflow_basin : rank-3 array('i') with bounds (nbpt,8 * nbxmax_in,8 * nbxmax_in)
 ====================================================================
-fetch_basin_min,fetch_basin_max = fetch(basin_count,basin_area,basin_id,basin_hierarchy,basin_fac,outflow_grid,outflow_basin,fetch_basin,[nbpt,nwbas])
+basin_area_new,outflow_grid_new = fetch(basin_count,basin_area,basin_id,basin_hierarchy,basin_fac,outflow_grid,outflow_basin,fetch_basin,[nbpt,nwbas])
 
 Wrapper for ``fetch``.
 
 Parameters
 ----------
 basin_count : input rank-1 array('i') with bounds (nbpt)
-basin_area : in/output rank-2 array('f') with bounds (nbpt,nwbas)
+basin_area : input rank-2 array('f') with bounds (nbpt,nwbas)
 basin_id : input rank-2 array('i') with bounds (nbpt,nwbas)
 basin_hierarchy : input rank-2 array('f') with bounds (nbpt,nwbas)
 basin_fac : input rank-2 array('f') with bounds (nbpt,nwbas)
@@ -210,10 +210,10 @@ nwbas : input int, optional
 
 Returns
 -------
-fetch_basin_min : float
-fetch_basin_max : float
+basin_area_new : rank-2 array('f') with bounds (nbpt,nwbas)
+outflow_grid_new : rank-2 array('i') with bounds (nbpt,nwbas)
 ====================================================================
-routing_area,routing_cg,topo_resid,route_togrid,route_tobasin,route_nbintobas,global_basinid,route_outlet,route_type,origin_nbintobas = truncate(nbasmax,basin_count,basin_notrun,basin_area,basin_cg,basin_topoind,fetch_basin,basin_id,basin_coor,basin_type,basin_flowdir,outflow_grid,outflow_basin,inflow_number,inflow_grid,inflow_basin,[nbpt,nbxmax_in,nwbas])
+routing_area,routing_cg,topo_resid,route_togrid,route_tobasin,route_nbintobas,global_basinid,route_outlet,route_type,origin_nbintobas,routing_fetch = truncate(nbasmax,basin_count,basin_notrun,basin_area,basin_cg,basin_topoind,fetch_basin,basin_id,basin_coor,basin_type,basin_flowdir,outflow_grid,outflow_basin,inflow_number,inflow_grid,inflow_basin,[nbpt,nbxmax_in,nwbas])
 
 Wrapper for ``truncate``.
 
@@ -257,3 +257,4 @@ global_basinid : rank-2 array('i') with bounds (nbpt,nbasmax)
 route_outlet : rank-3 array('f') with bounds (nbpt,nbasmax,2)
 route_type : rank-2 array('f') with bounds (nbpt,nbasmax)
 origin_nbintobas : rank-1 array('i') with bounds (nbpt)
+routing_fetch : rank-2 array('f') with bounds (nbpt,nbasmax)

DomainDecompTests/DomainDecomp.py

@@ -129,8 +129,8 @@ if rank == 3 :
     print yy
 
 print rank, "Sum over Core regions :", part.sumcore(fx,order='F')/len(part.landcorelist)
-    
-part.zerocore(fx,order='F')
+
+fx = part.zerocore(fx,order='F')
 
 if np.sum(part.sumcore(fx,order='F')) > 0 :
     print "Error the core region was not set to zero on processor ", rank

F90subroutines/routing_interface.f90

@@ -287,7 +287,7 @@ SUBROUTINE linkup(nbpt, nbxmax_in, nwbas, basin_count, basin_area, basin_id, bas
 END SUBROUTINE linkup
 
 SUBROUTINE fetch(nbpt, nwbas, basin_count, basin_area, basin_id, basin_hierarchy, basin_fac,&
-     & outflow_grid, outflow_basin, fetch_basin, fetch_basin_min, fetch_basin_max)
+     & outflow_grid, outflow_basin, fetch_basin, basin_area_new, outflow_grid_new)
   !
   USE ioipsl
   USE grid
@@ -300,7 +300,7 @@ SUBROUTINE fetch(nbpt, nwbas, basin_count, basin_area, basin_id, basin_hierarchy
   INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
   INTEGER(i_std), INTENT (in) :: nwbas !!
   INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count !!
-  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_area !!
+  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_area !!
   INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_id !!
   REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in)       :: basin_hierarchy
   REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in)       :: basin_fac
@@ -310,12 +310,14 @@ SUBROUTINE fetch(nbpt, nwbas, basin_count, basin_area, basin_id, basin_hierarchy
   !! IN-OUTPUT VARIABLES
   REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: fetch_basin !!
   ! Output
-  REAL(r_std), INTENT(out) :: fetch_basin_min, fetch_basin_max
-  
-  CALL routing_reg_fetch(nbpt, area, contfrac, nwbas, basin_count, basin_area, basin_id, basin_hierarchy, basin_fac,&
-       & outflow_grid, outflow_basin, fetch_basin)
-  fetch_basin_min = MINVAL(fetch_basin)
-  fetch_basin_max = MAXVAL(fetch_basin)
+  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(out) :: basin_area_new !!
+  INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(out) :: outflow_grid_new !! Type of outflow on the grid box (unitless)
+  !
+  basin_area_new(:,:) = basin_area(:,:)
+  outflow_grid_new(:,:) = outflow_grid(:,:)
+  !
+  CALL routing_reg_fetch(nbpt, area, contfrac, nwbas, basin_count, basin_area_new, basin_id, basin_hierarchy, basin_fac,&
+       & outflow_grid_new, outflow_basin, fetch_basin)
   
   !
 END SUBROUTINE fetch
@@ -324,7 +326,7 @@ END SUBROUTINE fetch
 SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, basin_count, basin_notrun, basin_area, basin_cg, basin_topoind,&
      & fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, inflow_number,&
      & inflow_grid, inflow_basin, routing_area, routing_cg, topo_resid, route_togrid, route_tobasin, route_nbintobas, &
-     & global_basinid, route_outlet, route_type, origin_nbintobas)
+     & global_basinid, route_outlet, route_type, origin_nbintobas, routing_fetch)
   !
   USE ioipsl
   USE grid
@@ -360,6 +362,7 @@ SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, basin_count, basin_notrun,
   ! Output variables
   !
   REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: routing_area !! Surface of basin (m^2)
+  REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: routing_fetch !! Upstream are flowing into HTU (m^2)
   REAL(r_std), DIMENSION(nbpt,nbasmax,2), INTENT(out)  :: routing_cg !! Centre of gravity of HTU (Lat, Lon)
   REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: topo_resid !! Topographic index of the retention time (m)
   INTEGER(i_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: route_togrid !! Grid into which the basin flows (unitless)
@@ -384,6 +387,9 @@ SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, basin_count, basin_notrun,
   topo_resid(:,:) = topo_resid_glo(:,:) 
   route_togrid(:,:) = route_togrid_glo(:,:)
   route_tobasin(:,:) = route_tobasin_glo(:,:)
+
+  routing_fetch(:,:) = route_fetch_glo(:,:)
+  
   route_nbintobas(:) = route_nbintobas_glo(:)
   global_basinid(:,:) = global_basinid_glo(:,:)
   route_outlet(:,:,:) = route_outlet_glo(:,:,:)

F90subroutines/routing_reg.f90

@@ -25,6 +25,7 @@ MODULE routing_reg
   INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: global_basinid_glo !! ID of basin (unitless)
   REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: route_outlet_glo !! Coordinate of outlet (-)
   REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: route_type_glo !! Coordinate of outlet (-)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: route_fetch_glo !! Upstream area at each HTU
   
 CONTAINS
 !! ================================================================================================================================
@@ -3279,6 +3280,7 @@ SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, basin_
        route_togrid_glo(ib,:) = ib
        route_tobasin_glo(ib,:) = 0
        routing_area_glo(ib,:) = zero
+       route_fetch_glo(ib,:) = zero
        route_nbintobas_glo(ib) = basin_count(ib)
        origin_nbintobas_glo(ib) = basin_notrun(ib)
        !
@@ -3296,6 +3298,8 @@ SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, basin_
           route_outlet_glo(ib,ij,1) = basin_coor(ib,ij,1)
           route_outlet_glo(ib,ij,2) = basin_coor(ib,ij,2)
           !
+          route_fetch_glo(ib,ij) = fetch_basin(ib,ij)
+          !
           route_type_glo(ib,ij) = basin_type(ib,ij)
           !
           route_togrid_glo(ib,ij) = outflow_grid(ib,ij)
@@ -3631,6 +3635,7 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count,
     nbasmax_save = nbasmax
     
     ALLOCATE (routing_area_glo(nbpt,nbasmax), stat=ier)
+    ALLOCATE (route_fetch_glo(nbpt,nbasmax), stat=ier)
     ALLOCATE (routing_cg_glo(nbpt,nbasmax,2), stat=ier)
     ALLOCATE (topo_resid_glo(nbpt,nbasmax), stat=ier)
     ALLOCATE (route_togrid_glo(nbpt,nbasmax), stat=ier)

Interface.py

@@ -160,12 +160,11 @@ class HydroSuper :
                                                                        self.nbcoastal, self.coastal_basin, float(hydrodata.basinsmax))
         return
     #
-    def fetch(self) :
+    def fetch(self, part) :
         #
         self.fetch_basin = np.zeros(self.basin_area.shape, dtype=np.float32, order='F')
-        self.fetch_basin_min, self.fetch_basin_max = routing_interface.fetch(self.basin_count, self.basin_area, self.basin_id, self.basin_hierarchy, \
-                                                                             self.basin_fac, self.outflow_grid, self.outflow_basin, self.fetch_basin)
-        print("## Extrema of basin fetch :", self.fetch_basin_min, self.fetch_basin_max/1.0e6)
+        self.basin_area, self.outflow_grid  = routing_interface.fetch(self.basin_count, self.basin_area, self.basin_id, self.basin_hierarchy, \
+                                                                        self.basin_fac, self.outflow_grid, self.outflow_basin, self.fetch_basin)
         return
 #
 #
@@ -174,7 +173,7 @@ class HydroGraph :
     def __init__(self, nbasmax, hydrosuper) :
         self.nbasmax = nbasmax
         self.routing_area, self.routing_cg, self.topo_resid, self.route_togrid, self.route_tobasin, self.route_nbintobas, self.global_basinid, \
-            self.route_outlet, self.route_type, self.origin_nbintobas = \
+            self.route_outlet, self.route_type, self.origin_nbintobas, self.routing_fetch = \
                                     routing_interface.truncate(nbasmax, hydrosuper.basin_count, hydrosuper.basin_notrun, hydrosuper.basin_area, \
                                                                 hydrosuper.basin_cg, hydrosuper.basin_topoind, hydrosuper.fetch_basin, hydrosuper.basin_id, \
                                                                 hydrosuper.basin_outcoor, hydrosuper.basin_type, hydrosuper.basin_flowdir, \
@@ -407,6 +406,19 @@ class HydroGraph :
         CG_lon[:,:,:] = part.gather(cg[1,:,:,:])
         CG_lat[:,:,:] = part.gather(cg[0,:,:,:])
         #
+        # Save the fetch of each basin
+        #
+        fe =  procgrid.landscatter(self.routing_fetch[:,:], order='F')
+        fe = fe.astype(vtyp, copy=False)
+        fe[np.isnan(fe)] = NCFillValue
+        if part.rank == 0 :
+            fetch = outnf.createVariable("fetch", vtyp, ('htu','y','x'), fill_value=NCFillValue)
+            fetch.title = "Fetch contributing to each HTU"
+            fetch.units = "m^2"
+        else :
+            fetch = np.zeros((1,1,1))
+        fetch[:,:,:] = part.gather(fe)
+        #
         if part.rank == 0 :
             outnf.close()
         #

Partition.py

@@ -434,24 +434,47 @@ class partition :
     # Set to zero all points in the core
     #
     def zerocore(self, x, order='C') :
+        if order == 'C' :
+            if x.shape[-1] != self.nbland :
+                ERROR("Space dimension does not correspond to number of land points (case C)")
+                sys.exit()
+        elif order == 'F' :
+            if x.shape[0] != self.nbland :
+                ERROR("Space dimension does not correspond to number of land points (case F)")
+                sys.exit()
+        else :
+            ERROR("Unforessen order for variable to be set to zero over core region.")
+            sys.exit()
+        #
+        y = x
+        #
         if len(x.shape) == 1 :
-            x[self.landcorelist] = 0.0
+            y[self.landcorelist] = 0.0
         elif len(x.shape) == 2 :
             if order == 'C' :
-                x[:,self.landcorelist] = 0.0
+                y[:,self.landcorelist] = 0.0
             elif order == 'F' :
-                x[self.landcorelist,:] = 0.0
-            else :
-                ERROR("Unforessen order for variable to be set to zero over core region.")
-                sys.exit()
+                y[self.landcorelist,:] = 0.0
         else :
             ERROR("Unforessen rank for variable to summedibe set to zero over the core region.")
             sys.exit()
-        return
+        return y
     #
     # Function to sum over the core regions
     #
     def sumcore(self, x, order='C') :
+        if order == 'C' :
+            if x.shape[-1] != self.nbland :
+                ERROR("Space dimension does not correspond to number of land points (case C)")
+                sys.exit()
+        elif order == 'F' :
+            if x.shape[0] != self.nbland :
+                ERROR("Space dimension does not correspond to number of land points (case F)")
+                sys.exit()
+        else :
+            ERROR("Unforessen order for variable to summed over the core region.")
+            sys.exit()
+        #
         if len(x.shape) == 1 :
             y = np.nansum(x[self.landcorelist])
         elif len(x.shape) == 2 :
@@ -459,9 +482,6 @@ class partition :
                y =  np.nansum(x[:,self.landcorelist], axis=1)
             elif order == 'F' :
                y = np.nansum(x[self.landcorelist,:], axis=0)
-            else :
-                ERROR("Unforessen order for variable to summed over the core region.")
-                sys.exit()
         else :
             ERROR("Unforessen rank for variable to summed over the core region.")
             sys.exit()

RoutingPreProc.py

@@ -172,7 +172,7 @@ del hoverlap
 gc.collect()
 
 print("=================== Compute fetch ====================")
-hsuper.fetch()
+hsuper.fetch(part)
 
 hgraph = IF.HydroGraph(nbasmax, hsuper)