1. Fixed a few minor bugs on Monin-Obukhov and changed the Massman scheme to have a drag
   profile. 
2. Re-wrote the two-stream model, and but I haven't tested. Moving disks because I can't 
   compile on /n/gstore
1. Final fixes on the updated two-stream model.  Next commit will contain the old two-stream
   for comparison, plus some new fixes on my own paths
1.  The new two-stream model is available and working.  This method looks a lot more like
    the other scheme (Zhao and Qualls), and it seems to accomodate higher/more reasonable
    LAI without the crown model, which is inconsistent in the original formulation since 
    it only affects the transmissivity of direct radiation, violating the principle that 
    direct and diffuse radiation have the same optical properties for a beam of same 
    wavelength coming from the same direction.  I think in the future we should have only
    one "two-stream" model, but for the time being ICANRAD=2 uses the new two-stream and
    ICANRAD=0. 
2.  Various minor modifications/bug fixes to the code:
    a.  Changes to snow/temporary surface water model:
        i.   Snow liquid/ice ratio from precipitation is now continuous.  This is just a
             very minor modification from the original scheme between Ta=2.0C and Ta=2.5C.
        ii.  The fraction of ground covered by water had no obvious reference, so I 
             switched it to the same scheme used in CLM-4.
        iii. The fraction of ground covered by temporary surface water is now applied
             consistently throughout the code.  Before this was only used to determine 
             the roughness, but it is now applied to radiation and heat and water fluxes
             with the canopy air space as well.  This means that the ground flux to CAS
             can have contribution from both the soil and the temporary surface water.
        iv.  The temporary surface water albedo is now interpolated from the top soil
             albedo to the fresh snow albedo.  The top soil albedo is dependent on soil 
             colour (it has been like this for a while, but it wasn't consistent with
             the temporary  water albedo).
        v.   Wish list: snow should change albedo and density as it ages, CLM-4 and ECHAM-5
             have some relatively simple formulae to account for this.  Also, snow 
             radiation should be solved through two-stream, at least for PAR and NIR.
    b.  Changes in Massman (1997, ICANTURB = 2) and leaf/wood boundary layer model:
        i.   The drag coefficient divided by the sheltering factor are modelled by 
             Wohfahrt and Cernusca (2002), following their argument that CD/Ps are not
             really separable.
        ii.  Vegetation wind speed is now the average between the bottom and the top of
             the crown as opposed to the wind at the mid-crown height.  The mid-point 
             bias the wind speed (too low), since the profile follows an exponential
             decay.
        iii. Kinematic viscosity, thermal diffusivity are now functions of the canopy air 
             space temperature, following Monteith and Unsworth (2008).  Thermal
             expansion is assumed to be 1/T, which is the coefficient assuming perfect 
             gas (Dufour and van Mieghem, 1975). 
    c.  Radiation:  for long-wave radiation, the layer temperature is now a weighted 
        average of T^4, and the weights account for both area index and the emissivity.
    d.  Photosynthesis:  do not let GPP be negative, which causes some problems with 
        patch dynamics when running with the debugger.
    e.  Surface layer model:  The dimensionless scale had a "zeta0" term missing. 
    f.  Soil thermal conductivity: the original method came from LEAF-3 but had no 
        reference whatsoever (it wasn't in LEAF-2 or Bob Walko's paper).  I changed to 
        the method suggested by Monteith and Unsworth (2008), which is closer to what 
        SiB does: this means that all parameters are determined from slxsand and slxclay.
    g.  Patch fusion/fission: averaged quantities are not fused during initialisation, 
        they are all zero (and derived quantities such as temperature cannot be determined,
        due to divisions by zero). 
3.  RYAN'S MEMORY LEAK FIX! Thanks, Ryan!
1.  Two minor bug fixes in the calculation of vegetation conductance.  The integral should
    be done on the heat diffusivity, not momentum , and the top boundary condition should
    account for non-neutral conditions, hence the extra term
    (1 - zetatop * dpsimdzeta(zetatop,stable))
2.  Minor fix in make_joborder.r, it had a couple of lines for one specific run.
1. Minor bug fixes in the coupled model, these appeared after the latest modifications on
   ED-2 biophysics and were mostly typos or forgotten use statements.
1.  Tuned soil decomposition parameters for decomp_scheme = 2.  Although this model yields 
    much more reasonable soil carbon, long-term runs were underestimating soil carbon for 
    forests and overestimating soil carbon for savannahs. Still testing it, though.
