Fixes to carbon accounting and lr2 machine files

Merge branch 'rgknox-leafnpp-diag-fix'

This pull request mostly deals with some bug fixes to carbon
accounting. @jenniferholm noticed that the diagnostics for leaf_npp
were periodically showing negative values. Through #164 we identified
that this occurred because daily carbon balances during the allocation
sequences were sometimes negative to to high respiration and low
gpp. The model was correctly using storage carbon to "pay" the
negative daily carbon balance and allow maintenance respiration to
occur, it was not however correctly diagnosing this flow. The
accounting of this process is a little complicated because storage can
be used to pay for maintenance respiration as well as maintenance
turnover demand.

During the process of verifying that carbon accounting errors were
low, I triggered spurious values of output variables that triggered
netcdf write errors, this lead to the identification that
cohort%npp_accum was not being properly copied during copying of
cohorts during patch fission.

A fix was needed to lawrencium machine files for its lr2 partition for
serial runs. While these changes are unrelated to carbon accounting,
they are trivial and simple, so I bundled them here.

Fixes: #164, #169, #168 and possibly #154

User interface changes?: no

Code review: requesting @jenniferholm and @serbinsh for evaluation in
their science algorithms.

Testing:
  rgknox:
    Test suite: lawrencium lr3 (baseline) and lawrencium-lr2 (non-baseline) edTest, Rapid Science Check tool (single site multi-decadal analysis)
    Test baseline: 5c5928f
    Test namelist changes: none
    Test answer changes:
    Test summary: all PASS

  andre:
    Test suite: ed - yellowstone gnu, intel, pgi
                     hobart nag
    Test baseline: 30f84d7
    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

biogeochem/EDCanopyStructureMod.F90

@@ -1067,7 +1067,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si)
                    currentCohort%c_area/currentPatch%total_canopy_area)
              currentPatch%layer_height_profile(L,ft,iv) = currentPatch%layer_height_profile(L,ft,iv) + (remainder * fleaf * &
                    currentCohort%c_area/currentPatch%total_canopy_area*(layer_top_hite+layer_bottom_hite)/2.0_r8)
-             write(fates_log(), *) 'LHP', currentPatch%layer_height_profile(L,ft,iv)
+             if ( DEBUG ) write(fates_log(), *) 'LHP', currentPatch%layer_height_profile(L,ft,iv)
              if(currentCohort%dbh <= 0._r8.or.currentCohort%n == 0._r8)then
                 write(fates_log(), *) 'ED: dbh or n is zero in clmedlink', currentCohort%dbh,currentCohort%n
              endif

biogeochem/EDCohortDynamicsMod.F90

@@ -186,6 +186,9 @@ subroutine allocate_live_biomass(cc_p,mode)
                               ! accounts for the fact that live biomass may decline in the off-season, 
                               ! making leaf_memory unrealistic.    
     real(r8)  :: ratio_balive ! ratio between root+shoot biomass now and root+shoot biomass when leaves fell off. 
+    real(r8)  :: new_bl
+    real(r8)  :: new_br
+    real(r8)  :: new_bsw
 
     integer   :: ft           ! functional type
     integer   :: leaves_off_switch
@@ -220,68 +223,69 @@ subroutine allocate_live_biomass(cc_p,mode)
     ! 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)/freq_day
-       end if
 
-       currentcohort%bl = currentcohort%balive*leaf_frac
+       new_bl = currentcohort%balive*leaf_frac
+
+       new_br = pftcon%froot_leaf(ft) * (currentcohort%balive + currentcohort%laimemory) * leaf_frac
+
+       new_bsw = EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(currentcohort%balive + &
+            currentcohort%laimemory)*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
+       if(mode==1)then 
 
+          currentcohort%npp_leaf = currentcohort%npp_leaf + &
+                max(0.0_r8,new_bl - currentcohort%bl) / freq_day
+
           currentcohort%npp_froot = currentcohort%npp_froot + &
