diff --git a/plotting_script.py b/plotting_script.py
new file mode 100644
index 0000000..2105866
--- /dev/null
+++ b/plotting_script.py
@@ -0,0 +1,20 @@
+import xarray as xr
+import matplotlib.pyplot as plt
+import numpy as np
+import cartopy.crs as ccrs
+from cartopy.util import add_cyclic_point
+
+import sys
+sys.path.insert(0, '/scratch/mp586/git_repos/ACDC/tools')
+from ds_plots import *
+from netCDF4 import Dataset
+
+# load dataset 1
+nc = Dataset('/scratch/mp586/git_repos/ACDC/Data/monthlyclimo_'+input('Which experiment ? ')+'.nc')
+lats = nc.variables['lat'][:]
+lons = nc.variables['lon'][:]
+
+P = xr.DataArray(nc.variables['precipitation'][:])*86400
+
+mmap(np.asarray(P.mean('dim_0')),lats,lons,show = 'Yes')
+
diff --git a/tools/ds_plots.py b/tools/ds_plots.py
index 3f1ab99..bb94667 100644
--- a/tools/ds_plots.py
+++ b/tools/ds_plots.py
@@ -110,7 +110,8 @@ def hmap(mapdata,lat,cmap=None,clim=None,title=None,
 def mmap(mapdata,lat,lon,title=None,cb_ttl=None,cmap=None,clim=None,
               filepath=None,show=None,sigmask=None,p=None,
               do_zonal=None,newfig=True,
-              outside_color = None, outside_val = None):
+              outside_color = None, outside_val=None,
+              use_pcolor=False):
     # Make a map (robinson projection) of mapdata over lat/lon
     # Inputs
     #   mapdata : 2D field to be plotted
@@ -122,7 +123,14 @@ def mmap(mapdata,lat,lon,title=None,cb_ttl=None,cmap=None,clim=None,
     #   cb_ttl  : colorbar title
     #   cmap    : special colormap
     #   clim    : colorbar limits, as a 2 element array
+    #   use_pcolor: default is to use contourf. Set True for pcolormesh
     
+    # mp added a few lines to convert potential xarray inputs to numpy arrays
+    # so that add_cyclic_point works
+    lon = np.asarray(lon)
+    lat = np.asarray(lat)
+    mapdata = np.asarray(mapdata)
+
     if newfig:
         if do_zonal==1:
             # make 2 subplots, one for map, one for zonal average next to it
@@ -140,7 +148,12 @@ def mmap(mapdata,lat,lon,title=None,cb_ttl=None,cmap=None,clim=None,
     cyclic_data, cyclic_lons = add_cyclic_point(mapdata,coord=lon)
 
     # Make the plot
-    cs = plt.pcolormesh(cyclic_lons,lat,cyclic_data,transform=ccrs.PlateCarree())
+    # Default: contourf
+    if not use_pcolor:
+        cs = plt.contourf(cyclic_lons,lat,cyclic_data,transform=ccrs.PlateCarree())
+    else:
+        cs = plt.pcolormesh(cyclic_lons,lat,cyclic_data,transform=ccrs.PlateCarree())
+
     
     # set colormap
     if cmap: