This came from edrun/ in https://github.com/pism/green-hydro. I am co…

…mmitted to maintaining it here for stable0.7 and onward.

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

examples/std-greenland/hydro/.gitignore

@@ -0,0 +1,3 @@
+basemapfigs.py
+Greenland_5km_v1.1.nc
+pism_Greenland_5km_v1.1.nc

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

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
+

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
+

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
+

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
+