diff --git a/F90subroutines/routing_reg.f90 b/F90subroutines/routing_reg.f90
index 13bb6509effcbccdc8e0f00f1cb464e937da0a81..52464f9a8c1eb1d45de28180f9394870d08f5001 100644
--- a/F90subroutines/routing_reg.f90
+++ b/F90subroutines/routing_reg.f90
@@ -2073,10 +2073,6 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
! At this point any of the outflows is designated by a negative values in
! outflow_grid
!
- WRITE(numout,*) "*****"
- WRITE(numout,*) "Linkup 0 - sp, sb = ", sp, sb
- WRITE(numout,*) "Linkup 0 - outflow_grid = ", outflow_grid(sp,sb)
- !
found = 0
IF ( outflow_grid(sp,sb) == 0 ) THEN
found = 1
@@ -2130,14 +2126,14 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
& nbpt, nwbas, inp, basin_count, basin_id, basin_hierarchy, basin_fac, &
& basin_flowdir, nbcoastal, coastal_basin, bop, bopqual)
!
- IF ( bop .LT. undef_int ) THEN
+ IF ( bop .LT. undef_int .AND. bop .NE. sb) THEN
!
CALL routing_updateflow(sp, sb, inp, bop, nbpt, nwbas, inflowmax, outflow_grid, outflow_basin, &
& inflow_number, inflow_grid, inflow_basin)
IF ( outflow_basin(sp,sb) == bop ) THEN
solved(sp,1) = solved(sp,1) + 1
found = 1
- WRITE(numout,*) "Solution found in the original outflow_grid"
+ WRITE(numout,*) "Solution found in the original outflow_grid", sb, bop
ENDIF
!
ENDIF
@@ -2190,7 +2186,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
! we know it's wrong but it will serve to make decision later
angp1 = routing_anglediff_g(sp, dp1, basin_flowdir(sp,sb))
! Check if grid point
- IF (dp1 .GT. 0) THEN
+ IF (dp1 .GT. 0 .AND. dp1 .NE. sp) THEN
CALL routing_reg_bestsubbasin(sp, sb, basin_id(sp,sb), basin_hierarchy(sp,sb), basin_fac(sp,sb), &
& basin_flowdir(sp,sb), invented_basins, &
& nbpt, nwbas, dp1, basin_count, basin_id, basin_hierarchy, basin_fac, &
@@ -2205,7 +2201,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
dm1 = neighbours_g(sp,order_ref(nb*2+1))
! we know it's wrong but it will serve to make decision later
angm1 = routing_anglediff_g(sp, dm1, basin_flowdir(sp,sb))
- IF (dm1 .GT. 0) THEN
+ IF (dm1 .GT. 0 .AND. dm1 .NE. sp) THEN
CALL routing_reg_bestsubbasin(sp, sb, basin_id(sp,sb), basin_hierarchy(sp,sb), basin_fac(sp,sb), &
& basin_flowdir(sp,sb), invented_basins, &
& nbpt, nwbas, dm1, basin_count, basin_id, basin_hierarchy, basin_fac, &
@@ -2256,7 +2252,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
IF ( outflow_basin(sp,sb) == bop ) THEN
solved(sp,2) = solved(sp,2) + 1
found = 1
- WRITE (numout,*) "Neighbours, output found at:" , nb, "level"
+ WRITE (numout,*) "Neighbours, output found at:" , nb, "level. sp,sb=", sp, sb, " dop,bop=", dop, bop
ELSE
nb = nb+1
ENDIF
@@ -2340,7 +2336,8 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
IF (( neighbours_g(sp,order_ref(nb)) == minhloc(basin_id(sp,sb),1)) .AND. (found == 0)) THEN
dop = minhloc(basin_id(sp,sb),1)
bop = minhloc(basin_id(sp,sb),2)
- IF ((dop < undef_int ) .AND. (dop .GT. 0) .AND. (bop < undef_int) .AND. (bop .GT. 0)) THEN
+ IF ((dop < undef_int ) .AND. (dop .GT. 0) .AND. (bop < undef_int) .AND. (bop .GT. 0) &
+ & .AND. (dop .NE. sp) .AND. (bop .NE. sb) ) THEN
CALL routing_updateflow(sp, sb, dop, bop, nbpt,nwbas, inflowmax, outflow_grid, outflow_basin, &
& inflow_number, inflow_grid, inflow_basin)
! It is possible that the lowest hierarchy is in two grid cells
@@ -2352,6 +2349,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
ELSE
WRITE(numout,*) "Lowest hierarchy may be in two different grid points"
END IF
+ WRITE(numout,*) "sp,sb = ", sp,sb, " dop,bop = ", dop,bop
WRITE(numout,*) "hierarch of (sp,sb)", basin_hierarchy(sp,sb)
WRITE(numout,*) "Lowest basinid hierarch", basin_hierarchy(dop,bop)
END IF
diff --git a/HydroGrid.py b/HydroGrid.py
index 755383651291c81c66ebe8629365460b54ef8785..da9dc9c341f855107268b44a38a76a87c4939eb1 100644
--- a/HydroGrid.py
+++ b/HydroGrid.py
@@ -83,7 +83,7 @@ class HydroGrid :
DEBUG("# Range Lon :"+str(np.min(self.lon))+" -- "+str(np.max(self.lon)))
DEBUG("# Range Lat :"+str(np.min(self.lat))+" -- "+str(np.max(self.lat)))
#
- if not os.path.exists(wfile):
+ if wfile is None or not os.path.exists(wfile):
self.polyll, self.polylist, self.centers, self.radius, self.index = corners(self.lon, self.lat, hdlon, hdlat)
def select(self, c, r) :
diff --git a/Interface.py b/Interface.py
index d954856dd14f40dfdcbb9341be4d5add6344cc9e..4d32f1c245ca4c7fe8fcad2925d6faa423f5525f 100644
--- a/Interface.py
+++ b/Interface.py
@@ -88,7 +88,7 @@ def addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorne
lon.units="grid box centre degrees east"
lon.title="Longitude"
lon.axis="X"
- lon[:,:] = longitude[:,:]
+ lon[:,:] = globalgrid.lonmat[:,:]
#
# Latitude
latitude = part.gather(procgrid.lat_full)
@@ -98,7 +98,7 @@ def addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorne
lat.standard_name="grid latitude"
lat.title="Latitude"
lat.axis="Y"
- lat[:] = latitude[:,:]
+ lat[:,:] = globalgrid.latmat[:,:]
#
# Bounds of grid box
#
@@ -206,7 +206,7 @@ def finalfetch(part, routing_area, basin_count, route_togrid, route_tobasin, fet
maxdiff_sorted = 0.0
else :
maxdiff_sorted = np.max(np.abs(sorted_outareas[0:largest_pos]-old_sorted))
- old_sorted[:] = sorted_outareas[0:largest_pos]
+ old_sorted[:] = sorted_outareas[0:largest_pos]
iter_count += 1
#
@@ -402,12 +402,12 @@ class HydroSuper :
# Precision in m^2 of the upstream areas when sorting.
sorted_outareas = (np.unique(np.rint(np.array(xtmp)/prec))*prec)[::-1]
# If mono-proc no need to iterate as fetch produces the full result.
+ l = min(sorted_outareas.shape[0],largest_pos)
if part.size == 1 :
maxdiff_sorted = 0.0
else :
- l = min(sorted_outareas.shape[0],largest_pos)
maxdiff_sorted = np.max(np.abs(sorted_outareas[0:largest_pos]-old_sorted[0:l]))
- old_sorted[:l] = sorted_outareas[0:largest_pos]
+ old_sorted[:l] = sorted_outareas[0:largest_pos]
iter_count += 1
self.fetch_basin = np.copy(fetch_basin)
@@ -653,7 +653,7 @@ class HydroGraph :
outnf.createDimension('x', globalgrid.ni)
outnf.createDimension('y', globalgrid.nj)
outnf.createDimension('land', len(procgrid.area))
- outnf.createDimension('htu', self.nbasmax)
+ outnf.createDimension('z', self.nbasmax)
outnf.createDimension('bnd', nbcorners)
outnf.createDimension('inflow', self.max_inflow)
else :
@@ -689,18 +689,18 @@ 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, ('htu', 'y','x'), "routetogrid", "Grid into which the basin flows", "-", grgrid, vtyp, "int")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "routetogrid", "Grid into which the basin flows", "-", grgrid, vtyp, "int")
# route to basin
- self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "routetobasin", "Basin in to which the water goes", "-", self.route_tobasin[:,:], vtyp, "int")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "routetobasin", "Basin in to which the water goes", "-", self.route_tobasin[:,:], vtyp, "int")
# basin id
- self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "basinid", "ID of basin", "-", self.global_basinid[:,:], vtyp, "int")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basinid", "ID of basin", "-", self.global_basinid[:,:], vtyp, "int")
#
# basin area
- self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "basin_area", "area of basin", "m^2", self.routing_area[:,:], vtyp, "float")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basin_area", "area of basin", "m^2", self.routing_area[:,:], vtyp, "float")
- self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "basin_orog", "Mean orography", "m", self.routing_orog[:,:], vtyp, "float")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basin_orog", "Mean orography", "m", self.routing_orog[:,:], vtyp, "float")
- self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "basin_floodp", "area of floodplains", "m^2", self.routing_floodp[:,:], vtyp, "float")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basin_floodp", "area of floodplains", "m^2", self.routing_floodp[:,:], vtyp, "float")
# route number into basin
self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), "routenbintobas", "Number of basin into current one", "-", self.route_nbintobas[:], vtyp, "int")
@@ -709,37 +709,37 @@ class HydroGraph :
self.add_variable(outnf, procgrid, NCFillValue, part, ( 'y','x'), "originnbintobas", "Number of sub-grid basin into current one before truncation", "-", self.origin_nbintobas[:], vtyp, "int")
#
# latitude of outlet
- self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "outletlat", "Latitude of Outlet", "degrees north", self.route_outlet[:,:,0], vtyp, "float")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "outletlat", "Latitude of Outlet", "degrees north", self.route_outlet[:,:,0], vtyp, "float")
# longitude of outlet
- self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "outletlon", "Longitude of Outlet", "degrees east", self.route_outlet[:,:,1], vtyp, "float")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "outletlon", "Longitude of Outlet", "degrees east", self.route_outlet[:,:,1], vtyp, "float")
# type of outlet
- self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "outlettype", "Type of outlet", "code", self.route_type[:,:], vtyp, "float")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "outlettype", "Type of outlet", "code", self.route_type[:,:], vtyp, "float")
#
# topographic index
- self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "topoindex", "Topographic index of the retention time", "m", self.topo_resid[:,:], vtyp, "float")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "topoindex", "Topographic index of the retention time", "m", self.topo_resid[:,:], vtyp, "float")
#
# Inflow number
- self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "route_innum", "Number of inflow", "-", self.route_innum[:,:], vtyp, "int")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "route_innum", "Number of inflow", "-", self.route_innum[:,:], vtyp, "int")
#
# Inflow grid
#gingrid = part.l2glandindex(self.inflow_grid[:,:,:inflow_size])
- self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow', 'htu','y','x'), "route_ingrid", "Grid from which the water flows", "-", self.route_ingrid[:,:,:], vtyp, "int")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow', 'z','y','x'), "route_ingrid", "Grid from which the water flows", "-", self.route_ingrid[:,:,:], vtyp, "int")
#
# Inflow basin
- self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow', 'htu','y','x'), "route_inbasin", "Basin from which the water flows", "-", self.route_inbasin[:,:,:], vtyp, "int")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow', 'z','y','x'), "route_inbasin", "Basin from which the water flows", "-", self.route_inbasin[:,:,:], vtyp, "int")
#
# Save centre of gravity of HTU
#
# Check if it works
- self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "CG_lon", "Longitude of centre of gravity of HTU", "degrees east", self.routing_cg[:,:,1], vtyp, "float")
+ 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, "float")
- self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "CG_lat", "Latitude of centre of gravity of HTU", "degrees north", self.routing_cg[:,:,0], vtyp, "float")
+ 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, "float")
#
# Save the fetch of each basin
#
- self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "fetch", "Fetch contributing to each HTU", "m^2", self.routing_fetch[:,:], vtyp, "float")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "fetch", "Fetch contributing to each HTU", "m^2", self.routing_fetch[:,:], vtyp, "float")
#
# Close the file
if part.rank == 0:
diff --git a/ModelGrid.py b/ModelGrid.py
index a9e01654155b3241796822b828b0bd3b3ae1f45e..f5e3269a3832303ca4379bcdbdeeac498f99ce1c 100644
--- a/ModelGrid.py
+++ b/ModelGrid.py
@@ -94,7 +94,7 @@ def corners(indF, proj, istart, jstart, lon, lat, res_lon, res_lat) :
# Extract the coordinates of the grid file. If ni < 0 and nj < 0 then the
# full grid is read.
#
-def getcoordinates(geo, istart, ni, jstart, nj, res = False) :
+def getcoordinates(geo, istart, ni, jstart, nj) :
#
# Guess the type of grid
#
@@ -138,8 +138,12 @@ def getcoordinates(geo, istart, ni, jstart, nj, res = False) :
#
lon_full=np.copy(geo.variables["XLONG_M"][0,jstart:jstart+nj,istart:istart+ni])
lat_full=np.copy(geo.variables["XLAT_M"][0,jstart:jstart+nj,istart:istart+ni])
- return griddesc, lon_full, lat_full
-
+ #
+ # Compute resolutions
+ #
+ res_lon = np.mean(np.gradient(lon_full, axis=1))
+ res_lat = np.mean(np.gradient(lat_full, axis=0))
+ #
elif griddesc['type'] == "RegLonLat" :
# We have a regulat lat/lon grid
nbi_g = geo.variables["lon"].shape[0]
@@ -168,14 +172,11 @@ def getcoordinates(geo, istart, ni, jstart, nj, res = False) :
lat_full=np.transpose(np.tile(np.copy(geo.variables["lat"][jstart:jstart+nj]),(ni,1)))
griddesc['inilat'] = np.copy(geo.variables["lat"][0])
#
- if res:
- return griddesc, lon_full, lat_full, res_lon, res_lat
- else:
- return griddesc, lon_full, lat_full
-
else :
ERROR("Unknown grid type")
sys.exit()
+ #
+ return griddesc, lon_full, lat_full, res_lon, res_lat
#
#
#
@@ -183,6 +184,8 @@ def getland (geo, ist, ni, jst, nj) :
vn=list(v.name for v in geo.variables.values())
if "LANDMASK" in vn :
land=geo.variables["LANDMASK"][0,jst:jst+nj,ist:ist+ni]
+ elif "Contfrac" in vn :
+ land=geo.variables["Contfrac"][jst:jst+nj,ist:ist+ni]
elif "elevation" in vn :
land=geo.variables["elevation"][jst:jst+nj,ist:ist+ni]
if "missing_value" in geo.variables["elevation"].ncattrs() :
@@ -210,20 +213,16 @@ def getland (geo, ist, ni, jst, nj) :
#
class ModelGrid :
def __init__(self, istart, ni, jstart, nj) :
- # Uncompatibility with Partition over Iberia_regular
- """
- if ni < 2 or nj < 2 :
- INFO("Found impossibleDomain too small for ModelGrid to work : "+str(ni)+str(nj))
- ERROR("Domain too small")
- sys.exit()
- """
+ #
+ #
#
self.source = config.get("OverAll", "ModelGridFile")
geo=Dataset(self.source,'r')
#
# Get the coordinates from the grid file.
#
- griddesc, self.lon_full, self.lat_full = getcoordinates(geo, istart, ni, jstart, nj)
+ griddesc, self.lon_full, self.lat_full, self.res_lon, self.res_lat = getcoordinates(geo, istart, ni, jstart, nj)
+ self.nj,self.ni = self.lon_full.shape
#
# Extract the land/ea mask.
#
@@ -232,15 +231,6 @@ class ModelGrid :
indFi=[[i+1] for j in range(nj) for i in range(ni)]
indFj=[[j+1] for j in range(nj) for i in range(ni)]
#
- # Define some grid variables.
- #
- if griddesc['type'] == "RegLonLat":
- griddesc, self.lon_full, self.lat_full, self.res_lon, self.res_lat = getcoordinates(geo, istart, ni, jstart, nj, res = True)
- else:
- self.res_lon = np.mean(np.gradient(self.lon_full, axis=1))
- self.res_lat = np.mean(np.gradient(self.lat_full, axis=0))
- self.nj,self.ni = self.lon_full.shape
- #
# Gather all land points
#
self.nbland, self.coordll,self.neighbours,self.indP,indF = gatherland(self.lon_full,self.lat_full,self.land,ind,\
@@ -400,7 +390,7 @@ class GlobalGrid :
INFO("Opening : "+self.source)
geo=Dataset(self.source,'r')
#
- griddesc, lon_full, lat_full = getcoordinates(geo, 0, -1, 0, -1)
+ griddesc, lon_full, lat_full, res_lon, res_lat = getcoordinates(geo, 0, -1, 0, -1)
#
# Default behaviour if global is requested in the configuration file.
#
@@ -409,6 +399,8 @@ class GlobalGrid :
self.nj, self.ni = lon_full.shape
self.jgstart = 0
self.igstart = 0
+ self.lonmat = lon_full
+ self.latmat = lat_full
else :
dist=np.sqrt(np.power(lon_full-np.min(lonrange),2) + np.power(lat_full-np.min(latrange),2))
jmin,imin=np.unravel_index(np.argmin(dist, axis=None), dist.shape)
@@ -419,7 +411,10 @@ class GlobalGrid :
self.jgstart = jmin
self.ni = imax-imin+1
self.igstart = imin
- print("Shape : ", lon_full.shape, "N = ", self.nj, self.ni, " Start :", self.jgstart, self.igstart)
+ self.lonmat = lon_full[self.jgstart:self.jgstart+self.nj,self.igstart:self.igstart+self.ni]
+ self.latmat = lat_full[self.jgstart:self.jgstart+self.nj,self.igstart:self.igstart+self.ni]
+ INFO("Shape of full grid kept in GlobalGrid : "+str(self.lonmat.shape))
+
self.land = getland(geo, self.igstart, self.ni, self.jgstart, self.nj)
self.nbland = int(np.sum(self.land))
diff --git a/Partition.py b/Partition.py
index 8038e5ae9aafca5e4e2452c85629b7c742eff72c..db43384748809ea6b1b404d1e8b8c25b57b18eff 100644
--- a/Partition.py
+++ b/Partition.py
@@ -20,7 +20,8 @@ def halfpartition(partin, land) :
if dom["nbi"] > dom["nbj"] :
hh = [h for h in range(dom["nbi"])]
- nb = np.array([np.ma.sum(np.ma.filled(land[dom["jstart"]:dom["jstart"]+dom["nbj"],dom["istart"]:dom["istart"]+h],0)) - dom['nbland']/2. for h in hh])
+ nb = np.array([np.ma.sum(np.ma.filled(land[dom["jstart"]:dom["jstart"]+dom["nbj"],\
+ dom["istart"]:dom["istart"]+h],0)) - dom['nbland']/2. for h in hh])
i = np.nanargmin(np.abs(nb))
new_dom = {"nbi":dom["nbi"]-hh[i],"nbj":dom["nbj"], \
"istart":dom["istart"]+hh[i],"jstart":dom["jstart"]}
@@ -28,13 +29,15 @@ def halfpartition(partin, land) :
else :
hh = [h for h in range(dom["nbj"])]
- nb = np.array([np.ma.sum(np.ma.filled(land[dom["jstart"]:dom["jstart"]+h,dom["istart"]:dom["istart"]+dom["nbi"]],0)) - int(dom['nbland']/2.) for h in hh])
+ nb = np.array([np.ma.sum(np.ma.filled(land[dom["jstart"]:dom["jstart"]+h,\
+ dom["istart"]:dom["istart"]+dom["nbi"]],0)) - int(dom['nbland']/2.) for h in hh])
i = np.nanargmin(np.abs(nb))
new_dom = {"nbi":dom["nbi"],"nbj":dom["nbj"]-hh[i], \
"istart":dom["istart"],"jstart":dom["jstart"]+hh[i]}
dom["nbj"] = hh[i]
- new_dom['nbland'] = int(np.nansum(land[new_dom["jstart"]:new_dom["jstart"]+new_dom["nbj"],new_dom["istart"]:new_dom["istart"]+new_dom["nbi"]]))
+ new_dom['nbland'] = int(np.nansum(land[new_dom["jstart"]:new_dom["jstart"]+new_dom["nbj"],\
+ new_dom["istart"]:new_dom["istart"]+new_dom["nbi"]]))
dom['nbland'] = dom['nbland']-new_dom['nbland']
partout.append(dom)
@@ -47,23 +50,24 @@ def halfpartition(partin, land) :
def fit_partition(partin, land):
partout = []
for i, dom in enumerate(partin):
- l = land[dom["jstart"]:dom["jstart"]+dom["nbj"],dom["istart"]:dom["istart"]+dom["nbi"]]
- l1 = np.sum(l, axis = 0)
-
- i0 = np.where(l1>0)[0][0]
- i1 = np.where(l1>0)[0][-1]
+ if dom["nbland"] > 0 :
+ l = land[dom["jstart"]:dom["jstart"]+dom["nbj"],dom["istart"]:dom["istart"]+dom["nbi"]]
+ l1 = np.sum(l, axis = 0)
+
+ i0 = np.where(l1>0)[0][0]
+ i1 = np.where(l1>0)[0][-1]
- l2 = np.ma.sum(l, axis = 1)
- j0 = np.where(l2>0)[0][0]
- j1 = np.where(l2>0)[0][-1]
+ l2 = np.ma.sum(l, axis = 1)
+ j0 = np.where(l2>0)[0][0]
+ j1 = np.where(l2>0)[0][-1]
- dom["jstart"] = j0 + dom["jstart"]
- dom["nbj"] = j1-j0+1
- dom["istart"] = i0 + dom["istart"]
- dom["nbi"] = i1-i0+1
- dom["nbland"] = int(np.nansum(land[dom["jstart"]:dom["jstart"]+dom["nbj"],dom["istart"]:dom["istart"]+dom["nbi"]]))
+ dom["jstart"] = j0 + dom["jstart"]
+ dom["nbj"] = j1-j0+1
+ dom["istart"] = i0 + dom["istart"]
+ dom["nbi"] = i1-i0+1
+ dom["nbland"] = int(np.nansum(land[dom["jstart"]:dom["jstart"]+dom["nbj"],dom["istart"]:dom["istart"]+dom["nbi"]]))
- partout.append(dom)
+ partout.append(dom)
return partout
#
#
@@ -91,10 +95,6 @@ def adjustpart(partin, land, nbcore) :
nbptproc=np.array([dom['nbi']*dom['nbj'] for dom in partin])
nbland = [dom["nbland"] for dom in partin]
- if np.min(nbptproc) < 3 :
- ERROR("The smallest domain in the partition is of a size less than 5. This will probably not work.")
- sys.exit()
-
return partin
#
# Partition domain in two steps : 1) halving, 2) adjusting by nb land points
@@ -362,6 +362,12 @@ class partition :
#
self.allistart = []
self.alljstart = []
+ #
+ if self.size != len(part) :
+ ERROR("There are too many processors for the size of the domain.")
+ ERROR(str(self.size)+" processors. But partition could only achieve "+str(len(part))+" domain with land points")
+ sys.exit()
+ #
for i in range(self.size) :
self.allistart.append(part[i]["istart"])
self.alljstart.append(part[i]["jstart"])
diff --git a/RPPtools.py b/RPPtools.py
index 31b2053f90ab34271f19d400fc61c7d507b90760..fb483b835e5d5afc97f1c020aa323f6d34f57208 100644
--- a/RPPtools.py
+++ b/RPPtools.py
@@ -18,7 +18,7 @@ import time
import configparser
config=configparser.ConfigParser()
config.read("run.def")
-saveweights=config.get("OverAll", "SaveWeights", fallback="true")
+saveweights=config.get("OverAll", "WeightFile", fallback=None)
EarthRadius=config.getfloat("OverAll", "EarthRadius", fallback=6370000.0)
#
# Logging in MPI : https://groups.google.com/forum/#!topic/mpi4py/SaNzc8bdj6U
@@ -75,7 +75,7 @@ def maxradius(cent, lon, lat) :
radius=[]
for lx, ly in zip(lon, lat) :
radius.append(np.sqrt((cent[0]-lx)**2+(cent[1]-ly)**2))
- return max(radius)
+ return np.max(np.array(radius))
#
# Simple routine to dump a field into a file.
#
@@ -138,7 +138,7 @@ class compweights :
#
# Only compute the weights if the file does not exist.
#
- if not os.path.exists(wfile) :
+ if wfile is None or not os.path.exists(wfile) :
INFO("Computing weights")
nbpts = len(modelgrid.centers)
i = 0
@@ -195,9 +195,9 @@ class compweights :
INFO(str(part.rank)+" Area Overlap error : "+str(min(overlaperr))+"-"+str(max(overlaperr))+\
" Nb of Hydro grid overlaping : "+str(min(self.hpts))+"-"+str(max(self.hpts)) )
#
- self.dumpweights(wfile, part, modelgrid, hydrogrid)
+ if wfile is not None :
+ self.dumpweights(wfile, part, modelgrid, hydrogrid)
#
- ##self.printtest([15,19], part, modelgrid)
#
else :
#
diff --git a/RoutingPreProc.py b/RoutingPreProc.py
index 3ebfcc81f6b680e9ece7b2e2c299ea26e24554d8..2ae15648e82c20162e2344b6151048add10fc654 100644
--- a/RoutingPreProc.py
+++ b/RoutingPreProc.py
@@ -21,7 +21,7 @@ nbasmax=config.getint("OverAll", "nbasmax")
numop=config.getint("OverAll", "numop", fallback=100)
OutGraphFile=config.get("OverAll","GraphFile")
DumpHydroSuper=config.getboolean("OverAll","DumpHydroSuper",fallback=False)
-wfile=config.get("OverAll","WeightFile",fallback="Weights.nc")
+wfile=config.get("OverAll","WeightFile",fallback=None)
#
import ModelGrid as MG
import Partition as PA
diff --git a/Tools/AddExtremes.py b/Tools/AddExtremes.py
new file mode 100644
index 0000000000000000000000000000000000000000..7ae4e4a4e9322f49667c07c351dad375027c2e60
--- /dev/null
+++ b/Tools/AddExtremes.py
@@ -0,0 +1,45 @@
+#
+#
+import numpy as np
+import os
+from netCDF4 import Dataset
+#
+filename="routing_MED.nc"
+#
+#################################################################################
+#
+def ncvarextrema(varin) :
+ expval = None
+ for (a,v) in varin.__dict__.items():
+ if ( a == "_FillValue" ) :
+ expval = v
+ elif (a == "missing_value" ) :
+ expval = v
+
+ if expval == None :
+ minval = np.min(varin[:,:])
+ maxval = np.max(varin[:,:])
+ else :
+ x = varin[:,:]
+ minval = np.min(x[np.where(x < expval)])
+ maxval = np.max(x[np.where(x < expval)])
+
+ return minval, maxval
+#
+#
+#
+def addextrema(varin, vmin, vmax) :
+ exec('varin.%s = "%s"' % ('min',vmin))
+ exec('varin.%s = "%s"' % ('max',vmax))
+#
+#
+#
+nf=Dataset(filename, 'a', format='NETCDF4_CLASSIC')
+#
+vars=["basins", "topoind_h", "topoind", "fac", "disto"]
+
+for v in vars :
+ vmin, vmax = ncvarextrema(nf.variables[v])
+ addextrema(nf.variables[v], vmin, vmax)
+
+nf.close()
diff --git a/Tools/AddFacDisto.py b/Tools/AddFacDisto.py
new file mode 100644
index 0000000000000000000000000000000000000000..22a557232a590e56043ef6221593d98ff2759703
--- /dev/null
+++ b/Tools/AddFacDisto.py
@@ -0,0 +1,193 @@
+#
+# Code to add flow accumulation and distance to ocean to the old routing.nc file
+#
+# fac is probably the accumulation of upstream grid boxes which flow through a given point.
+#
+import numpy as np
+import os, sys
+from netCDF4 import Dataset
+#
+rose=np.array([[0,0],[-1,0],[-1,+1],[0,+1],[+1,+1],[+1,0],[+1,-1],[0,-1],[-1,-1]])
+debug=False
+FillValue=1.0e+20
+#
+infile="routing.nc"
+outfile="routing_GLOB.nc"
+#
+# Copy variable
+def copyvar (innc, outnc, var, MinMax=False):
+ #
+ if ( debug ) : print("Working ", var)
+ shape=innc.variables[var].shape
+ #
+ # Scalar
+ #
+ if (len(shape) == 0 ) :
+ dtyp=innc.variables[var].datatype
+ outvar = outnc.createVariable(var,dtyp,())
+ outvar[0] = innc.variables[var][:]
+ #
+ # Vector
+ #
+ elif (len(shape) == 1 ) :
+ dtyp=innc.variables[var].datatype
+ ddim=innc.variables[var].dimensions
+ outvar = outnc.createVariable(var,dtyp,ddim)
+ outvar[:] = innc.variables[var][:]
+ #
+ # 2D matrix
+ #
+ elif (len(shape) == 2 ) :
+ dtyp=innc.variables[var].datatype
+ ddim=innc.variables[var].dimensions
+ outvar = outnc.createVariable(var,dtyp,ddim,fill_value=FillValue)
+ x = np.copy(innc.variables[var][:,:])
+ x[np.greater_equal(x, FillValue, where=~np.isnan(x))] = np.nan
+ outvar[:,:] = x
+ #
+ if MinMax :
+ outvar.min = np.nanmin(x)
+ outvar.max = np.nanmax(x)
+ #
+ else :
+ print('A 3D var is not yet foreseen in copyvar')
+ sys.exit(2)
+ #
+ # Transfer attributes
+ #
+ copyatt(innc.variables[var], outvar)
+ #
+ # Corrections to attributes
+ #
+ outvar.axis="YX"
+ outvar.coordinates="nav_lat, nav_lon"
+#
+# Function to copy attributes
+#
+#################################################################################
+
+# Function to copy attributes
+#
+#################################################################################
+#
+# Function which copied the attributes from one file netCDF file to another.
+#
+def copyatt(varin, varout):
+ for (a,v) in varin.__dict__.items():
+ if ( a != "_FillValue" ) :
+ if isinstance(v, str) :
+ vv=v.replace('\n','').replace('_',' ')
+ else :
+ vv=v
+ if ( debug ) : print("Attribute ", a, "--", vv)
+ exec('varout.%s = "%s"' % (a,vv))
+#
+# Function to copy global attributes
+#
+def copyglobatt(innc, outnc):
+ for (a,v) in innc.__dict__.items():
+ if ( a != "_FillValue" ) :
+ if isinstance(v, str) :
+ vv=v.replace('\n','').replace('_',' ')
+ else :
+ vv=v
+ if ( debug ) : print("Global Attribute ", a, "--", vv, type(vv))
+ setattr(outnc, a, vv)
+#
+#
+###################################################################################
+#
+#
+#
+###################################################################################
+#
+inf=Dataset(infile,'r')
+out=Dataset(outfile, 'w', format='NETCDF4_CLASSIC')
+#
+jdim,idim = inf.variables["nav_lon"].shape
+#
+# Creat output dimensions
+#
+xdim = out.createDimension('x', idim)
+ydim = out.createDimension('y', jdim)
+#
+#
+#
+lon=np.copy(inf.variables["nav_lon"][:,:])
+copyvar(inf, out, "nav_lon")
+lat=np.copy(inf.variables["nav_lat"][:,:])
+copyvar(inf, out, "nav_lat")
+#
+trip=np.copy(inf.variables["trip"][:,:])
+copyvar(inf, out, "trip")
+#
+basins=np.copy(inf.variables["basins"][:,:])
+copyvar(inf, out, "basins", MinMax=True)
+#
+topoind=np.copy(inf.variables["topoind"][:,:])
+copyvar(inf, out, "topoind", MinMax=True)
+#
+hdiff=np.copy(inf.variables["hdiff"][:,:])
+copyvar(inf, out, "hdiff", MinMax=True)
+#
+riverl=np.copy(inf.variables["riverl"][:,:])
+copyvar(inf, out, "riverl", MinMax=True)
+#
+orog=np.copy(inf.variables["orog"][:,:])
+copyvar(inf, out, "orog", MinMax=True)
+#
+varorog=np.copy(inf.variables["varorog"][:,:])
+copyvar(inf, out, "varorog", MinMax=True)
+#
+#
+# Compute the new variables
+#
+trip[np.greater_equal(trip, FillValue, where=~np.isnan(trip))] = np.nan
+fac=np.zeros((jdim,idim))
+disto=np.zeros((jdim,idim))
+disto[:,:]=np.nan
+for i in range(idim):
+ for j in range(jdim):
+ iflow=i
+ jflow=j
+ distotmp=0
+ while trip[jflow,iflow] > 0 and ~np.isnan(trip[jflow,iflow]) and trip[jflow,iflow] < 97:
+ # Add grid to flow accumulation
+ fac[jflow,iflow]=fac[jflow,iflow]+1
+ distotmp=distotmp+riverl[jflow,iflow]
+ # Move on
+ inci=rose[int(trip[jflow,iflow])][1]
+ incj=rose[int(trip[jflow,iflow])][0]
+ iflow=(iflow+inci)%idim
+ jflow=jflow+incj
+ #
+ # Save distance to ocean
+ #
+ if ~np.isnan(trip[j,i]) :
+ disto[j,i] = distotmp
+#
+# Add new fields to file
+#
+facnc = out.createVariable('fac',np.float,('y','x'),fill_value=FillValue)
+facnc.long_name="Flow accumulation"
+facnc.associate="nav_lat nav_lon"
+facnc.coordinates="nav_lat, nav_lon"
+facnc.units=""
+facnc.min=np.nanmin(fac)
+facnc.max=np.nanmax(fac)
+facnc[:,:] = fac[:,:]
+distonc = out.createVariable('disto',np.float,('y','x'),fill_value=FillValue)
+distonc.long_name="Distance to the ocean"
+distonc.associate="nav_lat nav_lon"
+distonc.coordinates="nav_lat, nav_lon"
+distonc.axis="YX"
+distonc.units="m"
+distonc.min=np.nanmin(disto)
+distonc.max=np.nanmax(disto)
+distonc[:,:] = disto[:,:]
+#
+# Close files
+#
+inf.close()
+out.close()
+
diff --git a/tests/Mallorca/BuildHTUs_Mallorca.pbs b/tests/Mallorca/BuildHTUs_Mallorca.pbs
index dacfbb50a98f6327ab0feff81535d652c881d00e..9cb1c5b62e0eb741bb5970da10bb96d808f92443 100644
--- a/tests/Mallorca/BuildHTUs_Mallorca.pbs
+++ b/tests/Mallorca/BuildHTUs_Mallorca.pbs
@@ -3,7 +3,7 @@
#PBS -N BuildHTUs_Mallorca
#
#PBS -j oe
-#PBS -l nodes=1:ppn=2
+#PBS -l nodes=1:ppn=32
#PBS -l walltime=4:00:00
#PBS -l mem=80gb
#PBS -l vmem=80gb
@@ -17,42 +17,57 @@ source ../../Environment
#
# Clean-up
#
-/bin/rm -f DocumentationInterface *.nc *.txt
+/bin/rm -f DocumentationInterface *.txt
#
-# Force RoutingPreProc to recompute the weights.
#
-/bin/rm -rf Weights
-#
-# 1 Proc
+# 2 Proc
#
-mpirun -n 1 python ../../RoutingPreProc.py
+/bin/rm -f run.def
+cp run_MEDCORDEX.def run.def
+mpirun -n 2 python ../../RoutingPreProc.py
if [ $? -gt 0 ] ; then
- echo "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
- echo "X Run on 1 Proc failed X"
- echo "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
+ echo "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
+ echo "X Run MECORDEX on 2 Proc failed X"
+ echo "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
exit
else
- echo "============================"
- echo "= Run on 1 Proc successful ="
- echo "============================"
- mv MEDCORDEX_test_graph.nc MEDCORDEX_test_graph_n1.nc
- mv MEDCORDEX_test_graph_HydroSuper.nc MEDCORDEX_test_graph_HydroSuper_n1.nc
+ echo "======================================"
+ echo "= Run MEDCORDEX on 2 Proc successful ="
+ echo "======================================"
+ ls -lt
fi
#
#
-# 2 Proc
+# 4 Proc
#
-mpirun -n 2 python ../../RoutingPreProc.py
+/bin/rm -f run.def *.txt
+cp run_E2OFD.def run.def
+mpirun -n 3 python ../../RoutingPreProc.py
+if [ $? -gt 0 ] ; then
+ echo "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
+ echo "X Run E2OFD on 3 Proc failed X"
+ echo "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
+ exit
+else
+ echo "=================================="
+ echo "= Run E2OFD on 3 Proc successful ="
+ echo "=================================="
+ ls -lt
+fi
+#
+# More points on 32 proc
+#
+/bin/rm -f run.def *.txt
+cp run_EuroSW.def run.def
+mpirun -n 32 python ../../RoutingPreProc.py
if [ $? -gt 0 ] ; then
- echo "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
- echo "X Run on 2 Proc failed X"
- echo "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
+ echo "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
+ echo "X Run EuroSW on 32 Proc failed X"
+ echo "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
else
- echo "============================"
- echo "= Run on 2 Proc successful ="
- echo "============================"
- mv MEDCORDEX_test_graph.nc MEDCORDEX_test_graph_n2.nc
- mv MEDCORDEX_test_graph_HydroSuper.nc MEDCORDEX_test_graph_HydroSuper_n2.nc
+ echo "===================================="
+ echo "= Run EuroSW on 32 Proc successful ="
+ echo "===================================="
+ ls -lt
fi
-ls -l
diff --git a/tests/Mallorca/run_E2OFD.def b/tests/Mallorca/run_E2OFD.def
new file mode 100644
index 0000000000000000000000000000000000000000..544545cad9fcef02f7c037be7f28da6d86f37efa
--- /dev/null
+++ b/tests/Mallorca/run_E2OFD.def
@@ -0,0 +1,27 @@
+[OverAll]
+#
+#
+EarthRadius = 6370000.
+#
+ModelGridFile = /bdd/MEDI/workspaces/polcher/NewRouting/MED_E2OFD_Contfrac.nc
+HydroFile = /bdd/MEDI/workspaces/polcher/NewRouting/routing_MED.nc
+# Mallorca
+WEST_EAST = 2.3, 3.75
+SOUTH_NORTH = 39.00, 40.1
+#
+# FORTRAN interface parameters
+#
+Documentation = true
+#
+# Configuration for the graph to be generated
+#
+nbasmax = 35
+#
+# Number of operation of simplification performed together
+#
+numop = 200
+#
+# Output
+#
+GraphFile = E2OFD_test_graph.nc
+#
diff --git a/tests/Mallorca/run_EuroSW.def b/tests/Mallorca/run_EuroSW.def
new file mode 100644
index 0000000000000000000000000000000000000000..d0186a0e682e47b296fe4ea7454eadb7debb49a3
--- /dev/null
+++ b/tests/Mallorca/run_EuroSW.def
@@ -0,0 +1,26 @@
+[OverAll]
+#
+#
+#
+EarthRadius = 6370000.
+#
+ModelGridFile = /bdd/MEDI/workspaces/polcher/NewRouting/geo_EuroSW.nc
+HydroFile = /bdd/MEDI/workspaces/polcher/NewRouting/routing_MED.nc
+# Mallorca
+WEST_EAST = 2.2, 3.6
+SOUTH_NORTH = 39.20, 40.10
+#
+WeightFile = EuroSW_Weights.nc
+#
+# FORTRAN interface parameters
+#
+Documentation = true
+#
+# Configuration for the graph to be generated
+#
+nbasmax = 35
+#
+# Output
+#
+GraphFile = EuroSW_test_graph.nc
+#
\ No newline at end of file
diff --git a/tests/Mallorca/run.def b/tests/Mallorca/run_MEDCORDEX.def
similarity index 80%
rename from tests/Mallorca/run.def
rename to tests/Mallorca/run_MEDCORDEX.def
index 981458b2511f0ee507e5cd69801ae568cb179186..97ce714e11aacef36d0a3ae939fdd19eea9af352 100644
--- a/tests/Mallorca/run.def
+++ b/tests/Mallorca/run_MEDCORDEX.def
@@ -4,10 +4,10 @@
EarthRadius = 6370000.
#
ModelGridFile = /bdd/MEDI/workspaces/polcher/NewRouting/geo_MEDCORDEX.nc
+HydroFile = /bdd/MEDI/workspaces/polcher/NewRouting/routing_MED.nc
# Mallorca
WEST_EAST = 2.3, 3.5
SOUTH_NORTH = 39.00, 40.1
-HydroFile = /bdd/MEDI/workspaces/polcher/NewRouting/routing_MED.nc
#
# FORTRAN interface parameters
#
@@ -25,8 +25,3 @@ numop = 200
#
GraphFile = MEDCORDEX_test_graph.nc
#
-# Diagnostics
-# You need to provide an interval in longitude and Latitude.
-#
-DiagLon = 3.9, 3.9
-DiagLat = 40.0, 40.0