size and composition structured history file diagnostics

Merge branch 'jenniferholm/io/sclass'

This feature merge adds size and composition structured output in
history type files. The dimensions of these new variables are
(time,levscpf,lndgrid). Patch level, instead of grid level,
diagnostics would also be desirable, but are not yet
available. levscpf is a combined size class and pft index compressed
into one dimension. The index of levscpf maps onto scls_levscpf, which
indicates the diameter size class of the lower bound [cm], and
pft_levscpf which indicates the pft. New variables have the form:
ED_NPP_[VAR]_GD_SCPF, and contain attribute data providing units and
long names. Notable additions are NPP fluxes partitioned by pool,
mortality partitioned by type, growth increment, basal area, GPP and
total autotrophic respiration.

Note: Science changing results come in the form of a new filter
preventing newly recruited cohorts from fusing with non-new
recruits. Recruits that have not experienced 1 day of growth yet do
not have many initialized rates, and therefore should not be fused
with other non-recruits. This prevents contamination of diagnostics.

Test suite: ed, lawrencium-lr2; ed yellowstone intel, gnu, pgi
Test baseline: none (science changing results)
Test namelist changes: none
Test status: PASS on all except the f09_g16 and f19_g16 grids, which
were expected fails.

Code reviewed by rgknox and jenniferholm, with discussions by rfisher,
cdkoven and bandre.

README

@@ -111,3 +111,4 @@ components/clm/src/biogeophys -- Biogeophysics (Hydrology)
                                      # (NOTE: ./xmlchange RESUBMIT=10 to set RESUBMIT to number
                                      # #  of times to automatically resubmit -- 10 in this example)
 
+END

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

@@ -12,6 +12,7 @@ module EDCohortDynamicsMod
   use EDTypesMod            , only : ed_site_type, ed_patch_type, ed_cohort_type
   use EDTypesMod            , only : fusetol, nclmax
   use EDtypesMod            , only : ncwd, numcohortsperpatch, udata
+  use EDtypesMod            , only : sclass_ed,nlevsclass_ed,AREA
   !
   implicit none
   private
