From f729ce163e99633d3fb4e0fc6540c030d6967092 Mon Sep 17 00:00:00 2001
From: Jan Polcher <jan.polcher@lmd.jussieu.fr>
Date: Sat, 22 Jun 2019 10:22:44 +0200
Subject: [PATCH] Version which writes the fetch of each HTU to the netCDF file
---
DocumentationInterface | 11 ++++----
DomainDecompTests/DomainDecomp.py | 4 +--
F90subroutines/routing_interface.f90 | 24 ++++++++++-------
F90subroutines/routing_reg.f90 | 5 ++++
Interface.py | 22 +++++++++++----
Partition.py | 40 +++++++++++++++++++++-------
RoutingPreProc.py | 2 +-
7 files changed, 76 insertions(+), 32 deletions(-)
diff --git a/DocumentationInterface b/DocumentationInterface
index 1d2bdf4..58d11e5 100644
--- a/DocumentationInterface
+++ b/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)
diff --git a/DomainDecompTests/DomainDecomp.py b/DomainDecompTests/DomainDecomp.py
index ccfbddf..8281e5c 100644
--- a/DomainDecompTests/DomainDecomp.py
+++ b/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
diff --git a/F90subroutines/routing_interface.f90 b/F90subroutines/routing_interface.f90
index 3af02f2..1d6a2e8 100644
--- a/F90subroutines/routing_interface.f90
+++ b/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(:,:,:)
diff --git a/F90subroutines/routing_reg.f90 b/F90subroutines/routing_reg.f90
index cb9d66e..c2f1c3f 100644
--- a/F90subroutines/routing_reg.f90
+++ b/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)
diff --git a/Interface.py b/Interface.py
index e5e7340..32f1917 100644
--- a/Interface.py
+++ b/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()
#
diff --git a/Partition.py b/Partition.py
index ec0333d..d7dc356 100644
--- a/Partition.py
+++ b/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()
diff --git a/RoutingPreProc.py b/RoutingPreProc.py
index 63aa0a0..0c6440a 100644
--- a/RoutingPreProc.py
+++ b/RoutingPreProc.py
@@ -172,7 +172,7 @@ del hoverlap
gc.collect()
print("=================== Compute fetch ====================")
-hsuper.fetch()
+hsuper.fetch(part)
hgraph = IF.HydroGraph(nbasmax, hsuper)
--
GitLab