Some clean-up in ModelGrid.py and tests that it works with decomposition in small partitions.

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.
         #

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

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
 #