-               max(0._r8,pftcon%froot_leaf(ft)*(currentcohort%balive+currentcohort%laimemory)*leaf_frac - currentcohort%br) / &
-               freq_day
+                max(0._r8,new_br - currentcohort%br) / freq_day
 
-          currentcohort%npp_bsw = max(0._r8,EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(currentcohort%balive + &
-                currentcohort%laimemory)*leaf_frac - currentcohort%bsw)/freq_day
+          currentcohort%npp_bsw =  max(0.0_r8, new_bsw - currentcohort%bsw)/freq_day
 
           currentcohort%npp_bdead =  currentCohort%dbdeaddt
 
        end if
+
+       currentcohort%bl = new_bl
+       currentcohort%br = new_br
+       currentcohort%bsw = new_bsw
 
-       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
 
+    else ! Leaves are off (leaves_off_switch==1)
 
-       currentcohort%bl  = 0.0_r8
+       !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
+
        ideal_balive      = currentcohort%laimemory * pftcon%froot_leaf(ft) +  &
             currentcohort%laimemory*  EDecophyscon%sapwood_ratio(ft) * currentcohort%hite
-       currentcohort%br  = pftcon%froot_leaf(ft) * (ideal_balive + currentcohort%laimemory) * leaf_frac
-       currentcohort%bsw = EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(ideal_balive + &
-            currentcohort%laimemory)*leaf_frac 
-
-       ratio_balive           = currentcohort%balive / ideal_balive
-       currentcohort%br       = currentcohort%br  * ratio_balive
-       currentcohort%bsw      = currentcohort%bsw * ratio_balive
+       ratio_balive      = currentcohort%balive / ideal_balive
+
+       new_br  = pftcon%froot_leaf(ft) * (ideal_balive + currentcohort%laimemory) * &
+             leaf_frac *  ratio_balive
+       new_bsw = EDecophyscon%sapwood_ratio(ft) * currentcohort%hite * &
+             (ideal_balive + currentcohort%laimemory) * leaf_frac * 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)/freq_day
+                max(0.0_r8,new_br-currentcohort%br)/freq_day
 
-          currentcohort%npp_bsw = max(0.0_r8,EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(ideal_balive + &
-                currentcohort%laimemory)*leaf_frac*ratio_balive - currentcohort%bsw)/freq_day
+          currentcohort%npp_bsw = max(0.0_r8, new_bsw-currentcohort%bsw)/freq_day
 
           currentcohort%npp_bdead =  currentCohort%dbdeaddt
 
        end if
 
+       currentcohort%bl  = 0.0_r8
+       currentcohort%br  = new_br
+       currentcohort%bsw = new_bsw
+
     endif
 
     if (abs(currentcohort%balive -currentcohort%bl- currentcohort%br - currentcohort%bsw)>1e-12) then
@@ -1041,12 +1045,11 @@ subroutine copy_cohort( currentCohort,copyc )
 
     n%npp_acc_hold    = o%npp_acc_hold
     n%npp_tstep       = o%npp_tstep
+    n%npp_acc         = o%npp_acc
 
     if ( DEBUG ) write(fates_log(),*) 'EDcohortDyn Ia ',o%npp_acc
     if ( DEBUG ) write(fates_log(),*) 'EDcohortDyn Ib ',o%resp_acc
 
-    n%npp_acc         = o%npp_acc
-
     n%resp_tstep      = o%resp_tstep
     n%resp_acc        = o%resp_acc
     n%resp_acc_hold   = o%resp_acc_hold

biogeochem/EDPhysiologyMod.F90

@@ -21,6 +21,12 @@ module EDPhysiologyMod
   use EDTypesMod          , only : ncwd, cp_nlevcan, numpft_ed, senes
   use EDTypesMod          , only : ed_site_type, ed_patch_type, ed_cohort_type
 
+  use shr_log_mod           , only : errMsg => shr_log_errMsg
+  use abortutils            , only : endrun
+  use FatesGlobals          , only : fates_log
+
+
+
   implicit none
   private
 
@@ -40,7 +46,8 @@ module EDPhysiologyMod
   public :: flux_into_litter_pools
 
   logical, parameter :: DEBUG  = .false. ! local debug flag
