Debugged the enthalpy solver and added new enthalpy diffusivity options

In consultation with Hari Rajaram and others at CU, I made several changes to ensure
that enthalpy is conserved and to prevent oscillations in thickness, velocity and temperature.
These changes apply to config option temperature = TEMP_ENTHALPY = 3.

Answers are unchanged for the standard option (temperature = 1), as confirmed by Mac LIVV.

I modified glissade.F90 to transport enthalpy instead of transporting
temperature and waterfrac separately.  Specific enthalpy is given by

        enth = (1 - phi_w) * rhoi * ci * T
                  + phi_w  * rhow * (L + ci * Tpmp)

where phi_w is the water fraction and Tpmp is the pressure melting point temperature.
Since this expression is slightly nonlinear in phi_w and T (it contains a term proportional
to phi_w*T), energy is not quite conserved if temperature and enthalpy are passed separately.
Passing enthalpy requires a setup call to temp2enth and a follow-up call to enth2temp.

In glissade_transport.F90, I removed some code associated with transporting waterfrac.
Now by default, only one tracer is transported (either temp or enthalpy).

I made several changes in glissade_temp.F90 and glissade_enthalpy.F90:

(1) I modified enth2temp to make sure waterfrac = 0 whenever T < Tpmp.

(2) I added code (based on what was used for standard temperature option) to check for energy
    conservation in the enthalpy scheme.  I verified that energy is conserved in the EISMiNT-2 test
    with the enthalpy solver.
    (Specifically, the total energy input to the ice in a time step, given by the sum of the internal
    dissipation and the conductive fluxes at the interfaces, is equal to the final energy
    minus the initial energy.)

(3) As part of this calculation, I made sure that the conductive flux to the bed is computed correctly,
    as required for the basal melt calculation.

(4) I added a thickness term and corrected the time units in the basal melt calculation;
    previously the computed basal melt was off by several orders of magnitude.

(5) I changed the way enthalpy diffusivity is computed at layer interfaces.  This was
    the change required to remove oscillations in thickness, temperature, velocity, and
    related variables.  These oscillations were occuring with a period of several thousand
    years, most pronounced near the ice margin, with a magnitude of ~10 m in H, ~2 m/yr in v,
    and ~1 degree in T.  They were caused by sudden large changes in diffusivity when
    a layer changed from cold to temperate or vice versa.

    My solution is as follows:

    Let d(enth)/dz = the gradient of enthalpy
    We can write
         d(enth)/dz = d(enth_T)/dz + d(enth_w)/dz
    where
         enth_T = (1-phi_w) * rhoi*ci*T,
         enth_w =    phi_w  * rhow*(L + ci*Tpmp).
    Now let f = d(enth_T)/dz / d(enth)/dz
         (with the constraint 0 <= f <= 1).
    For cold ice, f = 1 and alpha = alphai = ki/(rhoi*ci)
    For temperate ice, f ~ 0 and alpha << alphai.
    At the interface between cold and temperate ice:
     f ~ 0 if the temperate ice has large phi_w and the cold ice has T close to Tpmp;
     f ~ 1 if the temperate ice has close to zero phi_w and the cold ice has T << Tpmp.
    There are two options for averaging diffusivity between temperature nodes based on f:
    (1) arithmetic mean:   alpha = f*alphai + (1-f)*alpha0
    (2) harmonic mean:     alpha = 1 / (f/alphai + (1-f)/alpha0)
    Both means converge to the desired values for f = 0 and f = 1, but for
     intermediate values of f, the arithmetic mean gives a greater diffusivity.
    Thus the harmonic mean allows large temperature gradients to develop between cold and
     temperate ice, relative to the arithmetic mean.  I am not sure which average
     gives a better fit to observations; this remains to be studied.

For improved diagnostics, I added enthalpy as an I/O variable in glide_vars.def.  It belongs
to the glide_temper derived type.

I amoved the 3D dissipation variable, dissip, from glide_tempwk to glide_temper, since it is a useful
diagnostic and not just a work variable.  The vertical staggering of dissip depends on the dycore
(Glide v. Glissade), like temp and flwa, so I added code in glimmer_ncio.F90 to make the switch
between dissip and dissipstag if desired.

