-  Reproduction, seed_dispersal fix:  nseed_maygo was incorrectly weighted

-  Change phenology so that if phenology_status flag is updated every day.

-  Terminate_patches subroutine called in vegetation_dyanamics

NOTE: THIS REVISION HAS A NON-WORKING VERSION OF BIG-LEAF ED AS AN OPTION
IN ED2IN.  THIS FEATURE IS UNDER-DEVELOPMENT AND WILL BE FULLY FUNCTIONAL
SOON

R Knox Changes:
 ------ Hybrid Solver ----
 This commit encapsulates hybrid time stepping.  This uses an euler
 forward step for everything but leaf and canopy temperature. Leaf 
 and canopy temperature use a bdf2 implicit step.  Boundary
 conditions for the implicit step use the following:
 -forward step leaf water and canopy humidity.
 -current step ground temperature, radiation, leaf liquid state,
          atmospheric temperature, aerodynamic conductances

 Future versions of the hybrid solver could work out using forward step
 radiation (cosz approx), and solving for the soil surface water
balancing
 prior to the implicit step.  This is done after the implicit step,
prohibiting
 a good estimate of forward state surface water and therefore
disallowing
 forward step ground temperature.

 Two new files are added, their names are self explanatory.
 hybrid_driver.f90
 bdf2_solver.f90

 Hybrid time-stepping is now option number 3 in the ED2IN card.

 Next step:  incorporating hybrid stepping in the coupled model.

 Added some text to the ED2IN describing the hybrid option
 Amended bdf2_solver to solve for wood also.  This is untested and
needs discussion.  Use    ibranch_thermo=0 with hybrid stepping until this
is tested.

-also added option to read ED runs done with the late model (from Ke)
 but removing the water sites


MD Longo Changes (mdl rv 289):
1. Bug fix in energy budget: the internal energy that left the soil
layers to be transpired was removed from the soil but not added
anywhere else. Because phase change happens at the leaves, I added
this energy to the leaves so the energy budget is closed

2. Bug fix in the mixed branch thermodynamics (IBRANCH_THERMO=1). In
case one part was resolvable and the other wasnt, the model was
forcing the non-resolvable part to be in thermal equilibrium with the
canopy air space when it left the integrator. This was causing
instabilities, because inside the integrator both pools are bound
together. Now we always impose thermal equilibrium between leaves
and wood, and make them in thermal equilibrium with canopy air space
only when neither the leaves nor the wood are resolvable.

3. Switched the prognostic variable of the canopy air space to
specific enthalpy, which now includes the effect of moisture. There
were a few reasons for this change:
   a. We assume constant pressure within one integration call. From
the first law of thermodynamics change in enthalpy is equal to heat
exchange when pressure remains constant, which makes enthalpy a
natural prognostic variable. The caveat is that enthalpy is not
preserved between time steps, so we still track potential temperature,
which is conserved when no phase change and no energy is exchanged,
but pressure changes, like between two consecutive calls.
   b. The energy associated with water vapor was not part of the budget
so transpiration, and evaporation from ground and vegetation were both
treated as sinks of energy, while condensation was treated as source.
This is not true because the energy remains in the system, at the
canopy air space.
   c. The amount of energy transferred between the surfaces and the CAS
was incorrect. Phase change requires latent heat, but in ED once
water evaporates from the surfaces, not only does it change phase, but
it also moves from one pool to the other, so the total energy
exchanged should include the loss of internal energy, and by using
enthalpy we can track the total energy exchanged with simpler equations.

4. Because of the enthalpy, we redefined the latent heat of
sublimation and evaporation. They used to be constants, which is
thermodynamically inconsistent. Because we assume the specific heat
to be constants, the latent heat must be assumed a linear function of
temperature so latent heat becomes the difference in enthalpy between
the two phases involved. The original terms were renamed as alvl3 and
alvi3, to represent that they are values at the triple point, and alvl
and alvi became functions

5. I moved all conversions between enthalpy/energy and temperature to
therm_lib and therm_lib8, and renamed some of the functions (qtk,
qwtk) to better represent their thermodynamic meaning.

Minor changes;
- nonlocal_dispersal wasnt assigned for PFTs (12-15)
- I just added the 2 variables to control the detailed output, a flag
to tell which types of output the user wants (IDETAILED) and
PATCH_KEEP, a flag to define which patch remains (used only when
IDETAILED is not 0). Most detailed output works only for single
polygon runs.
- Ramspost: now the buffers have enough size to deal with 4-D cumulus
variables, and cumulus variables seem to be working fine.\
- Reduced the photosynthesis tolerance to the same as therm_lib8
- Bug fix in canopy_structure_dynamics so that Massaman+branch thermo works
-  IALLOM parameters for dbh2bl to the way they were originally intended
-  Changed vt_vector from pointer to allocatable array
-  Reduced the number of pointer variables to 1100
-  Removed obsolete variables like wood projected area