Commit 728fde56 authored by Anthony's avatar Anthony
Browse files

Solving conflict

parents d8071e83 3aaa4036
......@@ -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
......
......@@ -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) :
......
......@@ -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:
......
......@@ -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))
......
......@@ -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"])
......
......@@ -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 :
#
......
......@@ -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
......
#
#
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()
#
# 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
#
#######################################################