Add more history dims

Merge branch 'track_cwd_in_fluxes' from ckoven

Added yet more history dims: CWD, fuel, canopy, canopy/leaf, and
canopy/leaf/pft

Added more dimensions to facilitate introspection of the model via the
history files. Two major areas here:

* 2 dimensions for the dead stuff: a CWD size-class dimension, and a
  fuel dimension (which is the CWD sizes plus 2 spots for leaves and
  grass).

* 3 new dimensions across the canopy vertical gradients to facilitate
  figuring out the light environment. most of the raw variables here
  are indexed by canopy level, leaf level, and pft, so I made a
  gigantic dimensions which multiplexes all three of these into a
  single dimension. since this is big (2x40x2 =160 currently with
  ed_npfts = 2, but will grow rapidly once we add more PFTs), I made
  two smaller dimensions as well: one that is canopy layer by leaf
  layer (so length 2x40 = 80) that averages the light environment
  across PFTs (and ought to be sufficient/overkill for most needs in
  the long run) and one that is just canopy dimension (length 2) for
  looking at things like the light environment of the top level of
  leaves for each canopy level, and which could also be used to
  collapse some of the canopy-level diagnostics that I had previously
  defined just as 1d onto so that they use 1 variable instead of 2.

Also added more diagnostic fluxes to track litter input fluxes due to
disturbance and termination, which were previously being passed
straight from live to dead pools without tracking the fluxes, and
which are now still passed directly from live to dead pools, but,
separately are being tracked and added to the litter_in and more
resolved fluxes. separately, i fixed unit errors on some of the
litter_in, litter_out, and other fluxes, but I've reverted the ones to
litter_in and litter_out to pass regression tests (since those two
variables are in the test) and will submit a second PR with just the
fixes to that shortly.

Added a bunch of variables on these new dimensions:

* for the cwd we're tracking the stocks, input and output fluxes on
the cwdsc dimensio, and for the fuel we are tracking fuel moisture on
the fuel dimension

* for the vertical canopy stuff we are tracking fabd, fabi, par and
lai for sun and shade leaves across various slices of the new
dims. for now these are all default on, but we can turn the larger
ones off when we have a better sense of having fixed the canopy light
environment.

Renamed nlevcan to nlevleaf, because it is a better descriptor for
what that variable represents -- the leaf sublayers within each canopy
level -- but also mainly because CLM already has a variable nlevcan
which collided with FATES's during the specification of the array
lengths in histfilemod.

updated 3/17 to note that @rgknox added one final dimension, which is
the 'classic' ED dimension of cohort size bin multiplexed by patch age
bin, and also did a bunch of cleanup, gave better names to the
multiplexed dimension names, and replaced a bunch of inline unit
changes with actual named unit conversion constants

Added a named constant to clm history module to set length of
dimension names.

User interface changes?: no

Code review: Some review already by @rgknox and @rosiealice, could do
with some more. note that this is on top of PR #184 so best to wait
until that is pulled so that these new changes are more apparent.

Testing:
  ckoven:
    Test suite: Ed suite on lawrencium
    Test baseline: a5dc8da
    Test namelist changes:
    Test answer changes: bit for bit

    Test summary: all pass.

  andre:
    Test suite: ed - yellowstone gnu, intel, pgi
                     hobart nag
    Test baseline: e99fcf9
    Test namelist changes: none
    Test answer changes: bit for bit
    Test summary: all tests pass

    Test suite: clm_short - yellowstone gnu, intel, pgi
    Test baseline: clm4_5_12_r195
    Test namelist changes: none
    Test answer changes: bit for bit
    Test summary: all tests pass

components/clm/src/ED/biogeochem/EDCanopyStructureMod.F90

@@ -12,9 +12,11 @@ module EDCanopyStructureMod
   use EDCohortDynamicsMod   , only : copy_cohort, terminate_cohorts, fuse_cohorts
   use EDtypesMod            , only : ed_site_type, ed_patch_type, ed_cohort_type, ncwd
   use EDTypesMod            , only : nclmax
-  use EDTypesMod            , only : nlevcan
+  use EDTypesMod            , only : nlevleaf
   use EDTypesMod            , only : numpft_ed
+  use EDtypesMod            , only : AREA
   use FatesGlobals          , only : endrun => fates_endrun
+  use FatesInterfaceMod     , only : hlm_days_per_year
 
   ! CIME Globals
   use shr_log_mod           , only : errMsg => shr_log_errMsg
@@ -96,10 +98,10 @@ subroutine canopy_structure( currentSite )
     real(r8) :: cc_loss
     real(r8) :: lossarea
     real(r8) :: newarea
