Implemented an efficient depression-filling algorithm

This commit adds an efficient algorithm for filling depressions in the head field
when running the flux-routing hydrology scheme.

The old scheme was very slow to converge on large grids such as the 4km
Northern Hemisphere grid.

The new scheme is based on the algorithm of Planchon and Darboux (2001).
The basic idea is:

* Initially, set phi = phi_in on the boundary, and set phi = a large number elsewhere
  (where phi is the head field).

* Loop through the domain.  For each cell c, with value phi(c) not yet fixed to a known value,
  compute phi_min8(n), the current minimum of phi in the 8 neighbor cells.

  - If phi_in(c) > phi_min8(n) + eps, then set phi(c) = phi_in(c) and mark that cell as having a known value,
    since phi(c) cannot go any lower.
  - If phi_in(c) < phi_min8(n) + eps, but phi(c) > phi_min8(c) + eps, set phi(c) = phi_min8(n) + eps.
    Do not mark the cell as having a known value, because it might be lowered further.

* Continue until no further lowering of phi is possible.  At that point, phi = phi_out.

Here, eps is a small number greater than zero, which ensures that there are no flat surfaces
when the depression-filling is done.  Thus, it is no longer necessary to call fix_flats.

On the 4km N.H. grid, the number of depression-fill iterations is reduced from several hundred
per time step to ~10.

libglide/glide_types.F90

@@ -1945,7 +1945,7 @@ module glide_types
                                                  !> Where u > umax, let u = umax when evaluating beta(u)
 
      ! Note: A basal process model is not currently supported, but a specified mintauf can be passed to subroutine calcbeta
-     !       to simulate a plastic bed..
+     !       to simulate a plastic bed.
      real(dp),dimension(:,:)  ,pointer :: mintauf => null() ! Bed strength (yield stress) calculated with basal process model
 
   end type glide_basal_physics

libglissade/glissade.F90

@@ -737,10 +737,6 @@ subroutine glissade_initialise(model, evolve_ice)
 
     if (model%climate%overwrite_acab_value /= 0 .and. model%options%is_restart == RESTART_FALSE) then
 
-       !WHL - debug
-       if (main_task) print*, 'overwrite_acab value (m/yr):', &
-            model%climate%overwrite_acab_value * scyr*thk0/tim0
-
        call glissade_overwrite_acab_mask(model%options%overwrite_acab,          &
                                          model%climate%acab,                    &
                                          model%geometry%thck,                   &
@@ -1987,6 +1983,11 @@ subroutine glissade_thermal_solve(model, dt)
           endwhere
        endif
 
+       !WHL - debug - Set mask = 0 where thck = 0 for dome test
+!       where (model%geometry%thck == 0)
+!          bwat_mask = 0
+!       endwhere
+
        call parallel_halo(bwat_mask, parallel)
 
        ! Compute bwat based on a steady-state flux routing scheme