Final version of L1L2 for TC paper by Robinson et al.

This commit streamlines the 3D velocity subroutine for L1L1 and makes
a small but important change in the algorithm.

The change is to evaluate the membrane stresses in tau_xz and tau_yz
at layer midpoints instead of lower layer boundaries, consistent
with where the SIA terms are evaluated.
I found that a small vertical offset can disrupt the balance between these
two terms and make the L1L2 solver unstable for the slab problem at resolutions
finer than ~200 m.  With the fix, the L1L2 stability curve is parallel
to the SIA stability curve at fine resolution, with L1L2 being slightly
more stable than SIA.

I also removed the alternate L1L2 discretization methods that were added
in a recent commit.  The alternate strategies evaluated some terms
at edges instead of vertices, and did not change the results significantly.
With this commit, all terms are evaluated at either cell centers or vertices.

This is the code version used for the slab tests shown for CISM
in Section 3 of the paper by Robinson, Goldberg & Lipscomb (2021, TC).
For these tests, I temporarily changed the energy conservation criterion
in glissade_therm.F90 to avoid false non-conservation alarms when
running at very fine resolution.  See the comments in that module.

This commit is answer-changing for L1L2, in a good way.

libglissade/glissade_therm.F90

@@ -1122,9 +1122,10 @@ subroutine glissade_therm_driver(whichtemp,                         &
              if (abs((efinal-einit-delta_e)/dttem) > 1.0d-7) then
              ! WHL: For stability tests with a very short time step (e.g., < 1.d-6 year),
              !      the energy-conservation error can be triggered by machine roundoff.
-             !      For these tests, I uncommented the line below, which compares the
-             !      error to the total amount of energy.  The latter criterion is less likely
-             !      to give false positives, but might be more likely to give false negatives.
+             !      For the tests in Robinson et al. (2021), I replaced the line above
+             !      with the line below, which compares the error to the total energy.
+             !      The latter criterion is less likely to give false positives,
+             !       but might be more likely to give false negatives.
 !!             if (abs((efinal-einit-delta_e)/(efinal)) > 1.0d-8) then
 
                 if (verbose_column) then