I'm not sure what you are doing wrong with your own interpolation/reprojecting, but it can be done even simpler.
The contourf
accepts a latlon
keyword which, when true, accepts lat/lon inputs and converts it to the map projection automatically. So:
datain = msk_array
fig = plt.figure(figsize=(12,5))
map = Basemap(projection='robin',resolution='c',lat_0=0,lon_0=0)
ny, nx = datain.shape
xin = np.linspace(map.xmin,map.xmax,nx) # nx is the number of x points on the grid
yin = np.linspace(map.ymin,map.ymax,ny) # ny in the number of y points on the grid
lons = np.arange(-180,180,0.25) #from raster.GetGeoTransform()
lats = np.arange(90,-90,-0.25)
lons, lats = np.meshgrid(lons,lats)
xx, yy = m(lons,lats)
levels = [-1000,-800,-600,-400,-200,0,200,400,600,800,1000]
cntr = map.contourf(xx, yy,datain, levels,cmap=cm.RdBu)
cbar = map.colorbar(cntr,location='bottom',pad='15%')
# Add some more info to the map
cstl = map.drawcoastlines(linewidth=.5)
meri = map.drawmeridians(np.arange(0,360,60), linewidth=.2, labels=[1,0,0,1], labelstyle='+/-', color='grey' )
para = map.drawparallels(np.arange(-90,90,30), linewidth=.2, labels=[1,0,0,1], labelstyle='+/-', color='grey')
boun = map.drawmapboundary(linewidth=0.5, color='grey')
Note that i changed the lats
definition in order to remove the flipping of your input raster, just a personal preference.