I also incorporated a couple of changes made by Ryan Haagenson (CU):
(1) He changed the time stepping from Crank-Nicolson to backward Euler for improved stability.
(2) He removed the dependence of flwa on waterfrac.  This will be restored later, after more testing.

And I added a number of debug print statements that will be cleaned up in a future commit.
Restored the dependence of flwa on waterfrac for the enthalpy solver

Now, flwa can increase by up to a factor of 2.8 with waterfrac = 0.01.
This dependence had previously been commented out by Ryan H. for testing.

I verified that with this dependence included, the EISMINT-2 ice sheet
converges toward a faster, thinner steady state, as expected.

The evolution of thickness, temperature, etc. is still smooth, as with
the previous commit.

No answer changes for standard config settings (temperature = 1).
Created a new glissade_therm module that combines glissade_temp and glissade_enthalpy.

The files glissade_temp.F90 and glissade_enthalpy.F90 are still in the repo,
but later will be deleted.

Answers are unchanged for the standard dome evolution test (temperature = 1).
Answers changed at roundoff level for the EISMINT-2 enthalpy test (temperature = 3).

After combining glissade_temp and glissade_enthalpy into glissade_therm,
I did a lot of cleanup:
(1) I removed the old tempwk variables and some other unused variables.
(2) I removed 'model' from the argument list for all subroutines
    other than init and driver.
(3) I made the enthalpy code more closely integrated with the standard temperature
    code, so as to reduce enthalpy-specific code to a minimum.
(4) I switched to using Glissade ice and floating masks.

More cleanup still to come.
More cleanup of glissade_therm.F90

These are further changes toward the goal of having a Glissade temperature module
that is easy to understand and can readily be translated to MPAS.

Answers are unchanged from the previous commit on Nov. 11 (based on the evolving dome
with standard temperature and EISMINT-2 with the enthalpy scheme).

Some changes:
(1) I renamed some subroutines.
(2) I cleaned up the pressure melting point subroutines.
(3) I removed many uses of the model derived type syntax from inside glissade_velo_driver.
    I replaced these with pointers to components of the derived type.  This makes the code
    easier to read and later translate to MPAS.

Note: I had trouble getting temp to point correctly to model%temper%temp (because of the range
starting with 0 instead of 1?) and had to implement a workaround.

There are still some model derived types in the code, but I'm hoping to eliminate these next
by moving some subroutines elsewhere (e.g., the dissipation and basal friction subroutines).
Moved heat dissipation and basal friction calculations; replaced pointers with subroutine arguments

Answers for standard problems are unchanged (verified by Mac LIVV for evolving dome).

Previously, the interior heat dissipation was computed in glissade_therm for both Glam and Glissade.
Basal friction was computed in glissade_therm for Glam (but was already part of the velocity solution
for Glissade).

With this commit, the call to the heat dissipation subroutine moves from glissade_therm.F90 to glissade.F90,
just after the velocity solve.  This means that for exact restart, the 3D dissip field is now required,
but the 3D fields of tau_eff and efvs are not.  I changed the restart variable list accordingly
(in glide_setup.F90) and verified that the model restarts exactly for dome and EISMINT-2 cases.
In glide_vars.def, dissip and dissipstag are now loadable; tau_eff and efvs are not.
The subroutines that compute heat dissipation remain in glissade_therm.

The basal friction calculation remains in the velocity solver for Glissade.
For Glam, it has moved to glissade.F90, just after the velocity solve.
The relevant Glam subroutine, glam_basal_friction, is now in glam_velo.F90.

The units of dissip are now deg/s, consistently throughout the code.  Previously the units
switched between deg and deg/s.

In glam_velo.F90, I noticed that the call to glide_calcstrsstr (required for computing tau_eff)
took place after the Picard nonlinear solve but not the JFNK solve.  I moved the call so that
it takes place for either nonlinear solver.  This changes the answers for Glam evolving problems
using JFNK.

I decided to pass a number of arguments in and out of glissade_therm_driver rather than use pointers
to the model derived types.  This means	that glissade_therm_driver can use SI units without
having to rescale fields within	  the subroutine. (Input arguments can be passed with SI units.)
For now, glissade_therm_driver still uses model	units for time and thickness, but I will
switch entirely	to SI with a future commit (which will change answers at roundoff level).
Removed scale factors thk0 and tim0 from glissade_therm (answer-changing commit)

