diff --git a/ModelGrid.py b/ModelGrid.py index 24e4a358728e05e3e30ea9464084116e666f812e..7661a43ff69ba8ea0ade4f465cb50999a3f04563 100644 --- a/ModelGrid.py +++ b/ModelGrid.py @@ -93,7 +93,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 # @@ -137,8 +137,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] @@ -167,14 +171,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 # # # @@ -209,20 +210,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. # @@ -231,15 +228,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,\ @@ -399,7 +387,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. # diff --git a/Partition.py b/Partition.py index 8038e5ae9aafca5e4e2452c85629b7c742eff72c..93165133451ee5dce7e52033e2a744c98f409f49 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) @@ -91,10 +94,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 diff --git a/tests/Mallorca/run_EuroSW.def b/tests/Mallorca/run_EuroSW.def index d0186a0e682e47b296fe4ea7454eadb7debb49a3..2fba5f8f6b15f8a2e090529162caf636b4829277 100644 --- a/tests/Mallorca/run_EuroSW.def +++ b/tests/Mallorca/run_EuroSW.def @@ -10,7 +10,7 @@ HydroFile = /bdd/MEDI/workspaces/polcher/NewRouting/routing_MED.nc WEST_EAST = 2.2, 3.6 SOUTH_NORTH = 39.20, 40.10 # -WeightFile = EuroSW_Weights.nc +##WeightFile = EuroSW_Weights.nc # # FORTRAN interface parameters #