Committing the BlatterStressBalance code.

This version is isothermal (so far) and does not compute the
volumetric strain heating (yet).

Also, we still use the no-slip boundary condition at the base
(although the code should work without it).

To implement this, I had to change some things elsewhere in PISM:

- the IceFlowLaw::effective_viscosity() interface (takes second
  invariant "gamma" as an argument; this way we can use the same
  effective_voscosity() code in 2D and 3D cases

- moved the computation of the basal frictional heating from SSA to
  ShallowStressBalance; this way we can re-use it in
  BlatterStressBalance.

- renamed "D2" -> "strain_heating_contribution" (for clarity).

To try it out, run

./ssa_testj -Mx 120 -My 120 -snes_monitor
mpiexec -n 4 ./blatter_test -i ssa_test_j.nc -snes_monitor -ksp_monitor -blatter_Mz 10 -Lz 800 -o_order zyx

This is *without* the multigrid method, which also works. Try adding
"-pc_type mg -pc_mg_levels 2".
one-letter change ... hmmm, perhaps the original code was written on an OS with a case-insensitive file system ... must cause lots of Jobs for people on other OSes ...
Moved the 'extern "C"' before the class declaration.

This should fix the compile-time error on Ubuntu.
Got rid of compiler warnings.

(gcc with "-pedantic" uses the C89 (C90?) standard, and this code uses some C99 features.)
Renamed Sigma -> strain_heating throughout.
refactoring so that K is out of flux expressions and in its own routine, and adding power beta; does not compile
Renamed Sigma -> strain_heating.
Stronger emphasis on the fact that we require doxygen and graphviz to re-build source code browsers.
Renamed Sigma -> strain_heating throughout.
(cherry-picked from e4e9ee32188ab64a869196b9d9186e9160b001a4)
This form builds and runs.  Only tested with beta=2 so far.  Results from Test P are slightly different in value and slower (perhaps a lot slower).  I think the change that happened in the beta = 2 case is the change in the max(K W) code in PISMRoutingHydrology::adaptive_for_W_evolution().
Put computation of max(KW) in right spot.  Avoid unnecessary computation in beta=2 case.  Recovers full speed and exact regression on Test P.
Clarify units of k.  Switch Nbreen example over to Hewitt et al (2012) Table 1 parameter values for alpha, beta, k, W_r = h_r.  Looks good ... though its not clear what that means.
minor clarification of choices relative to literature
minor: edit of comment to match notation
I had gotten confused over the units of, and the size of, the hydraulic conductivity constant.  As a result I was working with really low-conductivity hydrology which implied unreasonably high W values (water thicknesses) and very fast run times.  With this commit the code (but not all the test cases) fully match the Schoof et al [2012] and Hewitt et al [2012] conventions.  The examples/nbreen/run.sh script and its override parameters in examples/nbreen/nbreen_config.cdl also match the Schoof/Hewitt work.  The model results for Nbreen then seem very reasonable relative the Schoof/Hewitt work.  (The results are more expensive computationally, but also more interesting physically!  The documentation in the private hydrolakes/ repo has not been fully updated to match yet.)  The Test P example cannot be used until it is inspected and updated; I think all we need to do is shorten the run.
Added a Python PDD modeling tool

PyPDD is a small Python PDD model that I wrote in the last days. For short documentation see:

https://github.com/jsegu/pypdd

I believe it does what PISM should do after issue #133 is fixed, except that it handles only monthly data, does not (yet) uses the Calov & Greve integrand, and integration is very simple (numpy.sum).
Added a PDD model test (details bel w)

The test runs:

    pypdd.py -b -o pypdd_smb.nc
    pisms -Lx 1 -Ly 1 -Mx 21 -My 21 -y 0 -o i.nc
    pclimate -i i.nc -o pism_smb.nc -times 0:1:1 -atmosphere given
      -atmosphere_given_file atm.nc -atmosphere_given_period 1
      -surface pdd -pdd_annualize -pdd_std_dev 0.01
    nccmp.py -v climatic_mass_balance pypdd_smb.nc pism_smb.nc

This produces four files:

* atm.nc, containing a fake climate generated by pypdd.py;
* i.nc, containing input for pclimate;
* pypdd_smb.nc, with the resulting mass-balance from pypdd.py;
* pism_smb.nc, with output from pclimate.

The fake climate in atm.nc is defined by:

    temp(t, x, y) = -10 * (y + cos(2*pi*t))
    prec(t, x, y) = x * (sign(x) - cos(2*pi*t))

With x, y in km and t in yr. You can think of it as a squared geographical domain where:

* temperature amplitude is the same (10°C) everywhere
* temperature always peaks in summer
* mean temperature increases from North (-10°C) to South (10°C)
* min precipitation is zero everywhere
* mean precipitation is zero in the center line of the domain, and increases to 1 m/yr in the western and eastern edges.
* precipitation peaks in winter in the western (maritime) domain half and in summer in the eastern (continental) domain half.

Expectedly, smb.nc and pism_smb.nc differ. Notably, pism_smb.nc is y-symetric, which reflects PISM's insensitivity to precipitation seasonality. This is most obvious in the accumulation field and consistent with issue #133.
Adds Dukowicz, Price, and Lipscomb (Consistent approximations...)
Fixes the PDD code. See #133. Needs testing.

The new code keeps track of the snow depth, but re-sets it to zero
once a year (on pdd_balance_year_start_day, October 1 by default).

For restartability the new field (snow_depth) is saved to output
files. It is not required for re-starting, though: if it is not found
in an input file, snow depth is set to zero.

Both -pdd_annualize and "regular" scheme does what I think they should
do. Using the -pdd_annualize option trades sub-annual ice-dynamics
feedback for a faster computation of the mass balance. In the absence
of such feedback the two schemes are equivalent (assuming that the run
started on October 1).

Also: the new code ignores negative precipitation, which we allowed
for historical reasons (as far as I can tell).
Disabled the PDD regression test.

It fails, but for the wrong reason (i.e. I don't think it helps
us find flaws in the new PDD code).