diff --git a/examples/std-greenland/basemapfigs.py b/examples/std-greenland/basemapfigs.py
index 5baa370d89..ee7076169d 100755
--- a/examples/std-greenland/basemapfigs.py
+++ b/examples/std-greenland/basemapfigs.py
@@ -88,6 +88,34 @@
   myvmin     = 0.0
   myvmax     = 4.0
   ticklist   = [0, 1, 2, 3, 4]
+elif field == 'bmelt':
+  fill       = -2.0e+09
+  logscale   = True
+  contour100 = False
+  myvmin     = 0.9e-4
+  myvmax     = 1.1
+  ticklist   = [0.0001, 0.001, 0.01, 0.1, 1.0]
+elif field == 'tillwat':
+  fill       = -2.0e+09
+  logscale   = False
+  contour100 = False
+  myvmin     = 0.0
+  myvmax     = 2.0
+  ticklist   = [0.0, 0.5, 1.0, 1.5, 2.0]
+elif field == 'bwat':
+  fill       = -2.0e+09
+  logscale   = True
+  contour100 = False
+  myvmin     = 0.9e-4
+  myvmax     = 1.1
+  ticklist   = [0.0001, 0.001, 0.01, 0.1, 1.0]
+elif field == 'bwprel':
+  fill       = -2.0e+09
+  logscale   = False
+  contour100 = False
+  myvmin     = 0.0
+  myvmax     = 1.0
+  ticklist   = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]
 else:
   print 'invalid choice for FIELD option'
   sys.exit(3)
@@ -110,11 +138,28 @@
 width  = x.max() - x.min()
 height = y.max() - y.min()
 
-# load data and mask out ice-free areas
-myvar = np.squeeze(nc.variables[field][:])
+# load data
+if field == 'bwprel':
+    thkvar = np.squeeze(nc.variables['thk'][:])
+    myvar = np.squeeze(nc.variables['bwp'][:])
+    myvar = np.ma.array(myvar, mask=(thkvar == 0.0))
+    thkvar = np.ma.array(thkvar, mask=(thkvar == 0.0))
+    myvar = myvar / (910.0 * 9.81 * thkvar)
+else:
+    myvar = np.squeeze(nc.variables[field][:])
+
+# mask out ice free etc.
 if field == 'surfvelmag':
     myvar = myvar.transpose()