@@ -108,7 +109,7 @@ subroutine create_cohort(patchptr, pft, nn, hite, dbh, &
     endif
 
     ! Calculate live biomass allocation
-    call allocate_live_biomass(new_cohort)
+    call allocate_live_biomass(new_cohort,0)
 
     ! Assign canopy extent and depth
     new_cohort%c_area  = c_area(new_cohort)
@@ -134,6 +135,13 @@ subroutine create_cohort(patchptr, pft, nn, hite, dbh, &
        patchptr%shortest => new_cohort 
     endif
 
+    ! Recuits do not have mortality rates, nor have they moved any
+    ! carbon when they are created.  They will bias our statistics
+    ! until they have experienced a full day.  We need a newly recruited flag.
+    ! This flag will be set to false after it has experienced 
+    ! growth, disturbance and mortality.
+    new_cohort%isnew = .true.
+
     call insert_cohort(new_cohort, patchptr%tallest, patchptr%shortest, tnull, snull, &
          storebigcohort, storesmallcohort)
 
@@ -143,7 +151,7 @@ subroutine create_cohort(patchptr, pft, nn, hite, dbh, &
   end subroutine create_cohort
 
   !-------------------------------------------------------------------------------------!
-  subroutine allocate_live_biomass(cc_p)
+  subroutine allocate_live_biomass(cc_p,mode)
     !
     ! !DESCRIPTION:
     ! Divide alive biomass between leaf, root and sapwood parts. 
@@ -153,6 +161,7 @@ subroutine allocate_live_biomass(cc_p)
     !
     ! !ARGUMENTS    
     type (ed_cohort_type), intent(inout), target  :: cc_p ! current cohort pointer
+    integer              , intent(in)             :: mode
     !
     ! !LOCAL VARIABLES:
     type (ed_cohort_type), pointer :: currentCohort
@@ -170,20 +179,19 @@ subroutine allocate_live_biomass(cc_p)
     ft = currentcohort%pft
     leaf_frac = 1.0_r8/(1.0_r8 + EDecophyscon%sapwood_ratio(ft) * currentcohort%hite + pftcon%froot_leaf(ft))     
 
-    currentcohort%bl = currentcohort%balive*leaf_frac    
-    ratio_balive = 1.0_r8
+    !currentcohort%bl = currentcohort%balive*leaf_frac    
     !for deciduous trees, there are no leaves  
 
     if (pftcon%evergreen(ft) == 1) then
        currentcohort%laimemory = 0._r8
        currentcohort%status_coh    = 2      
     endif    
 
-    !diagnore the root and stem biomass from the functional balance hypothesis. This is used when the leaves are 
+    ! iagnore the root and stem biomass from the functional balance hypothesis. This is used when the leaves are 
     !fully on. 
-    currentcohort%br  = pftcon%froot_leaf(ft) * (currentcohort%balive + currentcohort%laimemory) * leaf_frac
-    currentcohort%bsw = EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(currentcohort%balive + &
-         currentcohort%laimemory)*leaf_frac 
+    !currentcohort%br  = pftcon%froot_leaf(ft) * (currentcohort%balive + currentcohort%laimemory) * leaf_frac
+    !currentcohort%bsw = EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(currentcohort%balive + &
+    !     currentcohort%laimemory)*leaf_frac 
 
     leaves_off_switch = 0
     if (currentcohort%status_coh == 1.and.pftcon%stress_decid(ft) == 1.and.currentcohort%siteptr%dstatus==1) then !no leaves 
@@ -193,12 +201,46 @@ subroutine allocate_live_biomass(cc_p)
       leaves_off_switch = 1 !cold decid
     endif
 
-    if (leaves_off_switch==1) then
-       !the purpose of this section is to figure out the root and stem biomass when the leaves are off
-       !at this point, we know the former leaf mass (laimemory) and the current alive mass
-       !because balive may decline in the off-season, we need to adjust the root and stem biomass that are predicted
-       !from the laimemory, for the fact that we now might not have enough live biomass to support the hypothesized root mass
-       !thus, we use 'ratio_balive' to adjust br and bsw. Apologies that this is so complicated! RF
+    ! Use different proportions if the leaves are on vs off
+    if(leaves_off_switch==0)then
+
+       ! Tracking npp/gpp diagnostics only occur after growth derivatives is called
+       if(mode==1)then
+          ! it will not be able to put out as many leaves as it had previous timestep
+          currentcohort%npp_leaf = currentcohort%npp_leaf + &
+                max(0.0_r8,currentcohort%balive*leaf_frac - currentcohort%bl)/udata%deltat
+       end if
+
+       currentcohort%bl = currentcohort%balive*leaf_frac
+
+       !diagnose the root and stem biomass from the functional balance hypothesis. This is used when the leaves are 
+       !fully on. 
+       if(mode==1)then
+
+          currentcohort%npp_froot = currentcohort%npp_froot + &
+                max(0._r8,pftcon%froot_leaf(ft)*(currentcohort%balive+currentcohort%laimemory)*leaf_frac - currentcohort%br)/udata%deltat
+
+          currentcohort%npp_bsw = max(0._r8,EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(currentcohort%balive + &
+                currentcohort%laimemory)*leaf_frac - currentcohort%bsw)/udata%deltat
+
+          currentcohort%npp_bdead =  currentCohort%dbdeaddt
+
+       end if
+
+       currentcohort%br  = pftcon%froot_leaf(ft) * (currentcohort%balive + currentcohort%laimemory) * leaf_frac
+       currentcohort%bsw = EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(currentcohort%balive + &
+            currentcohort%laimemory)*leaf_frac
+
+
+    else ! Leaves are on (leaves_off_switch==1)
+
+    !the purpose of this section is to figure out the root and stem biomass when the leaves are off
+    !at this point, we know the former leaf mass (laimemory) and the current alive mass
+    !because balive may decline in the off-season, we need to adjust the root and stem biomass that are predicted
+    !from the laimemory, for the fact that we now might not have enough live biomass to support the hypothesized root mass
+    !thus, we use 'ratio_balive' to adjust br and bsw. Apologies that this is so complicated! RF
+
+
        currentcohort%bl  = 0.0_r8
        ideal_balive      = currentcohort%laimemory * pftcon%froot_leaf(ft) +  &
             currentcohort%laimemory*  EDecophyscon%sapwood_ratio(ft) * currentcohort%hite
@@ -208,7 +250,21 @@ subroutine allocate_live_biomass(cc_p)
 
        ratio_balive           = currentcohort%balive / ideal_balive
        currentcohort%br       = currentcohort%br  * ratio_balive
-       currentcohort%bsw      = currentcohort%bsw * ratio_balive               
+       currentcohort%bsw      = currentcohort%bsw * ratio_balive
+
+       ! Diagnostics
+       if(mode==1)then
+
+          currentcohort%npp_froot = currentcohort%npp_froot + &
+                max(0.0_r8,pftcon%froot_leaf(ft)*(ideal_balive + &
+                currentcohort%laimemory)*leaf_frac*ratio_balive-currentcohort%br)/udata%deltat
+
+          currentcohort%npp_bsw = max(0.0_r8,EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(ideal_balive + &
+                currentcohort%laimemory)*leaf_frac*ratio_balive - currentcohort%bsw)/udata%deltat
+
+          currentcohort%npp_bdead =  currentCohort%dbdeaddt
+
+       end if
 
     endif
 
@@ -294,6 +350,14 @@ subroutine nan_cohort(cc_p)
     currentCohort%resp_clm           = nan ! RESP: kgC/indiv/timestep
     currentCohort%resp_acc           = nan ! RESP: kGC/cohort/day
 
+    currentCohort%npp_leaf  = nan
+    currentCohort%npp_froot = nan
+    currentCohort%npp_bsw   = nan
+    currentCohort%npp_bdead = nan
+    currentCohort%npp_bseed = nan
+    currentCohort%npp_store = nan
+
+
     !RESPIRATION
     currentCohort%rd                 = nan
     currentCohort%resp_m             = nan ! Maintenance respiration.  kGC/cohort/year
@@ -387,6 +451,13 @@ subroutine zero_cohort(cc_p)
     currentcohort%gscan              = 0._r8 
     currentcohort%treesai            = 0._r8  
 
+    !    currentCohort%npp_leaf  = 0._r8
+    !    currentCohort%npp_froot = 0._r8
+    !    currentCohort%npp_bsw   = 0._r8
+    !    currentCohort%npp_bdead = 0._r8
+    !    currentCohort%npp_bseed = 0._r8
+    !    currentCohort%npp_store = 0._r8
+
   end subroutine zero_cohort
 
   !-------------------------------------------------------------------------------------!
@@ -533,7 +604,9 @@ subroutine fuse_cohorts(patchptr)
     iterate = 1
     fusion_took_place = 0   
     currentPatch => patchptr
-    maxcohorts =  currentPatch%NCL_p * numCohortsPerPatch  
+   ! maxcohorts =  currentPatch%NCL_p * numCohortsPerPatch
+    maxcohorts = numCohortsPerPatch
+
     !---------------------------------------------------------------------!
     !  Keep doing this until nocohorts <= maxcohorts                         !
     !---------------------------------------------------------------------!
@@ -559,7 +632,17 @@ subroutine fuse_cohorts(patchptr)
                    if (currentCohort%pft == nextc%pft) then              
 
                       ! check cohorts in same c. layer. before fusing
-                      if (currentCohort%canopy_layer == nextc%canopy_layer) then
+
+                      if (currentCohort%canopy_layer == nextc%canopy_layer) then 
+
+		      	 ! Note: because newly recruited cohorts that have not experienced
+			 ! a day yet will have un-known flux quantities or change rates
+			 ! we don't want them fusing with non-new cohorts.  We allow them
+			 ! to fuse with other new cohorts to keep the total number of cohorts
+			 ! down.
+
+                         if( (.not.(currentCohort%isnew) .and. .not.(nextc%isnew) ) .or. & 
+			          ( currentCohort%isnew .and. nextc%isnew ) ) then
 
                          fusion_took_place = 1         
                          newn = currentCohort%n + nextc%n    ! sum individuals in both cohorts.     
@@ -613,10 +696,25 @@ subroutine fuse_cohorts(patchptr)
                          currentCohort%npp         = (currentCohort%n*currentCohort%npp         + nextc%n*nextc%npp)/newn
                          currentCohort%gpp         = (currentCohort%n*currentCohort%gpp         + nextc%n*nextc%gpp)/newn
                          currentCohort%canopy_trim = (currentCohort%n*currentCohort%canopy_trim + nextc%n*nextc%canopy_trim)/newn
-                         currentCohort%dmort       = (currentCohort%n*currentCohort%dmort       + nextc%n*nextc%dmort)/newn
+			 currentCohort%dmort       = (currentCohort%n*currentCohort%dmort       + nextc%n*nextc%dmort)/newn
                          currentCohort%fire_mort   = (currentCohort%n*currentCohort%fire_mort   + nextc%n*nextc%fire_mort)/newn
                          currentCohort%leaf_litter = (currentCohort%n*currentCohort%leaf_litter + nextc%n*nextc%leaf_litter)/newn
 
+                         ! mortality diagnostics
+                         currentCohort%cmort = (currentCohort%n*currentCohort%cmort + nextc%n*nextc%cmort)/newn
+                         currentCohort%hmort = (currentCohort%n*currentCohort%hmort + nextc%n*nextc%hmort)/newn
+                         currentCohort%bmort = (currentCohort%n*currentCohort%bmort + nextc%n*nextc%bmort)/newn
+                         currentCohort%imort = (currentCohort%n*currentCohort%imort + nextc%n*nextc%imort)/newn
+                         currentCohort%fmort = (currentCohort%n*currentCohort%fmort + nextc%n*nextc%fmort)/newn
+
+                         ! npp diagnostics
+                         currentCohort%npp_leaf  = (currentCohort%n*currentCohort%npp_leaf  + nextc%n*nextc%npp_leaf)/newn
+                         currentCohort%npp_froot = (currentCohort%n*currentCohort%npp_froot + nextc%n*nextc%npp_froot)/newn
+                         currentCohort%npp_bsw   = (currentCohort%n*currentCohort%npp_bsw   + nextc%n*nextc%npp_bsw)/newn
+                         currentCohort%npp_bdead = (currentCohort%n*currentCohort%npp_bdead + nextc%n*nextc%npp_bdead)/newn
+                         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
+
                          do i=1, nlevcan_ed     
                             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))
@@ -639,6 +737,8 @@ subroutine fuse_cohorts(patchptr)
                             deallocate(nextc)            
                          endif
 
+                         endif ! Not a recruit
+
                       endif !canopy layer
                    endif !pft
                 endif  !index no. 
@@ -917,6 +1017,13 @@ subroutine copy_cohort( currentCohort,copyc )
     n%year_net_uptake = o%year_net_uptake
     n%ts_net_uptake   = o%ts_net_uptake
 
+    n%npp_leaf       = o%npp_leaf
+    n%npp_froot      = o%npp_froot
+    n%npp_bsw        = o%npp_bsw
+    n%npp_bdead      = o%npp_bdead
+    n%npp_bseed      = o%npp_bseed
+    n%npp_store      = o%npp_store
+
     !RESPIRATION
     n%rd              = o%rd
     n%resp_m          = o%resp_m
@@ -943,6 +1050,16 @@ subroutine copy_cohort( currentCohort,copyc )
     n%c_area          = o%c_area
     n%woody_turnover  = o%woody_turnover
 
+    ! Mortality diagnostics
+    n%cmort = o%cmort
+    n%bmort = o%bmort
+    n%imort = o%imort
+    n%fmort = o%fmort
+    n%hmort = o%hmort
+
+    ! Flags
+    n%isnew = o%isnew
+
     ! VARIABLES NEEDED FOR INTEGRATION 
     n%dndt            = o%dndt
     n%dhdt            = o%dhdt

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

@@ -322,15 +322,19 @@ real(r8) function dDbhdBl( cohort_in )
     dblddbh = 1.56_r8*0.0419_r8*(cohort_in%dbh**0.56_r8)*(EDecophyscon%wood_density(cohort_in%pft)**0.55_r8)
     dblddbh = dblddbh*cohort_in%canopy_trim
 
-    dDbhdBl = 1.0_r8/dblddbh
+    if( cohort_in%dbh<EDecophyscon%max_dbh(cohort_in%pft) ) then
+        dDbhdBl = 1.0_r8/dblddbh
+    else
+        dDbhdBl = 1.0d15  ! At maximum size, the leaf biomass is saturated, dbl=0
+    endif
 
     return
 
   end function dDbhdBl
 
 ! ============================================================================
 
-  real(r8) function mortality_rates( cohort_in )
+  subroutine mortality_rates( cohort_in,cmort,hmort,bmort )
 
     ! ============================================================================
     !  Calculate mortality rates as a function of carbon storage       
@@ -339,34 +343,39 @@ real(r8) function mortality_rates( cohort_in )
     use EDParamsMod,  only : ED_val_stress_mort
 
     type (ed_cohort_type), intent(in) :: cohort_in
+    real(r8),intent(out) :: bmort ! background mortality : Fraction per year
+    real(r8),intent(out) :: cmort  ! carbon starvation mortality
+    real(r8),intent(out) :: hmort  ! hydraulic failure mortality
 
     real(r8) :: frac  ! relativised stored carbohydrate
-    real(r8) :: smort ! stress mortality     : Fraction per year
-    real(r8) :: bmort ! background mortality : Fraction per year
 
     ! 'Background' mortality (can vary as a function of density as in ED1.0 and ED2.0, but doesn't here for tractability) 
     bmort = 0.014_r8 
 
     ! Proxy for hydraulic failure induced mortality. 
-    smort = 0.0_r8
     if(cohort_in%patchptr%btran_ft(cohort_in%pft) <= 0.000001_r8)then 
-       smort = smort + ED_val_stress_mort 
-    endif
+       hmort = ED_val_stress_mort
+     else
+       hmort = 0.0_r8
+     endif 
 
     ! Carbon Starvation induced mortality.
     if ( cohort_in%dbh  >  0._r8 ) then
-       if(Bleaf(cohort_in) > 0._r8.and.cohort_in%bstore <= Bleaf(cohort_in))then
+       if(Bleaf(cohort_in) > 0._r8 .and. cohort_in%bstore <= Bleaf(cohort_in))then
           frac = cohort_in%bstore/(Bleaf(cohort_in))
-          smort = smort + max(0.0_r8,ED_val_stress_mort*(1.0_r8 - frac))
+          cmort = max(0.0_r8,ED_val_stress_mort*(1.0_r8 - frac))
+        else
+          cmort = 0.0_r8
        endif
+
     else
        write(iulog,*) 'dbh problem in mortality_rates', &
             cohort_in%dbh,cohort_in%pft,cohort_in%n,cohort_in%canopy_layer,cohort_in%indexnumber
     endif
 
-    mortality_rates = smort + bmort
+    !mortality_rates = bmort + hmort + cmort
 
-  end function mortality_rates
+ end subroutine mortality_rates
 
 ! ============================================================================