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
Stronger emphasis on the fact that we require doxygen and graphviz to re-build source code browsers.
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.
Merge branch 'dev' into hydro-beta
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.
Added a PDD model test (details below)

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.
Merge branch 'dev' of github.com:pism/pism into dev
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).
Adds warnings about missing tools required to re-build docs.
enough of a start to break things, but not fix things; so it is a feature branch
Removes the -pdd_annualize option altogether.

Now that the PDD code is fixed (I hope) we don't need the "annualized"
version of the PDD. "Normal" version produces equivalent annual
average mass balance and does not interfere with time-stepping and
periodic forcing.