-    real(r8) :: arealayer(nlevcan) ! Amount of plant area currently in each canopy layer
-    real(r8) :: sumdiff(nlevcan)   ! The total of the exclusion weights for all cohorts in layer z 
+    real(r8) :: arealayer(nlevleaf) ! Amount of plant area currently in each canopy layer
+    real(r8) :: sumdiff(nlevleaf)   ! The total of the exclusion weights for all cohorts in layer z 
     real(r8) :: weight                ! The amount of the total lost area that comes from this cohort
-    real(r8) :: sum_weights(nlevcan)
+    real(r8) :: sum_weights(nlevleaf)
     real(r8) :: new_total_area_check
     real(r8) :: missing_area, promarea,cc_gain,sumgain
     integer  :: promswitch,lower_cohort_switch
@@ -229,14 +231,33 @@ subroutine canopy_structure( currentSite )
 
                             enddo
 
-                                currentPatch%leaf_litter(currentCohort%pft)  = &
-                                     currentPatch%leaf_litter(currentCohort%pft) + (currentCohort%bl)* &
-                                          currentCohort%n/currentPatch%area ! leaf litter flux per m2.
-
-                                currentPatch%root_litter(currentCohort%pft)  = &
-                                     currentPatch%root_litter(currentCohort%pft) + &
-                                     (currentCohort%br+currentCohort%bstore)*currentCohort%n/currentPatch%area
-
+                            currentPatch%leaf_litter(currentCohort%pft)  = &
+                                 currentPatch%leaf_litter(currentCohort%pft) + (currentCohort%bl)* &
+                                 currentCohort%n/currentPatch%area ! leaf litter flux per m2.
+
+                            currentPatch%root_litter(currentCohort%pft)  = &
+                                 currentPatch%root_litter(currentCohort%pft) + &
+                                 (currentCohort%br+currentCohort%bstore)*currentCohort%n/currentPatch%area
+
+                            ! keep track of the above fluxes at the site level as a CWD/litter input flux (in kg / site-m2 / yr)
+                            do c=1,ncwd
+                               currentSite%CWD_AG_diagnostic_input_carbonflux(c) = &
+                                    currentSite%CWD_AG_diagnostic_input_carbonflux(c) &
+                                    + currentCohort%n*(currentCohort%bdead+currentCohort%bsw) * &
+                                    SF_val_CWD_frac(c) * ED_val_ag_biomass * hlm_days_per_year / AREA
+                               currentSite%CWD_BG_diagnostic_input_carbonflux(c) = &
+                                    currentSite%CWD_BG_diagnostic_input_carbonflux(c) &
+                                    + currentCohort%n*(currentCohort%bdead+currentCohort%bsw) * &
+                                    SF_val_CWD_frac(c) * (1.0_r8 -  ED_val_ag_biomass)  * hlm_days_per_year / AREA
+                            enddo
+
+                            currentSite%leaf_litter_diagnostic_input_carbonflux(currentCohort%pft) = &
+                                 currentSite%leaf_litter_diagnostic_input_carbonflux(currentCohort%pft) +  &
+                                 currentCohort%n * (currentCohort%bl) * hlm_days_per_year  / AREA
+                            currentSite%root_litter_diagnostic_input_carbonflux(currentCohort%pft) = &
+                                 currentSite%root_litter_diagnostic_input_carbonflux(currentCohort%pft) + &
+                                 currentCohort%n * (currentCohort%br+currentCohort%bstore) * hlm_days_per_year  / AREA
+
                             currentCohort%n = 0.0_r8
                             currentCohort%c_area = 0._r8
                          else  
@@ -280,13 +301,33 @@ subroutine canopy_structure( currentSite )
 
                             enddo
 
