Add an output to check the area calculation

RPPtools.py

@@ -36,3 +36,43 @@ def boxit(cent, dlon, dlat, polyres) :
     # Close
     boxll.append(boxll[0])
     return boxll
+#
+# Simple routine to dump a field into a file.
+#
+def dumpfield(x, lon, lat, filename, varname) :
+    #
+    from netCDF4 import Dataset
+    NCFillValue=1.0e20
+    vtyp=np.float64
+    #
+    print("Dumping overlap into file :", filename)
+    i=np.nonzero(np.mean(x,axis=0))
+    j=np.nonzero(np.mean(x,axis=1))
+    ilen=np.max(i)-np.min(i)+1
+    jlen=np.max(j)-np.min(j)+1
+    print("Overlap")
+    outnf=Dataset(filename, 'w', format='NETCDF4_CLASSIC')
+    # Dimensions
+    outnf.createDimension('lon', ilen)
+    outnf.createDimension('lat', jlen)
+    lonnc = outnf.createVariable("lon", vtyp, ('lat','lon'), fill_value=NCFillValue)
+    lonnc.units="Longitude of grid box centre degrees east"
+    lonnc.standard_name="grid longitude"
+    lonnc.title="Longitude"
+    lonnc.axis="X"
+    lonnc[:,:] = lon[np.min(j):np.max(j)+1,np.min(i):np.max(i)+1]
+    latnc = outnf.createVariable("lat", vtyp, ('lat','lon'), fill_value=NCFillValue)
+    latnc.units="Latitude of grid box centre degrees north"
+    latnc.standard_name="grid latitude"
+    latnc.title="Latitude"
+    latnc.axis="Y"
+    latnc[:,:] = lat[np.min(j):np.max(j)+1,np.min(i):np.max(i)+1]
+    #
+    xnc = outnf.createVariable(varname, vtyp, ('lat','lon'), fill_value=NCFillValue)
+    xnc[:,:] = x[np.min(j):np.max(j)+1,np.min(i):np.max(i)+1]
+    outnf.close()
+    #
+    return
+#
+#
+#