-
+  character(len=*), parameter, private :: sourcefile = &
+        __FILE__
   ! ============================================================================
 
 contains
@@ -817,12 +824,12 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in)
     currentCohort%carbon_balance = currentCohort%npp_acc_hold - &
           currentCohort%md *  EDecophyscon%leaf_stor_priority(currentCohort%pft)
 
-    ! Allowing only carbon from NPP pool to account for npp flux into the maintenance pools
+    ! Allowing only carbon from NPP pool to account for npp flux into the maintenance turnover pools
     ! ie this does not include any use of storage carbon or balive to make up for missing carbon balance in the transfer
-    currentCohort%npp_leaf  = min(currentCohort%npp_acc_hold*currentCohort%leaf_md/currentCohort%md, &
-                                  currentCohort%leaf_md*EDecophyscon%leaf_stor_priority(currentCohort%pft))
-    currentCohort%npp_froot = min(currentCohort%npp_acc_hold*currentCohort%root_md/currentCohort%md, &
-                                  currentCohort%root_md*EDecophyscon%leaf_stor_priority(currentCohort%pft))
+    currentCohort%npp_leaf  = max(0.0_r8,min(currentCohort%npp_acc_hold*currentCohort%leaf_md/currentCohort%md, &
+                                  currentCohort%leaf_md*EDecophyscon%leaf_stor_priority(currentCohort%pft)))
+    currentCohort%npp_froot = max(0.0_r8,min(currentCohort%npp_acc_hold*currentCohort%root_md/currentCohort%md, &
+                                  currentCohort%root_md*EDecophyscon%leaf_stor_priority(currentCohort%pft)))
 
 
     if (Bleaf(currentCohort) > 0._r8)then
@@ -839,21 +846,26 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in)
           !what is the flux into the store?
           currentCohort%storage_flux = currentCohort%carbon_balance * f_store                     
 
+          currentCohort%npp_store = currentCohort%carbon_balance * f_store         
           if ( DEBUG ) write(fates_log(),*) 'EDphys B ',f_store
 
           !what is the tax on the carbon available for growth? 
           currentCohort%carbon_balance = currentCohort%carbon_balance * (1.0_r8 - f_store)  
        else  !cbalance is negative. Take C out of store to pay for maintenance respn.
+
           currentCohort%storage_flux = currentCohort%carbon_balance 
+
+          ! Note that npp_store only tracks the flux between NPP and storage.  Storage can 
+          ! also be drawn down to support some turnover demand.
+          currentCohort%npp_store = min(0.0_r8,currentCohort%npp_acc_hold)
+
           currentCohort%carbon_balance = 0._r8 
        endif
 
     else
 
-       currentCohort%storage_flux = 0._r8
-       currentCohort%carbon_balance = 0._r8
-       write(fates_log(),*) 'ED: no leaf area in gd',currentCohort%n,currentCohort%bdead, &
-             currentCohort%dbh,currentCohort%balive
+       write(fates_log(),*) 'No target leaf area in GrowthDerivs? Bleaf(cohort) <= 0?'
+       call endrun(msg=errMsg(sourcefile, __LINE__))
 
     endif
 
@@ -947,7 +959,6 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in)
     endif
 
     currentCohort%npp_bseed = currentCohort%seed_prod
-    currentCohort%npp_store = max(0.0_r8,currentCohort%storage_flux)
 
     ! calculate change in diameter and height 
     currentCohort%ddbhdt = currentCohort%dbdeaddt * dDbhdBd(currentCohort)

biogeophys/EDAccumulateFluxesMod.F90

@@ -9,14 +9,18 @@ module EDAccumulateFluxesMod
   ! Rosie Fisher. March 2014. 
   !
   ! !USES:
+  use abortutils, only        : endrun
+  use shr_log_mod , only      : errMsg => shr_log_errMsg
   implicit none
   private
   !
   public :: AccumulateFluxes_ED
 
   logical :: DEBUG = .false.  ! for debugging this module
-  !------------------------------------------------------------------------------
 