-                                currentPatch%leaf_litter(currentCohort%pft)  = &
-                                     currentPatch%leaf_litter(currentCohort%pft) + currentCohort%bl* &
-                                          currentCohort%n/currentPatch%area ! leaf litter flux per m2.
+                            currentPatch%leaf_litter(currentCohort%pft)  = &
+                                 currentPatch%leaf_litter(currentCohort%pft) + currentCohort%bl* &
+                                 currentCohort%n/currentPatch%area ! leaf litter flux per m2.
+
+                            currentPatch%root_litter(currentCohort%pft)  = &
+                                 currentPatch%root_litter(currentCohort%pft) + &
+                                 (currentCohort%br+currentCohort%bstore)*currentCohort%n/currentPatch%area
+
+                            ! keep track of the above fluxes at the site level as a CWD/litter input flux (in kg / site-m2 / yr)
+                            do c=1,ncwd
+                               currentSite%CWD_AG_diagnostic_input_carbonflux(c) = &
+                                    currentSite%CWD_AG_diagnostic_input_carbonflux(c) &
+                                    + currentCohort%n*(currentCohort%bdead+currentCohort%bsw) * &
+                                    SF_val_CWD_frac(c) * ED_val_ag_biomass * hlm_days_per_year / AREA
+                               currentSite%CWD_BG_diagnostic_input_carbonflux(c) = &
+                                    currentSite%CWD_BG_diagnostic_input_carbonflux(c) &
+                                    + currentCohort%n*(currentCohort%bdead+currentCohort%bsw) * &
+                                    SF_val_CWD_frac(c) * (1.0_r8 -  ED_val_ag_biomass)  * hlm_days_per_year / AREA
+                            enddo
+
+                            currentSite%leaf_litter_diagnostic_input_carbonflux(currentCohort%pft) = &
+                                 currentSite%leaf_litter_diagnostic_input_carbonflux(currentCohort%pft) +  &
+                                 currentCohort%n * (currentCohort%bl) * hlm_days_per_year  / AREA
+                            currentSite%root_litter_diagnostic_input_carbonflux(currentCohort%pft) = &
+                                 currentSite%root_litter_diagnostic_input_carbonflux(currentCohort%pft) + &
+                                 currentCohort%n * (currentCohort%br+currentCohort%bstore) * hlm_days_per_year  / AREA
 
-                                currentPatch%root_litter(currentCohort%pft)  = &
-                                     currentPatch%root_litter(currentCohort%pft) + &
-                                     (currentCohort%br+currentCohort%bstore)*currentCohort%n/currentPatch%area
                             currentCohort%n = 0.0_r8
                             currentCohort%c_area = 0._r8
 
@@ -635,7 +676,7 @@ subroutine canopy_spread( currentSite )
     ! !LOCAL VARIABLES:
     type (ed_cohort_type), pointer :: currentCohort
     type (ed_patch_type) , pointer :: currentPatch
-    real(r8) :: arealayer(nlevcan) ! Amount of canopy in each layer. 
+    real(r8) :: arealayer(nlevleaf) ! Amount of canopy in each layer. 
     real(r8) :: inc                   ! Arbitrary daily incremental change in canopy area 
     integer  :: z
     !----------------------------------------------------------------------
@@ -694,7 +735,7 @@ subroutine canopy_summarization( nsites, sites, bc_in )
     use FatesInterfaceMod    , only : bc_in_type
     use EDPatchDynamicsMod   , only : set_patchno
     use EDPatchDYnamicsMod   , only : set_root_fraction
-    use EDCohortDynamicsMod  , only : size_and_type_class_index
+    use EDTypesMod           , only : sizetype_class_index
     use EDGrowthFunctionsMod , only : tree_lai, c_area
     use EDEcophysConType     , only : EDecophyscon
     use EDtypesMod           , only : area
@@ -750,8 +791,8 @@ subroutine canopy_summarization( nsites, sites, bc_in )
 
              ! Update the cohort's index within the size bin classes
              ! Update the cohort's index within the SCPF classification system
-             call size_and_type_class_index(currentCohort%dbh,currentCohort%pft, &
-                                            currentCohort%size_class,currentCohort%size_by_pft_class)
+             call sizetype_class_index(currentCohort%dbh,currentCohort%pft, &
+                                       currentCohort%size_class,currentCohort%size_by_pft_class)
 
 
              currentCohort%b = currentCohort%balive+currentCohort%bdead+currentCohort%bstore
@@ -769,21 +810,25 @@ subroutine canopy_summarization( nsites, sites, bc_in )
 
              ! Check for erroneous zero values. 
              if(currentCohort%dbh <= 0._r8 .or. currentCohort%n == 0._r8)then
-                write(fates_log(),*) 'ED: dbh or n is zero in canopy_summarization', currentCohort%dbh,currentCohort%n
+                write(fates_log(),*) 'ED: dbh or n is zero in canopy_summarization', &
+                      currentCohort%dbh,currentCohort%n
              endif
              if(currentCohort%pft == 0.or.currentCohort%canopy_trim <= 0._r8)then
-                write(fates_log(),*) 'ED: PFT or trim is zero in canopy_summarization',currentCohort%pft,currentCohort%canopy_trim
+                write(fates_log(),*) 'ED: PFT or trim is zero in canopy_summarization', &
+                      currentCohort%pft,currentCohort%canopy_trim
              endif
              if(currentCohort%balive <= 0._r8)then
