diff --git a/DocumentationInterface b/DocumentationInterface
index d71e5bd0ef33711f1003a996c72e3d97b0e2e288..1d2bdf4278a66098288000d6af072d6d93eda769 100644
--- a/DocumentationInterface
+++ b/DocumentationInterface
@@ -186,7 +186,7 @@ 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 = fetch(basin_count,basin_area,basin_id,basin_hierarchy,basin_fac,outflow_grid,outflow_basin,[nbpt,nwbas])
+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])
Wrapper for ``fetch``.
@@ -199,6 +199,7 @@ basin_hierarchy : input rank-2 array('f') with bounds (nbpt,nwbas)
basin_fac : input rank-2 array('f') with bounds (nbpt,nwbas)
outflow_grid : in/output rank-2 array('i') with bounds (nbpt,nwbas)
outflow_basin : input rank-2 array('i') with bounds (nbpt,nwbas)
+fetch_basin : in/output rank-2 array('f') with bounds (nbpt,nwbas)
Other Parameters
----------------
@@ -209,7 +210,8 @@ nwbas : input int, optional
Returns
-------
-fetch_basin : rank-2 array('f') with bounds (nbpt,nwbas)
+fetch_basin_min : float
+fetch_basin_max : float
====================================================================
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])
diff --git a/DomainDecompTests/DomainDecomp.py b/DomainDecompTests/DomainDecomp.py
index ca0bacfc43cc81ed5a3423184f6c18bedecbeeeb..ccfbddffa119a104e1e6b9e992a55c76e9074ee7 100644
--- a/DomainDecompTests/DomainDecomp.py
+++ b/DomainDecompTests/DomainDecomp.py
@@ -109,12 +109,38 @@ x[jhoff:jhoff+part.nj,ihoff:ihoff+part.ni]=rank
lx = modelgrid.landgather(x)
part.landsendtohalo(lx)
xland = part.landscatgat(modelgrid, lx)
-
+
if rank == 0 :
xland[~np.isnan(xland)] = 1
xland[np.isnan(xland)] = 0
print("Difference between global and gathered version of land only :",np.sum(xland-gg.land))
+#
+# Test a 2D landsendtohalo with a FORTRAN ordering
+#
+fx=np.zeros((modelgrid.nbland,nbc), dtype=np.float32, order='F')
+fx[:,:] = rank
+
+part.landsendtohalo(fx, order='F')
+
+if rank == 3 :
+ print "Before Zerored core on proc 1"
+ yy=modelgrid.landscatter(fx[:,3])
+ print yy
+
+print rank, "Sum over Core regions :", part.sumcore(fx,order='F')/len(part.landcorelist)
+
+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
+ sys.exit()
+
+if rank == 3 :
+ print "Zerored core on proc 3"
+ yy=modelgrid.landscatter(fx[:,3])
+ print yy
+
print("len areas : ", len(modelgrid.area))
areas = part.landscatgat(modelgrid, np.array(modelgrid.area))
if rank == 0 :
@@ -125,4 +151,4 @@ y=np.zeros((nbc,modelgrid.nbland))
print("Test : ", y.shape, modelgrid.nbland)
y[rank,:]=rank
-savefile.dumpnetcdf("TestDomainDecomp.nc", gg, modelgrid, part, rland, x, y)
+savefile.dumpnetcdf("TestDomainDecomp.nc", gg, modelgrid, part, rland, x, y, fx)
diff --git a/DomainDecompTests/savefile.py b/DomainDecompTests/savefile.py
index 33a2dd303526240f3967dc6059ea34386d171eb6..2d328da07281896e13810a39dbd8050807861e57 100644
--- a/DomainDecompTests/savefile.py
+++ b/DomainDecompTests/savefile.py
@@ -4,7 +4,7 @@ from netCDF4 import Dataset
FillValue=1.0e20
-def dumpnetcdf(filename, globalgrid, procgrid, part, x1dland, x1dxy, x2dland) :
+def dumpnetcdf(filename, globalgrid, procgrid, part, x1dland, x1dxy, x2dland, fx) :
#
vtyp=np.float64
cornerind=[0,2,4,6]
@@ -87,7 +87,17 @@ def dumpnetcdf(filename, globalgrid, procgrid, part, x1dland, x1dxy, x2dland) :
x2dl.standard_name="Some 2D land variable"
x2dl.title="2D Land"
x2dl[:,:,:] = fullX2D[:,:,:]
-
+ #
+ # Test a Fortran gathered field
+ #
+ fg = procgrid.landscatter(fx[:,:], order='F')
+ if part.rank == 0 :
+ l2id = outnf.createVariable("land2d", vtyp, ('htu','lat','lon'), fill_value=FillValue)
+ l2id.title = "Land 2d in F order"
+ l2id.units = "-"
+ else :
+ l2id = np.zeros((1,1,1))
+ l2id[:,:,:] = part.gather(fg)
#
# Proc by land
if part.rank == 0 :
diff --git a/F90subroutines/routing_interface.f90 b/F90subroutines/routing_interface.f90
index f7bba7d9bb62f432b41f1725d3d4c8fb3f4ec26c..3af02f2d00b58db41fe458add451e00863f28bd4 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)
+ & outflow_grid, outflow_basin, fetch_basin, fetch_basin_min, fetch_basin_max)
!
USE ioipsl
USE grid
@@ -307,11 +307,16 @@ SUBROUTINE fetch(nbpt, nwbas, basin_count, basin_area, basin_id, basin_hierarchy
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid !! Type of outflow on the grid box (unitless)
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: outflow_basin !!
!
- !! OUTPUT VARIABLES
- REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(out) :: fetch_basin !!
+ !! 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)
+
!
END SUBROUTINE fetch
diff --git a/F90subroutines/routing_reg.f90 b/F90subroutines/routing_reg.f90
index 81a97d171fcfe206e68c2dfaba69ecea0d275e70..cb9d66e5d0b472e826d513240df9020d66a2ae4b 100644
--- a/F90subroutines/routing_reg.f90
+++ b/F90subroutines/routing_reg.f90
@@ -2738,10 +2738,6 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, basin_count, basin
!
ENDDO
!
- ! Compute the area upstream of each basin
- !
- fetch_basin(:,:) = zero
- !
!
DO ib=1,nbpt
!
diff --git a/Interface.py b/Interface.py
index be95db1e9adf49db9a004eacf79b9de21cb48c92..e5e73406a333bc19c25883b0c95d15c16477be30 100644
--- a/Interface.py
+++ b/Interface.py
@@ -161,8 +161,11 @@ class HydroSuper :
return
#
def fetch(self) :
- self.fetch_basin = 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 = 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)
return
#
#
diff --git a/ModelGrid.py b/ModelGrid.py
index 0f2036b37666432a79b55f37a814bf2229fa6601..a253481756013228f478c7e5f6274ff910a3c489 100644
--- a/ModelGrid.py
+++ b/ModelGrid.py
@@ -297,7 +297,7 @@ class ModelGrid :
#
def landscatter(self, var, order='C') :
#
- # Some arrays can be in FORTRAN convention and thus the dimension to be scattered is not the last but the first.
+ # Some arrays can be in FORTRAN convention and thus the dimension to be scattered is not the last but the first.
#
dims = var.shape
transpose=False
diff --git a/Partition.py b/Partition.py
index 09a8b0db4c6c80ee6540de5c9acb63b995d6f316..678432781bb2527f8602c227e93cf2dbca78f8fc 100644
--- a/Partition.py
+++ b/Partition.py
@@ -147,6 +147,18 @@ def coresendlist(innersend_map, istart, ihstart, ni, jstart, jhstart, nj) :
innersend[ic][1][:] = innersend[ic][1][:]+ihoff
return sendto, innersend
#
+# All points belonging to the core of the domain
+#
+def listcoreland(istart, ihstart, ni, jstart, jhstart, nj, landind) :
+ corelandlist=[]
+ ihoff=istart-ihstart
+ jhoff=jstart-jhstart
+ for i in range(ni) :
+ for j in range(nj) :
+ if landind[jhoff+j,ihoff+i] >= 0 :
+ corelandlist.append(landind[jhoff+j,ihoff+i])
+ return corelandlist
+#
# Get list of halo points which need to receive data
#
def haloreceivelist(halosource_map, rank) :
@@ -285,6 +297,10 @@ class partition :
for ic in range(self.nbsend) :
wunit.write("To "+str(self.sendto[ic])+" send innersend :"+str(self.innersend[ic])+"\n")
wunit.write("To "+str(self.sendto[ic])+" send landinnersend :"+str(self.landinnersend[ic])+"\n")
+ #
+ self.landcorelist = listcoreland(self.istart, self.ihstart, self.ni, self.jstart, self.jhstart, self.nj, landindmap)
+ #
+ return
#
# Send halo to the other procs
#
@@ -307,13 +323,21 @@ class partition :
#
# For field gathered in land point send halo to the other procs
#
- def landsendtohalo(self, x) :
+ def landsendtohalo(self, x, order='C') :
#
for i in range(max(self.nbreceive,self.nbsend)) :
if i < self.nbsend :
if len(self.landinnersend[i]) > 0 :
if len(x.shape) == 1 :
self.comm.send(x[self.landinnersend[i]], dest=self.sendto[i])
+ elif len(x.shape) == 2 :
+ if order == 'C' :
+ self.comm.send(x[:,self.landinnersend[i]], dest=self.sendto[i])
+ elif order == 'F' :
+ self.comm.send(x[self.landinnersend[i],:], dest=self.sendto[i])
+ else :
+ ERROR("Unforessen order of the variable to be sent to halo of land points")
+ sys.exit()
else :
ERROR("Unforessen rank of the variable to be sent to halo of land points")
sys.exit()
@@ -321,6 +345,14 @@ class partition :
if len(self.landhalosrc[i]) > 0 :
if len(x.shape) == 1 :
x[self.landhalosrc[i]] = self.comm.recv(source=self.receivefrom[i])
+ elif len(x.shape) == 2 :
+ if order == 'C' :
+ x[:,self.landhalosrc[i]] = self.comm.recv(source=self.receivefrom[i])
+ elif order == 'F' :
+ x[self.landhalosrc[i],:] = self.comm.recv(source=self.receivefrom[i])
+ else :
+ ERROR("Unforessen order of the variable to be sent to halo of land points")
+ sys.exit()
else :
ERROR("Unforessen rank of the variable to be received in halo of land points")
sys.exit()
@@ -398,3 +430,34 @@ class partition :
ERROR("The first dimension does not have the length of the number of land points")
sys.exit()
return y
+ #
+ # Set to zero all points in the core
+ #
+ def zerocore(self, x, order='C') :
+ if len(x.shape) == 1 :
+ x[self.landcorelist] = 0.0
+ elif len(x.shape) == 2 :
+ if order == 'C' :
+ x[:,self.landcorelist] = 0.0
+ elif order == 'F' :
+ x[self.landcorelist,:] = 0.0
+ else :
+ ERROR("Unforessen order for variable to zero")
+ sys.exit()
+
+ return
+ #
+ # Function to sum over the core regions
+ #
+ def sumcore(self, x, order='C') :
+ if len(x.shape) == 1 :
+ y = np.nansum(x[self.landcorelist])
+ elif len(x.shape) == 2 :
+ if order == 'C' :
+ 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 zero")
+ sys.exit()
+ return y