-myvar = np.ma.array(myvar, mask=(myvar == fill))
+    thkvar = np.squeeze(nc.variables['thk'][:]).transpose()
+    myvar = np.ma.array(myvar, mask=(thkvar == 0.0))
+elif (field == 'bmelt') | (field == 'bwat'):
+    myvar[myvar < myvmin] = myvmin
+    maskvar = np.squeeze(nc.variables['mask'][:])
+    myvar = np.ma.array(myvar, mask=(maskvar != 2))
+else:
+    maskvar = np.squeeze(nc.variables['mask'][:])
+    myvar = np.ma.array(myvar, mask=(maskvar != 2))
 
 m = Basemap(width=1.1*width,    # width in projection coordinates, in meters
             height=1.05*height,      # height
diff --git a/examples/std-greenland/hydro/.gitignore b/examples/std-greenland/hydro/.gitignore
new file mode 100644
index 0000000000..312e1da7bf
--- /dev/null
+++ b/examples/std-greenland/hydro/.gitignore
@@ -0,0 +1,3 @@
+basemapfigs.py
+Greenland_5km_v1.1.nc
+pism_Greenland_5km_v1.1.nc
diff --git a/examples/std-greenland/hydro/README.md b/examples/std-greenland/hydro/README.md
new file mode 100644
index 0000000000..8c8f2d7a67
--- /dev/null
+++ b/examples/std-greenland/hydro/README.md
@@ -0,0 +1,35 @@
+std-greenland/hydro/
+===========
+
+These subglacial hydrology runs are documented in the paper
+
+*  E. Bueler, W. Van Pelt (2014) _Mass-conserving subglacial hydrology in the_
+   _Parallel Ice Sheet Model_. Geoscientific Model Development Discussions
+   7 (4) pp. 4705–4775, doi:10.5194/gmdd-7-4705-2014
+
+These runs start with input data from files generated one level up (directory
+`examples/std-greenland/`), especially from script `gridseq.sh`.  Thus these
+runs are based on the well-documented runs in Chapter 1 of the PISM User's
+Manual, specifically those in the "grid sequencing" section.
+
+Do
+
+    $ ln -s ../pism_Greenland_5km_v1.1.nc
+    $ ln -s ../Greenland_5km_v1.1.nc
+
+Run a grid sequencing like in Chapter 1 of the User's Manual to get
+`g2km_gridseq.nc` (or similar).  Then do
+
+    $ ./run-decoupled.sh 5 g2km_gridseq.nc              # 5 year runs
+
+This run produces four files of importance, namely `routing-decoupled.nc`,
+`ex_routing-decoupled.nc`, `distributed-decoupled.nc`, `ex_distributed-decoupled.nc`.
+
+To generate figures from the paper do:
+
+    $ ln -s ../basemapfigs.py
+    $ ./allfigs.sh g2km_gridseq
+
+The second script `allfigs.sh` calls other figure-generating
+scripts: `basemapfigs.py`, `genGreenfig.sh`, `genscatfig.sh`,
+and `showPvsW.py`.  It also uses `mogrify` from the [Imagemagick](http://www.imagemagick.org/) tools.)
diff --git a/examples/std-greenland/hydro/allfigs.sh b/examples/std-greenland/hydro/allfigs.sh
new file mode 100755
index 0000000000..3831e9bdfc
--- /dev/null
+++ b/examples/std-greenland/hydro/allfigs.sh
@@ -0,0 +1,29 @@
+#!/bin/bash
+
+# generate all figures
+
+INITROOT=$1
+
+./genGreenfig.sh $INITROOT Greenland_5km_v1.1
+mv -v $INITROOT-velbase_mag.png $INITROOT-velbase-mag.png
+mv -v $INITROOT-velsurf_mag.png $INITROOT-velsurf-mag.png
+mv -v Greenland_5km_v1.1-surfvelmag.png Greenland-surfvelmag.png
+
+./basemapfigs.py routing-decoupled tillwat
+./basemapfigs.py routing-decoupled bwat
+
+./basemapfigs.py distributed-decoupled tillwat
+./basemapfigs.py distributed-decoupled bwat
+./basemapfigs.py distributed-decoupled bwprel
+
+./genscatfig.sh ex_distributed-decoupled.nc $INITROOT.png
+
+rm -rf listpng.txt
+ls *.png > listpng.txt
+
+for name in `cat listpng.txt`; do
+  echo "autocropping ${name} ..."
+  mogrify -trim +repage $name
+done
+rm -rf listpng.txt
+
diff --git a/examples/std-greenland/hydro/genGreenfig.sh b/examples/std-greenland/hydro/genGreenfig.sh
new file mode 100755
index 0000000000..a04f0fb563
--- /dev/null
+++ b/examples/std-greenland/hydro/genGreenfig.sh
@@ -0,0 +1,18 @@
+#!/bin/bash
+
+# these figures show "inputs" to the hydrology model, i.e. fields
+#   velsurf_mag, velbase_mag, bmelt
+# and surfvelmag from Greenland_5km_v1.1.nc
+# run as:
+#   $ ./genGreenfig.sh g2km-init Greenland_5km_v1.1
+# to generate figure files
+
+RESULTROOT=$1
+GREENROOT=$2
+
+for varname in bmelt velsurf_mag velbase_mag; do
+   ./basemapfigs.py $RESULTROOT $varname
+done
+
+./basemapfigs.py $GREENROOT surfvelmag
+
diff --git a/examples/std-greenland/hydro/genscatfig.sh b/examples/std-greenland/hydro/genscatfig.sh
new file mode 100755
index 0000000000..50a819e265
--- /dev/null
+++ b/examples/std-greenland/hydro/genscatfig.sh
@@ -0,0 +1,19 @@
+#!/bin/bash
+
+# these figures show velocity bins, with coloring by ice thickness,
+# where P(W) looks quite different in each bin
+# each bin is managably small for even 2km Greenland run, I think
+# run as:
+#   $ ./genscatfig.sh ex_distributed-decoupled.nc g2km.png
+# to generate figure files like
+#   bin*-g2km.png
+
+FILENAME=$1
+OUTROOT=$2
+
+./showPvsW.py -wmin 0.0 -wmax 0.15 -c thk -cmin 0 -cmax 2000 -s hydrovelbase_mag -smin 300  -smax 600 -o bin300-${OUTROOT} --colorbar $FILENAME
+./showPvsW.py -wmin 0.0 -wmax 0.15 -c thk -cmin 0 -cmax 2000 -s hydrovelbase_mag -smin 100  -smax 200  -o bin100-${OUTROOT}  $FILENAME
+./showPvsW.py -wmin 0.0 -wmax 0.15 -c thk -cmin 0 -cmax 2000 -s hydrovelbase_mag -smin 30   -smax 60  -o bin30-${OUTROOT}   $FILENAME
+./showPvsW.py -wmin 0.0 -wmax 0.15 -c thk -cmin 0 -cmax 2000 -s hydrovelbase_mag -smin 10   -smax 20   -o bin10-${OUTROOT}   $FILENAME
+./showPvsW.py -wmin 0.0 -wmax 0.15 -c thk -cmin 0 -cmax 2000 -s hydrovelbase_mag -smin 3    -smax 6   -o bin1-${OUTROOT}    $FILENAME
+
diff --git a/examples/std-greenland/hydro/run-decoupled.sh b/examples/std-greenland/hydro/run-decoupled.sh
new file mode 100755
index 0000000000..dbadd6e948
--- /dev/null
+++ b/examples/std-greenland/hydro/run-decoupled.sh
@@ -0,0 +1,51 @@
+#!/bin/bash
+
+# uses 2 files:
+#    pism_Greenland_5km_v1.1.nc  from examples/std-greenland/
+#    g2km_gridseq.nc             ditto; or similar
+# these files are documented in the PISM User's Manual, chapter 1
+#
+# do
+#   $ ./run-decoupled.sh 5 g2km_gridseq
+# for 5 year run
+
+set -e  # exit on error
+
+# check if env var PISM_DO was set (i.e. set PISM_DO=echo for a 'dry' run)
+if [ -n "${PISM_DO:+1}" ] ; then  # check if env var is already set
+  echo "#   PISM_DO = $PISM_DO  (already set)"
+else
+  PISM_DO=""
+fi
+
+DURATION=$1
+INNAME=$2
+
+MPIDO="mpiexec -n 6"
+
+CLIMATE="-surface given -surface_given_file pism_Greenland_5km_v1.1.nc"
+PHYS="-sia_e 3.0 -stress_balance ssa+sia -topg_to_phi 15.0,40.0,-300.0,700.0 -pseudo_plastic -pseudo_plastic_q 0.5 -till_effective_fraction_overburden 0.02 -tauc_slippery_grounding_lines"
+CALVING="-calving ocean_kill -ocean_kill_file pism_Greenland_5km_v1.1.nc"
+
+# run this to check for no shock: continue g2km_gridseq.nc run
+NAME=cont.nc
+cmd="$MPIDO pismr -i $INNAME -skip -skip_max 20 $CLIMATE $PHYS $CALVING -ts_file ts_$NAME -ts_times 0:yearly:$DURATION -y $DURATION -o $NAME"
+#$PISM_DO $cmd
+echo
+
+# suitable for -hydrology routing,distributed runs which are decoupled:
+EXVAR="mask,thk,topg,usurf,tillwat,bwat,hydrobmelt,bwatvel"
+EXVARDIST="${EXVAR},bwp,bwprel,hydrovelbase_mag"
+
+# -hydrology routing
+NAME=routing-decoupled.nc
+cmd="$MPIDO pismr -i $INNAME -no_mass -energy none -stress_balance none $CLIMATE -extra_file ex_$NAME -extra_times 0:monthly:$DURATION -extra_vars $EXVAR -ts_file ts_$NAME -ts_times 0:daily:$DURATION -hydrology routing -hydrology_bmelt_file $INNAME -report_mass_accounting -ys 0 -y $DURATION -max_dt 0.03 -o $NAME"
+$PISM_DO $cmd
+echo
+
+# -hydrology distributed
+NAME=distributed-decoupled.nc
+cmd="$MPIDO pismr -i $INNAME -no_mass -energy none -stress_balance none $CLIMATE -extra_file ex_$NAME -extra_times 0:monthly:$DURATION -extra_vars $EXVARDIST -ts_file ts_$NAME -ts_times 0:daily:$DURATION -hydrology distributed -hydrology_bmelt_file $INNAME -hydrology_velbase_mag_file $INNAME -report_mass_accounting -ys 0 -y $DURATION -max_dt 0.03 -o $NAME"
+$PISM_DO $cmd
+echo
+
diff --git a/examples/std-greenland/hydro/showPvsW.py b/examples/std-greenland/hydro/showPvsW.py
new file mode 100755
index 0000000000..958083507e
--- /dev/null
+++ b/examples/std-greenland/hydro/showPvsW.py
@@ -0,0 +1,195 @@
+#!/usr/bin/env python
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+import sys
+import argparse
+
+try:
+    from netCDF4 import Dataset as NC
+except:
+    print "netCDF4 is not installed!"
+    sys.exit(1)
+
+parser = argparse.ArgumentParser(description='show scatter plot P versus W from a PISM run')
+
+parser.add_argument('filename',
+                    help='file from which to get  P = bwprel  and  W = bwat')
+parser.add_argument('-o', default=None,
+                    help='output file for image, in a format matplotlib can write (e.g. .png, .pdf)')
+parser.add_argument('-d', type=int, default=-1,
+                    help='index of frame (default: last frame which is D=-1)')
+
+parser.add_argument('-c', default=None,
+                    help='name of variable to use to color points in scatter plot')
+parser.add_argument('-cmin', type=float, default=None,
+                    help='crop color values below at this value')
+parser.add_argument('-cmax', type=float, default=None,
+                    help='crop color values above at this value')
+parser.add_argument('--colorbar', action='store_true',
+                    help='whether to show colorbar() beside figure')
+
+parser.add_argument('-s', default=None,
+                    help='name of variable to use to select whether points appear in scatter plot')
+parser.add_argument('-smin', type=float, default=None,
+                    help='select minimum: if -c is used then below this value (of -s var) the points will not be plotted')
+parser.add_argument('-smax', type=float, default=None,
+                    help='select minimum: if -c is used then above this value (of -s var) the points will not be plotted')
+
+parser.add_argument('-wmin', type=float, default=None,
+                    help='lower limit on W axis')
+parser.add_argument('-wmax', type=float, default=None,
+                    help='upper limit on W axis')
+
+args = parser.parse_args()
+
+try:
+    nc = NC(args.filename, 'r')
+except:
+    print "ERROR: can't read from file ..."
+    sys.exit(1)
+
+print "  reading 'bwprel' field from file %s ..." % (args.filename)
+try:
+    bwprel = nc.variables["bwprel"]
+except:
+    print "ERROR: variable 'bwprel' not found ..."
+    sys.exit(2)
+
+print "  reading 'bwat' field from file %s ..." % (args.filename)
+try:
+    bwat = nc.variables["bwat"]
+except:
+    print "ERROR: variable 'bwat' not found ..."
+    sys.exit(3)
+
+if args.s != None:
+  print "  reading '%s' field, for point selection, from file %s ..." % (args.s, args.filename)
+  try:
+      ss = nc.variables[args.s]
+  except:
+      print "ERROR: variable '%s' not found ..." % (args.s)
+      sys.exit(5)
+  if args.smin != None:
+    print "    minimum value of '%s' for selection is %f" % (args.s, args.smin)
+  if args.smax != None:
+    print "    maximum value of '%s' for selection is %f" % (args.s, args.smax)
+
+if args.c != None:
+  print "  reading '%s' field, for point color, from file %s ..." % (args.c, args.filename)
+  try:
+      cc = nc.variables[args.c]
+  except:
+      print "ERROR: variable '%s' not found ..." % (args.c)
+      sys.exit(4)
+  if args.cmin != None:
+    print '    color minimum value is %f' % args.cmin
+  if args.cmax != None:
+    print '    color maximum value is %f' % args.cmax
+
+if args.d >= 0:
+   if np.shape(bwprel)[0] <= args.d:
+       print "ERROR: frame %d not available in variable bwprel" % args.d
+       sys.exit(4)
+   if np.shape(bwat)[0] <= args.d:
+       print "ERROR: frame %d not available in variable bwat" % args.d
+       sys.exit(5)
+   print "  using frame %d of %d frames" % (args.d, np.shape(bwat)[0])
+else:
+   args.d = -1
+   print "  reading last frame of %d frames" % (np.shape(bwat)[0])
+
+bwat   = np.asarray(np.squeeze(bwat[args.d,:,:])).flatten()
+bwprel = np.asarray(np.squeeze(bwprel[args.d,:,:])).flatten()
+
+bwprel[bwprel>1.0] = 1.0
+bwprel[bwprel<0.0] = 0.0
+
+if args.c != None:
+  ccc = np.asarray(np.squeeze(cc[args.d,:,:])).flatten()
+  if args.cmin != None:
+    ccc[ccc<args.cmin] = args.cmin
+  if args.cmax != None:
+    ccc[ccc>args.cmax] = args.cmax
+
+if args.s != None:
+  sss = np.asarray(np.squeeze(ss[args.d,:,:])).flatten()
+  if args.smin != None:
+    bwat = bwat[sss>=args.smin]
+    bwprel = bwprel[sss>=args.smin]
+    if args.c != None:
+      ccc = ccc[sss>=args.smin]
+    sss = sss[sss>=args.smin]
+  if args.smax != None:
+    bwat = bwat[sss<=args.smax]
+    bwprel = bwprel[sss<=args.smax]
+    if args.c != None:
+      ccc = ccc[sss<=args.smax]
+    sss = sss[sss<=args.smax]
+
+nc.close()
+
+# to reduce file size, remove zero water points
+if args.c != None:
+  ccc = ccc[bwat>0]
+bwprel = bwprel[bwat>0.0]
+bwat = bwat[bwat>0.0]
+
+# to reduce file size, remove zero pressure points
+if args.c != None:
+  ccc = ccc[bwprel>0]
+bwat = bwat[bwprel>0.0]
+bwprel = bwprel[bwprel>0.0]
+
+print "  scatter plot has %d points ...\n" % len(bwat)
+
+# W axis limits
+if args.wmin == None:
+  wwmin = min(bwat)
+else:
+  wwmin = args.wmin
+if args.wmax == None:
+  wwmax = max(bwat)
+else:
+  wwmax = args.wmax
+
+# color axis limits
+if args.c != None:
+  if args.cmin == None:
+    ccmin = min(ccc)
+  else:
+    ccmin = args.cmin
+  if args.cmax == None:
+    ccmax = max(ccc)
+  else:
+    ccmax = args.cmax
+
+import matplotlib
+if args.colorbar:
+  matplotlib.rcParams['figure.figsize'] = (3.8,3.0)
+else:
+  matplotlib.rcParams['figure.figsize'] = (3.0,3.0)
+
+plt.figure(1)
+if args.c != None:
+  plt.scatter(bwat,bwprel,c=ccc,vmin=ccmin,vmax=ccmax,linewidth=0.0,cmap='hsv')
+else:
+  plt.scatter(bwat,bwprel,c='k')
+plt.axis('tight')
+plt.gca().set_xlim((wwmin,wwmax))
+plt.gca().set_ylim((-0.02,1.05))
+plt.gca().set_xticks((0.0,0.05,0.10,0.15))
+plt.gca().set_xticklabels(('0','.05','.10','.15'))
+plt.text(0.07, 0.25, '%d' % int(args.smin) + r'$ < |\mathbf{v}_b| < $' + '%d' % int(args.smax))
+if args.colorbar:
+  plt.colorbar()
+plt.xlabel(r'$W$  (m)')
+plt.ylabel(r'$P/P_0$')
+
+if args.o == None:
+  plt.show()
+else:
+  print "  saving scatter plot in %s ...\n" % args.o
+  plt.savefig(args.o, dpi=300, bbox_inches='tight')
+