-                write(fates_log(),*) 'ED: balive is zero in canopy_summarization',currentCohort%balive
+                write(fates_log(),*) 'ED: balive is zero in canopy_summarization', &
+                      currentCohort%balive
              endif
 
              currentCohort => currentCohort%taller
 
           enddo ! ends 'do while(associated(currentCohort))
 
           if ( currentPatch%total_canopy_area-currentPatch%area > 0.000001_r8 ) then
-             write(fates_log(),*) 'ED: canopy area bigger than area',currentPatch%total_canopy_area ,currentPatch%area
+             write(fates_log(),*) 'ED: canopy area bigger than area', &
+                   currentPatch%total_canopy_area ,currentPatch%area
              currentPatch%total_canopy_area = currentPatch%area
           endif
 
@@ -1136,10 +1181,10 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si)
                          /currentPatch%tlai_profile(L,ft,iv)
                 enddo
 
-                currentPatch%tlai_profile(L,ft,currentPatch%nrad(L,ft)+1: nlevcan) = 0._r8
-                currentPatch%tsai_profile(L,ft,currentPatch%nrad(L,ft)+1: nlevcan) = 0._r8
-                currentPatch%elai_profile(L,ft,currentPatch%nrad(L,ft)+1: nlevcan) = 0._r8 
-                currentPatch%esai_profile(L,ft,currentPatch%nrad(L,ft)+1: nlevcan) = 0._r8
+                currentPatch%tlai_profile(L,ft,currentPatch%nrad(L,ft)+1: nlevleaf) = 0._r8
+                currentPatch%tsai_profile(L,ft,currentPatch%nrad(L,ft)+1: nlevleaf) = 0._r8
+                currentPatch%elai_profile(L,ft,currentPatch%nrad(L,ft)+1: nlevleaf) = 0._r8 
+                currentPatch%esai_profile(L,ft,currentPatch%nrad(L,ft)+1: nlevleaf) = 0._r8
 
              enddo
           enddo

components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90

@@ -9,6 +9,7 @@ module EDCohortDynamicsMod
   use FatesInterfaceMod     , only : hlm_freq_day
   use FatesConstantsMod     , only : r8 => fates_r8
   use FatesConstantsMod     , only : fates_unset_int
+  use FatesInterfaceMod     , only : hlm_days_per_year
   use EDPftvarcon           , only : EDPftvarcon_inst
   use EDEcophysContype      , only : EDecophyscon
   use EDGrowthFunctionsMod  , only : c_area, tree_lai
@@ -20,6 +21,7 @@ module EDCohortDynamicsMod
   use EDTypesMod            , only : sclass_ed,nlevsclass_ed,AREA
   use EDTypesMod            , only : min_npm2, min_nppatch
   use EDTypesMod            , only : min_n_safemath
+  use EDTypesMod            , only : sizetype_class_index
   ! CIME globals
   use shr_log_mod           , only : errMsg => shr_log_errMsg
   !
@@ -35,7 +37,6 @@ module EDCohortDynamicsMod
   public :: sort_cohorts
   public :: copy_cohort
   public :: count_cohorts
-  public :: size_and_type_class_index
   public :: allocate_live_biomass
 
   logical, parameter :: DEBUG  = .false. ! local debug flag