+  character(len=*), parameter, private :: sourcefile = &
+        __FILE__
+
 contains
 
   !------------------------------------------------------------------------------
@@ -33,6 +37,8 @@ subroutine AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time)
     use EDTypesMod        , only : ed_patch_type, ed_cohort_type, &
                                    ed_site_type, AREA
     use FatesInterfaceMod , only : bc_in_type,bc_out_type
+    use, intrinsic :: IEEE_ARITHMETIC
+
     !
     ! !ARGUMENTS    
     integer,            intent(in)            :: nsites
@@ -67,12 +73,11 @@ subroutine AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time)
                 ! _tstep fluxes are KgC/indiv/timestep _acc are KgC/indiv/day
 
                 if ( DEBUG ) then
-                   write(iulog,*) 'EDAccumFlux 64 ',ccohort%npp_acc, &
-                         ccohort%npp_tstep
+                   write(iulog,*) 'EDAccumFlux 64 ',ccohort%npp_tstep
                    write(iulog,*) 'EDAccumFlux 66 ',ccohort%gpp_tstep
                    write(iulog,*) 'EDAccumFlux 67 ',ccohort%resp_tstep
                 endif
-                
+
                 ccohort%npp_acc  = ccohort%npp_acc  + ccohort%npp_tstep 
                 ccohort%gpp_acc  = ccohort%gpp_acc  + ccohort%gpp_tstep 
                 ccohort%resp_acc = ccohort%resp_acc + ccohort%resp_tstep

biogeophys/EDBtranMod.F90

@@ -192,7 +192,7 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out)
                     end if
                  enddo
               enddo
-              
+
               !weight patch level output BTRAN for the
               bc_out(s)%btran_pa(ifp) = 0.0_r8
               do ft = 1,numpft_ed
@@ -203,10 +203,11 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out)
                     bc_out(s)%btran_pa(ifp)   = bc_out(s)%btran_pa(ifp) + cpatch%btran_ft(ft) * 1./numpft_ed
                  end if
               enddo
-
-              temprootr = sum(bc_out(s)%rootr_pagl(ifp,:))
+
+              ! While the in-pft root profiles summed to unity, averaging them weighted
+              ! by conductance, or not, will break sum to unity.  Thus, re-normalize.
+              temprootr = sum(bc_out(s)%rootr_pagl(ifp,1:cp_numlevgrnd))
               if(abs(1.0_r8-temprootr) > 1.0e-10_r8 .and. temprootr > 1.0e-10_r8)then
-                 write(iulog,*) 'error with rootr in canopy fluxes',temprootr,sum(pftgs),sum(cpatch%rootr_ft(1:2,:),dim=2)
                  do j = 1,cp_numlevgrnd
                     bc_out(s)%rootr_pagl(ifp,j) = bc_out(s)%rootr_pagl(ifp,j)/temprootr
                  enddo

main/ChecksBalancesMod.F90

@@ -206,7 +206,6 @@ subroutine FATES_BGC_Carbon_Balancecheck(nsites, sites, bc_in, is_beg_day, dtime
                   sites(s)%fire_c_to_atm*SHR_CONST_CDAY)
             sites(s)%cbal_err_fates = sites(s)%cbal_err_fates / SHR_CONST_CDAY
 
-
             sites(s)%cbal_err_bgc = sites(s)%totbgcc - &
                   sites(s)%totbgcc_old - &
                   (sites(s)%fates_to_bgc_last_ts*SHR_CONST_CDAY - &

main/FatesConstantsMod.F90

@@ -58,10 +58,16 @@ module FatesConstantsMod
   real(fates_r8), parameter :: t_water_freeze_k_triple = 273.16_fates_r8      
 
 
+  ! For numerical inquiry
+  real(fates_r8), parameter :: fates_huge = huge(g_per_kg)
+
+  real(fates_r8), parameter :: fates_tiny = tiny(g_per_kg)
+
   ! Geometric Constants
 
   ! PI
   real(fates_r8), parameter :: pi_const = 3.14159265359_fates_r8
 
 
+
 end module FatesConstantsMod