2.  Fixed some problems in the job management scripts.
1.  Big-leaf ED bug fixes:
    a.  Mortality due to fires were not incorporated in my version.  Now I use
        cpatch%mort_rate(5,ico) to find the change in nplant (disturbance mortality). 
        cpatch%mort_rate(5,ico) is used only when big leaf is on, otherwise it is just a
        diagnostic variable
    b.  Somewhere the update of energy was wrong in big leaf: it assumed that leaf/wood
        water was zero but it wasn't always zero, causing crashing in the thermodynamics.
        Now all places in the code use cmtl2uext to update energy, which is the built-in
        function that updates energy
    c.  Area was sometimes greater than one (!).  I changed rescale_patches so it runs
        even when there is a single PFT left, in case the area is not 1.
2.  Output fixes:
    a. LAI_PY and WAI_PY were incorrectly aggregated, removed the nplant that was 
       multiplying the cohort-level variable
    b. In radiate_driver.f90:  the code was trying to update some daily mean variables 
       even when IDOUTPUT/IMOUTPUT/IQOUTPUT were not requested.
3.  Added some fire survivorship parameters.  Currently they are all set to 0, so nothing
    changes, but in the future if we want to add fire-resistant PFTs they may be set 
    to something else.
4.  Blocked CROWN_MOD = 1 when ICANRAD = 1 or 2.  There is one term (mubar) that must be 
    correctly aggregated before using it.
5.  Few updates on shell scripts so regional runs can also use run_sitter.sh
1.  Bug fix in the new two stream radiation.  The bottom boundary condition was referring
    incorrectly using the direct radiation from the top of the canopy as opposed to the 
    bottom of the canopy.
2.  Added a new variable rad_profile, with the actual radiation profile (not the absorbed
    radiation).  The variable is a cpatch variable with dimensions (10,ncohorts)
     1 - Downward PAR, direct
     2 - Upward   PAR, direct (it is defined as zero for ICANRAD=1 and 2)
     3 - Downward PAR, diffuse
     4 - Upward   PAR, diffuse
     5 - Downward NIR, direct
     6 - Upward   NIR, direct (it is defined as zero for ICANRAD=1 and 2)
     7 - Downward NIR, diffuse
     8 - Upward   NIR, diffuse
     9 - Downward TIR, diffuse
    10 - Upward   TIR, diffuse
3.  Minor bug fixes in some output variables.
1. No change in EDBRAMS, only updates in the R scripts and shell scripts.
1. Bug fix in phenology.  Storage is now properly accounted for in the carbon balance.
2. Many minor changes in multiple R scripts. cannot remember all of them.
1.  Fixed the patch fusion, there was an extra logical flag that was forcing fusion even 
    when the number of patches was less than MAXPATCH.
2.  Added a cap in frost mortality rate, following Albani et al. (2006).
3.  Several updates and fixes in R scripts and sh.
1. Changes in the scripts to use sandy
2. Compiler settings for CentOS 6 (sandy's OS)
3. Few R script changes
Minor updates in some post-processing scripts
1.  Bug fix in output variable LAI_PY and WAI_PY (polygon-level).  They were incorrectly
    scaled by nplant (excessive faith on ctrl+c ctrl+v...)
2.  Added Ryan.s improvements on cumulus parametrisation
3.  Added a few new scripts that describe the different scenarios.
4.  Several minor updates in many R scripts (not directly linked to ED)
Quick update of a few R scripts (so it updates the back up)
1.  Bug fix in varf_read, one variable was declared with the wrong dimensions.
2.  Increased the number of isentropic analysis files, so it doesn't crash when 
    you try to run multiple years.
3.  Updates on R scripts (my stuff only)
1. Switched phenology_status = 2 to phenology_status = -2
   (it makes more sense to be negative)
2. Minor bug fixes on outgoing radiation in BRAMS; diagnostic only
3. Minor bug fixes for carbon balance when running with vegetation dynamics off.  Now it 
   also calculates growth rates even though it does not change size.
4. Several updates on R scripts.
1. Minor modification in utility scripts for BRAMS, now than can be used in .crontab
2. Minor updates on R scripts.  Compare_hourly.r now incorporates many of the comparisons
   from plot_eval_ed.r, but it puts everything in the same plot so it is easier to compare
   sites and simulations.
1. Means for growth, storage and "vleaf" respiration were incorrectly written in kgC/m2/day
   instead of kgC/m2/yr.  We should seriously standardise units in this model...
2. Updated some paths and some R scripts.
3. Turned off the annoying history file message.
1. Updated EFT periods for Paracou and Manaus.  Included the option of sub-hourly "scenario"