new mortality and canopy/understory diagnostics and option for seed h…

…omogenization

Merge branch 'new_mort_diagnostics'

added a bunch of new diagnostics and also a mode to prevent
competitive exclusion.

the new diagnostics cover:

* canopy vs. understory resolved population numbers, mortality, growth
  increment, GPP, and respiration fluxes so as to see what's going on
  in the understory separately from what's going on in the canopy

* now tracking cohort termination (due to both numerical truncation
  but more importantly also when cohorts' storage pools go negative)
  as a sixth mortality term, which has been occurring all along, but
  was not being tracked as a mortality term, and turns out to be
  important to the population budget

* now tracking recruitment, so that we can ensure conservation of
  individuals as dn/dt = recruitment - sum of all mortality terms

* redefined the canopy_area_by_patch_age diagnostic to use
  cohort-level rather than patch-level canopy area (since patch-level
  canopy area saturates at the patch area, whereas cohort-level canopy
  area doesn't so that the ratio of canopy_area_by_patch_age /
  patch_area_by_age becomes a useful metric of that space from
  pre-closure through multiple canopy levels

* the seed homogenization is default off but if set to be on prevents
  competitive exclusion, like seed rain, but without requiring a
  non-balanced carbon budget.

Fixes: #181

User interface changes?: No

Code review: Requesting @rosiealice and @rgknox review.

Testing:

  rgknox / ckoven

    Test suite: ed suite on lawrencium-lr3
    Test baseline: a5dc8da
    Test namelist changes:
    Test answer changes: bit for bit

    Test summary: all PASS

andre:

    Test suite: ed - yellowstone gnu, intel, pgi
                ed - hobart nag
    Test baseline: b0bceb0
    Test namelist changes: none?
    Test answer changes: bit for bit

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

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

@@ -107,10 +107,16 @@ subroutine canopy_structure( currentSite )
     integer  :: count_mi
     !----------------------------------------------------------------------
 
-    currentPatch => currentSite%oldest_patch
-
+    currentPatch => currentSite%oldest_patch    
+    ! 
+    ! zero site-level demotion / promotion tracking info
+    currentSite%demotion_rate(:) = 0._r8
+    currentSite%promotion_rate(:) = 0._r8
+    currentSite%demotion_carbonflux = 0._r8
+    currentSite%promotion_carbonflux = 0._r8
+    !
     ! Section 1: Check  total canopy area.    
-
+    !
     new_total_area_check = 0._r8
     do while (associated(currentPatch)) ! Patch loop    
 
@@ -200,6 +206,13 @@ subroutine canopy_structure( currentSite )
                          ! causing non-linearity issues with c_area.  is this really required? 
                          currentCohort%dbh = currentCohort%dbh 
                          copyc%dbh = copyc%dbh !+ 0.000000000001_r8
+
+                         ! keep track of number and biomass of demoted cohort
+                         currentSite%demotion_rate(currentCohort%size_class) = &
+                              currentSite%demotion_rate(currentCohort%size_class) + currentCohort%n
+                         currentSite%demotion_carbonflux = currentSite%demotion_carbonflux + &
+                              currentCohort%b * currentCohort%n
+
                          !kill the ones which go into canopy layers that are not allowed... (default nclmax=2) 
                          if(i+1 > nclmax)then 
                            !put the litter from the terminated cohorts into the fragmenting pools
@@ -246,8 +259,15 @@ subroutine canopy_structure( currentSite )
                          currentCohort%canopy_layer = i + 1 !the whole cohort becomes demoted
                          sumloss = sumloss + currentCohort%c_area 
 
+                         ! keep track of number and biomass of demoted cohort
+                         currentSite%demotion_rate(currentCohort%size_class) = &
+                              currentSite%demotion_rate(currentCohort%size_class) + currentCohort%n
+                         currentSite%demotion_carbonflux = currentSite%demotion_carbonflux + &
+                              currentCohort%b * currentCohort%n
+
                          !kill the ones which go into canopy layers that are not allowed... (default nclmax=2) 
                          if(i+1 > nclmax)then  
+
                            !put the litter from the terminated cohorts into the fragmenting pools
                             do c=1,ncwd
 
@@ -287,7 +307,7 @@ subroutine canopy_structure( currentSite )
 
                 enddo !currentCohort 
 
-                call terminate_cohorts(currentPatch)
+                call terminate_cohorts(currentSite, currentPatch)
                 arealayer(i) = arealayer(i) - sumloss
                 !Update arealayer for diff calculations of layer below. 
                 arealayer(i + 1) = arealayer(i + 1) + sumloss 
@@ -323,9 +343,9 @@ subroutine canopy_structure( currentSite )
 
        enddo !is there still excess area in any layer?      
 
-       call terminate_cohorts(currentPatch)
+       call terminate_cohorts(currentSite, currentPatch)
        call fuse_cohorts(currentPatch)
-       call terminate_cohorts(currentPatch)
+       call terminate_cohorts(currentSite, currentPatch)
 
        ! ----------- Check cohort area ------------------------------!
        do i = 1,z
@@ -368,6 +388,12 @@ subroutine canopy_structure( currentSite )
                          currentCohort%canopy_layer = i   
                          currentCohort%c_area = c_area(currentCohort)
 
+                         ! keep track of number and biomass of promoted cohort
+                         currentSite%promotion_rate(currentCohort%size_class) = &
+                              currentSite%promotion_rate(currentCohort%size_class) + currentCohort%n
+                         currentSite%promotion_carbonflux = currentSite%promotion_carbonflux + &
+                              currentCohort%b * currentCohort%n
+
                         ! write(fates_log(),*) 'promoting very small cohort', currentCohort%c_area,currentCohort%canopy_layer
                       endif
                       arealayer(currentCohort%canopy_layer) = arealayer(currentCohort%canopy_layer)+currentCohort%c_area 
@@ -427,12 +453,18 @@ subroutine canopy_structure( currentSite )
 
                          newarea = currentCohort%c_area - cc_gain !new area of existing cohort
                          copyc%n = currentCohort%n*cc_gain/currentCohort%c_area   !number of individuals in promoted cohort. 
-                         ! number of individuals in cohort remianing in understorey    
+                         ! number of individuals in cohort remaining in understorey    
                          currentCohort%n = currentCohort%n - (currentCohort%n*cc_gain/currentCohort%c_area) 
 
                          currentCohort%canopy_layer = i+1  !keep current cohort in the understory.        
                          copyc%canopy_layer = i ! promote copy to the higher canopy layer. 
 
+                         ! keep track of number and biomass of promoted cohort
+                         currentSite%promotion_rate(copyc%size_class) = &
+                              currentSite%promotion_rate(copyc%size_class) + copyc%n
+                         currentSite%promotion_carbonflux = currentSite%promotion_carbonflux + &
+                              copyc%b * copyc%n
+
                          ! seperate cohorts. 
                          ! needs to be a very small number to avoid causing non-linearity issues with c_area. 
                          ! is this really required? 
@@ -459,6 +491,12 @@ subroutine canopy_structure( currentSite )
                          ! if the upper canopy spread is smaller. this shold be dealt with by the 'excess area' loop.  
                          currentCohort%c_area = c_area(currentCohort) 
 
+                         ! keep track of number and biomass of promoted cohort
+                         currentSite%promotion_rate(currentCohort%size_class) = &
+                              currentSite%promotion_rate(currentCohort%size_class) + currentCohort%n
+                         currentSite%promotion_carbonflux = currentSite%promotion_carbonflux + &
+                              currentCohort%b * currentCohort%n
+
                          promswitch = 1
 
                         ! write(fates_log(),*) 'promoting whole cohort', currentCohort%c_area,cc_gain,currentCohort%canopy_layer, &
@@ -527,9 +565,9 @@ subroutine canopy_structure( currentSite )
           endif
        enddo !is there still not enough canopy area in any layer?         
 
-       call terminate_cohorts(currentPatch)
+       call terminate_cohorts(currentSite, currentPatch)
        call fuse_cohorts(currentPatch)
-       call terminate_cohorts(currentPatch)
+       call terminate_cohorts(currentSite, currentPatch)
 
        if(promswitch == 1)then
           !write(fates_log(),*) 'going into cohort check'

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

@@ -98,6 +98,7 @@ subroutine create_cohort(patchptr, pft, nn, hite, dbh, &
     new_cohort%dbh          = dbh
     new_cohort%canopy_trim  = ctrim
     new_cohort%canopy_layer = clayer
+    new_cohort%canopy_layer_yesterday = real(clayer, r8)
     new_cohort%laimemory    = laimemory
     new_cohort%bdead        = bdead
     new_cohort%balive       = balive
@@ -337,6 +338,7 @@ subroutine nan_cohort(cc_p)
     currentCohort%pft                = fates_unset_int  ! pft number                           
     currentCohort%indexnumber        = fates_unset_int  ! unique number for each cohort. (within clump?)
     currentCohort%canopy_layer       = fates_unset_int  ! canopy status of cohort (1 = canopy, 2 = understorey, etc.)   
+    currentCohort%canopy_layer_yesterday       = nan  ! recent canopy status of cohort (1 = canopy, 2 = understorey, etc.)   
     currentCohort%NV                 = fates_unset_int  ! Number of leaf layers: -
     currentCohort%status_coh         = fates_unset_int  ! growth status of plant  (2 = leaves on , 1 = leaves off)
     currentCohort%size_class         = fates_unset_int  ! size class index
@@ -483,7 +485,7 @@ subroutine zero_cohort(cc_p)
   end subroutine zero_cohort
 
   !-------------------------------------------------------------------------------------!
-  subroutine terminate_cohorts( patchptr )
+  subroutine terminate_cohorts( siteptr, patchptr )
     !
     ! !DESCRIPTION:
     ! terminates cohorts when they get too small      
@@ -493,6 +495,7 @@ subroutine terminate_cohorts( patchptr )
     use SFParamsMod, only : SF_val_CWD_frac
     !
     ! !ARGUMENTS    
+    type (ed_site_type), intent(inout), target :: siteptr
     type (ed_patch_type), intent(inout), target :: patchptr
     !
     ! !LOCAL VARIABLES:
@@ -501,6 +504,7 @@ subroutine terminate_cohorts( patchptr )
     type (ed_cohort_type) , pointer :: nextc
     integer :: terminate   ! do we terminate (1) or not (0) 
     integer :: c           ! counter for litter size class. 
+    integer :: levcan      ! canopy level
     !----------------------------------------------------------------------
 
     currentPatch  => patchptr
@@ -562,6 +566,17 @@ subroutine terminate_cohorts( patchptr )
        endif
 
        if (terminate == 1) then 
+          ! preserve a record of the to-be-terminated cohort for mortality accounting
+          if (currentCohort%canopy_layer .eq. 1) then
+             levcan = 1
+          else
+             levcan = 2
+          endif
+          siteptr%terminated_nindivs(currentCohort%size_class,currentCohort%pft,levcan) = &
+               siteptr%terminated_nindivs(currentCohort%size_class,currentCohort%pft,levcan) + currentCohort%n
+          !
+          siteptr%termination_carbonflux(levcan) = siteptr%termination_carbonflux(levcan) + &
+               currentCohort%n * currentCohort%b
           if (.not. associated(currentCohort%taller)) then
              currentPatch%tallest => currentCohort%shorter
           else 
@@ -757,6 +772,10 @@ subroutine fuse_cohorts(patchptr)
                          currentCohort%npp_bseed = (currentCohort%n*currentCohort%npp_bseed + nextc%n*nextc%npp_bseed)/newn
                          currentCohort%npp_store = (currentCohort%n*currentCohort%npp_store + nextc%n*nextc%npp_store)/newn
 
+                         ! recent canopy history
+                         currentCohort%canopy_layer_yesterday  = (currentCohort%n*currentCohort%canopy_layer_yesterday  + &
+                              nextc%n*nextc%canopy_layer_yesterday)/newn
+
                          do i=1, nlevcan     
                             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))
@@ -1034,6 +1053,7 @@ subroutine copy_cohort( currentCohort,copyc )
     n%gscan           = o%gscan
     n%leaf_cost       = o%leaf_cost
     n%canopy_layer    = o%canopy_layer
+    n%canopy_layer_yesterday    = o%canopy_layer_yesterday
     n%nv              = o%nv
     n%status_coh      = o%status_coh
     n%canopy_trim     = o%canopy_trim

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

@@ -273,6 +273,7 @@ subroutine spawn_patches( currentSite )
              !this is the case as the new patch probably doesn't have a closed canopy, and
              ! even if it does, that will be sorted out in canopy_structure. 
              nc%canopy_layer = 1 
+             nc%canopy_layer_yesterday = 1._r8 
 
              !mortality is dominant disturbance              
              if(currentPatch%disturbance_rates(1) > currentPatch%disturbance_rates(2))then 
@@ -395,7 +396,7 @@ subroutine spawn_patches( currentSite )
 
           !sort out the cohorts, since some of them may be so small as to need removing. 
           call fuse_cohorts(currentPatch)
-          call terminate_cohorts(currentPatch)
+          call terminate_cohorts(currentSite, currentPatch)
           call sort_cohorts(currentPatch)
 
           currentPatch => currentPatch%younger
@@ -412,7 +413,7 @@ subroutine spawn_patches( currentSite )
        currentSite%youngest_patch => new_patch
 
        call fuse_cohorts(new_patch)
-       call terminate_cohorts(new_patch)
+       call terminate_cohorts(currentSite, new_patch)
        call sort_cohorts(new_patch)
 
     endif !end new_patch area 

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

@@ -626,6 +626,7 @@ subroutine seeds_in( currentSite, cp_pnt )
     !
     ! !USES:
     use EDTypesMod, only : AREA
+    use EDTypesMod, only : homogenize_seed_pfts
     !
     ! !ARGUMENTS    
     type(ed_site_type), intent(inout), target  :: currentSite
@@ -635,12 +636,48 @@ subroutine seeds_in( currentSite, cp_pnt )
     type(ed_patch_type),  pointer :: currentPatch
     type(ed_cohort_type), pointer :: currentCohort
     integer :: p
+    logical :: pft_present(numpft_ed)
+    real(r8) :: npfts_present
     !----------------------------------------------------------------------
 
     currentPatch => cp_pnt
 
     currentPatch%seeds_in(:) = 0.0_r8
-
+
+    if ( homogenize_seed_pfts ) then
+       ! special mode to remove intergenerational filters on PFT existence: each PFT seeds all PFTs
+       ! first loop over all patches and cohorts to see what and how many PFTs are present on this site
+       pft_present(:) = .false.
+       npfts_present =  0._r8
+       currentPatch => currentSite%oldest_patch
+       do while(associated(currentPatch))
+          currentCohort => currentPatch%tallest
+          do while (associated(currentCohort))
+             p = currentCohort%pft
+             if (.not. pft_present(p)) then
+                pft_present(p) = .true.
+                npfts_present = npfts_present + 1._r8
+             endif
+             currentCohort => currentCohort%shorter
+          enddo !cohort loop                        
+          currentPatch => currentPatch%younger
+       enddo ! patch loop
+
+       ! now calculate the homogenized seed flux into each PFT pool
+       currentPatch => cp_pnt
+       currentCohort => currentPatch%tallest
+       do while (associated(currentCohort))
+          do p = 1, numpft_ed
+             if (pft_present(p)) then
+                currentPatch%seeds_in(p) = currentPatch%seeds_in(p) +  currentCohort%seed_prod * currentCohort%n / &
+                     (currentPatch%area * npfts_present)
+             endif
+          end do
+          currentCohort => currentCohort%shorter
+       enddo !cohort loop                  
+    else
+
+    ! normal case: each PFT seeds its own type
     currentCohort => currentPatch%tallest
     do while (associated(currentCohort))
        p = currentCohort%pft
@@ -649,6 +686,8 @@ subroutine seeds_in( currentSite, cp_pnt )
        currentCohort => currentCohort%shorter
     enddo !cohort loop
 
+    endif
+
     currentPatch => currentSite%oldest_patch
 
     do while(associated(currentPatch))
@@ -1040,6 +1079,9 @@ subroutine recruitment( t, currentSite, currentPatch )
            call create_cohort(currentPatch, temp_cohort%pft, temp_cohort%n, temp_cohort%hite, temp_cohort%dbh, &
                 temp_cohort%balive, temp_cohort%bdead, temp_cohort%bstore,  &
                 temp_cohort%laimemory, cohortstatus, temp_cohort%canopy_trim, currentPatch%NCL_p)
+
+           ! keep track of how many individuals were recruited for passing to history
+           currentSite%recruitment_rate(ft) = currentSite%recruitment_rate(ft) + temp_cohort%n
        endif
 
     enddo  !pft loop

components/clm/src/ED/main/EDInitMod.F90

@@ -81,6 +81,17 @@ subroutine zero_site( site_in )
     site_in%fates_to_bgc_this_ts = 0.0_r8
     site_in%fates_to_bgc_last_ts = 0.0_r8
 
+    ! termination and recruitment info
+    site_in%terminated_nindivs(:,:,:) = 0._r8
+    site_in%termination_carbonflux(:) = 0._r8
+    site_in%recruitment_rate(:) = 0._r8
+
+    ! demotion/promotion info
+    site_in%demotion_rate(:) = 0._r8
+    site_in%demotion_carbonflux = 0._r8
+    site_in%promotion_rate(:) = 0._r8
+    site_in%promotion_carbonflux = 0._r8
+
   end subroutine zero_site
 
   ! ============================================================================

components/clm/src/ED/main/EDMainMod.F90

@@ -118,7 +118,7 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in)
        call fuse_cohorts(currentPatch)            
 
        ! kills cohorts that are too small
-       call terminate_cohorts(currentPatch)
+       call terminate_cohorts(currentSite, currentPatch)
 
 
        currentPatch => currentPatch%younger
@@ -341,7 +341,7 @@ subroutine ed_update_site( currentSite, bc_in )
     currentPatch => currentSite%oldest_patch
     do while(associated(currentPatch))
 
-       call terminate_cohorts(currentPatch) 
+       call terminate_cohorts(currentSite, currentPatch) 
 
        ! FIX(SPM,040314) why is this needed for BFB restarts? Look into this at some point
        cohort_number = count_cohorts(currentPatch)