I modified glissade_therm_driver and associated subroutines to work entirely
in SI units.  The time step dttem now has units of seconds, and all thicknesses
have units of meters.  I removed thk0 and tim0 from these subroutines.
Input arguments (from glissade_tstep) are now passed in with SI units.

Answers change at roundoff level (around the 13th digit) for the evolving dome
and EISMINT-2 tests.  This is the first commit on this branch (glissade_thermal)
that changes answers for the standard evolving dome problem.

Note that the dissipation and flwa subroutines still work with scaled (non-SI) units.
That will probably change with a future commit.
Cleaned up glissade_therm subroutine that calculates flwa

I changed the name of glissade_therm_calcflwa to glissade_flow_factor and did
various cleanup.  Answers are unchanged for evolving dome (standard temp) and
EISMINT-2 (enthalpy scheme).

I rearranged the argument list for this subroutine and made some arguments optional.

I changed model%paramets%flow_factor to model%paramets%flow_fudge_factor, to avoid
confusion with the computed flow factor A.  (In the config file it is still called
'flow_factor', since I didn't want to edit all the SIA config files.)

I moved the call to glissade_flwa from glissade_tstep until after the halo updates for ice
geometry, to ensure that flwa will be accurate in halo cells.
Removed some scaling constants from glissade_therm.F90.

I converted three more subroutines (glissade_flow_factor and two heat
dissipation subroutines) to SI units.  This means that when these subroutines
are called from glissade_tstep, input arguments must be passed with SI units,
and in one case, output must be scaled back to model units.  The net effect
is to move scaling constants up a level, making glissade_tstep a bit messier
but glissade_therm cleaner.  Since glissade_therm is the module to be translated
to MPAS, I'm trying to make it as clean as possible.

With this commit, answers change at the machine roundoff level (around the 13th
decimal place) for the evolving dome and EISMINT-2 problems.
Cleaned up initialization routines in glissade_therm.F90

I modified subroutines glissade_init_therm and glissade_init_temp_column.
Subroutine glissade_init_therm now has a list of explicit arguments instead
the model derived type.

All subroutines in glissade_therm.F90 now assume SI units and have explicit
argument lists (instead of the model derived type).

I added a couple of halo calls in glissade_tstep before the call to glissade_init_therm,
to make sure that temperature is initialized correctly in all cells, including
the halo.

Temperature and some other variables are initialized to the unphysical value of -999.d0.
In glimmer_paramets.F90 I defined a new parameter, unphys_val = -999.d0.
Elsewhere in the code, variables are now set to unphys_val (instead of being set
to a literal constant).

Answers are unchanged relative to the previous commit.
Implemented the correct temperature basal BC for floating ice.

In glissade_therm.F90, I set the basal temperature for floating ice to tocnfrz,
a new parameter defined in glissade_physcon.F90.  I set tocnfrz = -1.8d0.

According to the documentation, this is the basal BC we use for ice shelves.
Until now, however, the shelf BC has not been treated correctly.
Implemented a small correction to basal melting for the enthalpy scheme.

Previously, the basal melt rate (bmlt) was always computed as

     bmlt = bflx / (rhoi*Li)

where bflx is the net heat flux available for melting, rhoi is ice density,
and Li is the latent heat of melting.

For the enthalpy scheme, I changed the denominator to (rhoi*Li - enth(upn)),
where enth(upn) is the specific enthalpy of the lowest ice layer.
This results in a slightly faster melt rate for temperate ice
(because of internal melting that has already occurred) and a
slightly slower rate for cold ice (because of the extra energy needed
to raise the temperature to T_pmp).  Note that enth = 0 for T = Tpmp
and waterfrac = 0, in which case the denominator reduces to rhoi*Li.

This change will make it easier to conserve energy exactly.
Some mods would be needed, however, to conserve energy in the rare case that
the basal melting extends beyond the lowest ice layer.

Basal melt rates for EISMINT-2 (with the enthalpy scheme) change by a fraction
of a percent. Since basal melt rates are ignored in the continuity equation by default,
ice state variables are unchanged.
Changes to reconcile old glissade_temp with new glissade_therm;	minor bmlt fix

In glissade_therm.F90, I added code to zero out bflx (also in the basal melt subroutine) in the case
that bflx is very slightly different from zero due to roundoff.  Without this code, the basal
temperature can be erroneously set to T = Tpmp - 0.01 on occasion.

Most of the changes made during the transition from glissade_temp to glissade_therm
have not been answer-changing except at the machine-roundoff level.  These are the exceptions:

(1) I fixed the indexing in a call to pmptemp in the basal melt subroutine.
    (Previously, the vertical index was off by one layer.)
(2) Zeroing out bflx as described above

In this commit I have added a logical variable called bmlt_bugfix in glissade_calcbmlt.
This is set to true by default.  If true, the two bug fixes are carried out and the
glissade_temp code should agree with the glissade_therm code (within machine roundoff)
when whichtemp = TEMP_PROGNOSIC = 1.  If false, the glissade_temp code should agree
(within machine roundoff) with the values obtained before I started working on
the temperature code.

Also, I recently moved the computations of dissip and bfricflx to glissade_tstep.
Here I am commenting out the old calls to the finddisp and calcbfric routines in glissade_temp,
to avoid double counting.  This switch can lead to roundoff-level changes relative
to the older location of the two calls.

In glissade_tstep I fixed the scaling of bmlt after the call to glissade_therm_driver.
I had forgotten	to multiply by tim0 in the conversion from m/s to model units,
resulting in erroneous diagnostics.  The correction fixes the bmlt diagnostics without
changing any other results (since bmlt by default does not feed	into the continuity equation).

In glide_bwater.F90 I added a comment.
Changes in parallel_slap; may already exist on 'develop' branch

I uncommented a 'stop' in subroutine parallel_stop so that serial runs will terminate correctly.
I think this change has already been made on 'develop', so it needn't be repeated.

I also added a TODO comment that probably is not on 'develop'.
Fixed a typo in the stream README file
Minor cleanup after bug fixes

In glissade_temp.F90, I removed the logical parameter that could be toggled
to run either with or without the bug fixes. Now, bug fixes are always applied.

Also, I verified that Glissade aborts if the user supplies a restart file
with 'temp' instead of 'tempstag' as an input field.  This didn't require
any code changes.

Answers are unchanged compared to the previous commit.
Removed or turned off some diagnostic print statements

The column diagnostics in glissade_therm_driver can be turned on by
setting verbose_temp = .true.
Minor mod in matrix factor; no answer change

I removed redundant factors of 2 in the numerator and denominator of
a matrix factor.  I thought this change would result in roundoff-level
answer changes, but it does not in my EISMINT-2 and dome test cases.
Changed some variable and subroutine names

In several files I changed "flow_fudge_factor" to "flow_enhancement_factor", given that
there are legitimate physical reasons for values greater than 1, so it is not simply a fudge.

In glissade_therm.F90, I changed the names of the two subroutines that
compute matrix elements for the vertical solve.

I removed a couple of debug print statements in glide_types.F90.
Allocate waterfrac and enthalpy whenever using the Glissade dycore.

Previously, waterfrace and enthalpy were allocated only when using the
option temperature = TEMP_ENTHALPY.  However, some subroutine shared
by the standard temperature and enthalpy schemes take waterfrac and/or
enthalpy as non-optional input arguments.

Rather than add logic for optional inputs, it was simpler to go ahead and
allocate the arrays (and initialize them to zero), especially given that
the enthalpy scheme may later be the default.
Implemented a depth-dependent freezing point of seawater as a shelf basal BC

The formula I'm using is

    tocnfrz = tocnfrz_sfc + dtocnfrz_dh * depth

where depth = thck * (rhoi/rhow)
      tocnfrz_sfc = -1.92 deg C
      dtocnfrz_dh = -7.53e-4 deg/m

These values are taken from the Ocean Water Freezing Point Calculator,
http://www.csgnetwork.com/h2ofreezecalc.html, assuming a salinity of
35 PSU and a pressure change of 1 kPa per meter depth.

This is not exact, but hopefully is good enough for government work.