Several improvements to treat soil freezing when using Richards equation

1. Introduction of a function to initialize the soil in hydraulic and thermal equilibrium upon first function call of Richards Equation, to be switched on by setting REQ_INITIALIZE_SOIL = TRUE in [SnowpackAdvanced]. If activated, the pressure head is forced into steady state to minimize water flow at the first time step. Then liquid water and ice content is partitioned to be in thermal equilibrium.
2. Making the method to describe the effect on ice on the hydraulic conductivity a configurable option via the key HYDRAULIC_CONDUCTIVITY_FROZEN_SOIL. Default would be IGNORE, which would be consistent with SNOWPACK thus far. Other options are OMEGA (Eq. 6 in Hansson et al. 2004) and LIQUIDPORESPACE. The latter simply calculates hydraulic conductivity using saturation calculation using the liquid water part only. It is inhibiting downward water percolation less, and thereby is numerically more stable.
3. A few small tweaks to the code for better treatment of internal energy. Richards equation solver should not touch the node temperatures, and use Qmf to describe total change in energy in a layer. PhaseChange.cc will then take care of adjusting the element and nodal temperatures to remain as consistent as possible. A few tweaks were made to make this behavior more consistent.