@@ -104,8 +105,8 @@ subroutine create_cohort(patchptr, pft, nn, hite, dbh, &
     new_cohort%balive       = balive
     new_cohort%bstore       = bstore
 
-    call size_and_type_class_index(new_cohort%dbh,new_cohort%pft, &
-                                   new_cohort%size_class,new_cohort%size_by_pft_class)
+    call sizetype_class_index(new_cohort%dbh,new_cohort%pft, &
+                              new_cohort%size_class,new_cohort%size_by_pft_class)
 
     if ( DEBUG ) write(fates_log(),*) 'EDCohortDyn I ',bstore
 
@@ -484,7 +485,7 @@ subroutine zero_cohort(cc_p)
   end subroutine zero_cohort
 
   !-------------------------------------------------------------------------------------!
-  subroutine terminate_cohorts( siteptr, patchptr )
+  subroutine terminate_cohorts( currentSite, patchptr )
     !
     ! !DESCRIPTION:
     ! terminates cohorts when they get too small      
@@ -494,7 +495,7 @@ subroutine terminate_cohorts( siteptr, patchptr )
     use SFParamsMod, only : SF_val_CWD_frac
     !
     ! !ARGUMENTS    
-    type (ed_site_type), intent(inout), target :: siteptr
+    type (ed_site_type) , intent(inout), target :: currentSite
     type (ed_patch_type), intent(inout), target :: patchptr
     !
     ! !LOCAL VARIABLES:
@@ -571,10 +572,10 @@ subroutine terminate_cohorts( siteptr, patchptr )
           else
              levcan = 2
           endif
-          siteptr%terminated_nindivs(currentCohort%size_class,currentCohort%pft,levcan) = &
-               siteptr%terminated_nindivs(currentCohort%size_class,currentCohort%pft,levcan) + currentCohort%n
+          currentSite%terminated_nindivs(currentCohort%size_class,currentCohort%pft,levcan) = &
+               currentSite%terminated_nindivs(currentCohort%size_class,currentCohort%pft,levcan) + currentCohort%n
           !
-          siteptr%termination_carbonflux(levcan) = siteptr%termination_carbonflux(levcan) + &
+          currentSite%termination_carbonflux(levcan) = currentSite%termination_carbonflux(levcan) + &
                currentCohort%n * currentCohort%b
           if (.not. associated(currentCohort%taller)) then
              currentPatch%tallest => currentCohort%shorter
@@ -604,6 +605,23 @@ subroutine terminate_cohorts( siteptr, patchptr )
              currentPatch%root_litter(currentCohort%pft) = currentPatch%root_litter(currentCohort%pft) + currentCohort%n* &
                   (currentCohort%br+currentCohort%bstore)/currentPatch%area 
 
+             ! keep track of the above fluxes at the site level as a CWD/litter input flux (in kg / site-m2 / yr)
+             do c=1,ncwd
+                currentSite%CWD_AG_diagnostic_input_carbonflux(c)  = currentSite%CWD_AG_diagnostic_input_carbonflux(c) &
+                     + currentCohort%n*(currentCohort%bdead+currentCohort%bsw) * &
+                     SF_val_CWD_frac(c) * ED_val_ag_biomass * hlm_days_per_year / AREA
+                currentSite%CWD_BG_diagnostic_input_carbonflux(c)  = currentSite%CWD_BG_diagnostic_input_carbonflux(c) &
+                     + currentCohort%n*(currentCohort%bdead+currentCohort%bsw) * &
+                     SF_val_CWD_frac(c) * (1.0_r8 -  ED_val_ag_biomass)  * hlm_days_per_year / AREA
+             enddo
+
+             currentSite%leaf_litter_diagnostic_input_carbonflux(currentCohort%pft) = &
+                  currentSite%leaf_litter_diagnostic_input_carbonflux(currentCohort%pft) +  &
+                  currentCohort%n * (currentCohort%bl) * hlm_days_per_year  / AREA
+             currentSite%root_litter_diagnostic_input_carbonflux(currentCohort%pft) = &
+                  currentSite%root_litter_diagnostic_input_carbonflux(currentCohort%pft) + &
+                  currentCohort%n * (currentCohort%br+currentCohort%bstore) * hlm_days_per_year  / AREA
+
              deallocate(currentCohort)     
           endif
        endif
@@ -619,7 +637,7 @@ subroutine fuse_cohorts(patchptr)
     ! Join similar cohorts to reduce total number            
     !
     ! !USES:
-    use EDTypesMod  , only :  nlevcan
+    use EDTypesMod  , only :  nlevleaf
     !
     ! !ARGUMENTS    
     type (ed_patch_type), intent(inout), target :: patchptr
@@ -775,7 +793,7 @@ subroutine fuse_cohorts(patchptr)
                          currentCohort%canopy_layer_yesterday  = (currentCohort%n*currentCohort%canopy_layer_yesterday  + &
                               nextc%n*nextc%canopy_layer_yesterday)/newn
 
-                         do i=1, nlevcan     
+                         do i=1, nlevleaf     
                             if (currentCohort%year_net_uptake(i) == 999._r8 .or. nextc%year_net_uptake(i) == 999._r8) then
                                currentCohort%year_net_uptake(i) = min(nextc%year_net_uptake(i),currentCohort%year_net_uptake(i))
                             else
@@ -1178,25 +1196,6 @@ function count_cohorts( currentPatch ) result ( backcount )
 
   end function count_cohorts
 
-  ! =====================================================================================
-
-  subroutine size_and_type_class_index(dbh,pft,size_class,size_by_pft_class)
-
-    use EDTypesMod, only: sclass_ed
-    use EDTypesMod, only: nlevsclass_ed
-
-    ! Arguments
-    real(r8),intent(in) :: dbh
-    integer,intent(in)  :: pft
-    integer,intent(out) :: size_class
-    integer,intent(out) :: size_by_pft_class
-
-    size_class        = count(dbh-sclass_ed.ge.0.0_r8)
-
-    size_by_pft_class = (pft-1)*nlevsclass_ed+size_class
-
-    return
-  end subroutine size_and_type_class_index