From e975b2e7785fe17d8dd341ef7d61a838a3fb5f91 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 4 Jan 2017 12:00:45 -0800 Subject: [PATCH 1/6] Potential fix to the negative leaf npp accounting problem. Still working on some diagnostics. Diagnostics are generating a floating point error. --- biogeochem/EDPhysiologyMod.F90 | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index fccd8c0843..2dd840be99 100755 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -839,12 +839,12 @@ subroutine Growth_Derivatives( currentSite, currentCohort) 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 @@ -860,19 +860,27 @@ subroutine Growth_Derivatives( currentSite, currentCohort) !what fraction of allocation do we divert to storage? !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(iulog,*) '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%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 write(iulog,*) 'ED: no leaf area in gd', currentCohort%indexnumber,currentCohort%n,currentCohort%bdead, & currentCohort%dbh,currentCohort%balive @@ -969,7 +977,6 @@ subroutine Growth_Derivatives( currentSite, currentCohort) 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) From 8de1e77630622ce0cc466f7328274864bdfe56c5 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 4 Jan 2017 13:30:47 -0800 Subject: [PATCH 2/6] CHanged the carbon accounting to when bleaf(cohort)<=0 to a failure. Added some print statements (temporarily). --- biogeochem/EDPhysiologyMod.F90 | 26 ++++++++++++++++++-------- main/ChecksBalancesMod.F90 | 16 ++++++++++++++++ 2 files changed, 34 insertions(+), 8 deletions(-) diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 2dd840be99..def99add7a 100755 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -20,6 +20,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 @@ -39,7 +45,8 @@ module EDPhysiologyMod public :: flux_into_litter_pools logical, parameter :: DEBUG = .false. ! local debug flag - + character(len=*), parameter, private :: sourcefile = & + __FILE__ ! ============================================================================ contains @@ -877,13 +884,16 @@ subroutine Growth_Derivatives( currentSite, currentCohort) else - 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 - write(iulog,*) 'ED: no leaf area in gd', currentCohort%indexnumber,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__)) + +! 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 +! write(iulog,*) 'ED: no leaf area in gd', currentCohort%indexnumber,currentCohort%n,currentCohort%bdead, & +! currentCohort%dbh,currentCohort%balive endif diff --git a/main/ChecksBalancesMod.F90 b/main/ChecksBalancesMod.F90 index 91d8d844fd..fa813b391e 100644 --- a/main/ChecksBalancesMod.F90 +++ b/main/ChecksBalancesMod.F90 @@ -119,6 +119,14 @@ subroutine SummarizeNetFluxes( nsites, sites, bc_in, is_beg_day ) sum(currentPatch%seed_decay) + sum(currentPatch%leaf_litter_out) + & sum(currentPatch%root_litter_out)) * & ( currentPatch%area/AREA ) * 1.e3_r8 / ( 365.0_r8*SHR_CONST_CDAY ) + print*,"TS TERMS:", sites(s)%fates_to_bgc_this_ts, & + sum(currentPatch%CWD_AG_out), & + sum(currentPatch%CWD_BG_out), & + sum(currentPatch%seed_decay), & + sum(currentPatch%leaf_litter_out), & + sum(currentPatch%root_litter_out), & + ( currentPatch%area/AREA ) * 1.e3_r8 / ( 365.0_r8*SHR_CONST_CDAY ) + ! sites(s)%tot_seed_rain_flux = sites(s)%tot_seed_rain_flux + & sum(sites(s)%seed_rain_flux) * 1.e3_r8 / ( 365.0_r8*SHR_CONST_CDAY ) @@ -206,6 +214,14 @@ 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 + print*,"ERR TERMS:",sites(s)%totfatesc, & + sites(s)%totfatesc_old, & + sites(s)%npp_timeintegrated, & + sites(s)%tot_seed_rain_flux*SHR_CONST_CDAY, & + sites(s)%fates_to_bgc_this_ts*SHR_CONST_CDAY, & + sites(s)%fire_c_to_atm*SHR_CONST_CDAY, & + sites(s)%cbal_err_fates + sites(s)%cbal_err_bgc = sites(s)%totbgcc - & sites(s)%totbgcc_old - & From 4869c73d947eba12095959d19986c77b90a94843 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 6 Jan 2017 12:46:13 -0800 Subject: [PATCH 3/6] Bug fix: in copy cohort, npp_acc was not being copied! Now it is. This was leading to NaNs later on in runs and causing lots of problems. I added some numerical inquiry intrinsics to catch bad values, but I am considering putting them in a different place. --- biogeochem/EDCohortDynamicsMod.F90 | 4 ++- biogeochem/EDPhysiologyMod.F90 | 9 +----- biogeophys/EDAccumulateFluxesMod.F90 | 46 ++++++++++++++++++++++++++-- main/ChecksBalancesMod.F90 | 17 ---------- main/FatesConstantsMod.F90 | 9 ++++++ 5 files changed, 56 insertions(+), 29 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 81143bd553..2ed5d80408 100755 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -1038,13 +1038,15 @@ subroutine copy_cohort( currentCohort,copyc ) n%gpp_acc_hold = o%gpp_acc_hold n%gpp_acc = o%gpp_acc n%gpp_tstep = o%gpp_tstep + 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_hold = o%npp_acc_hold + n%resp_tstep = o%resp_tstep n%resp_acc = o%resp_acc n%resp_acc_hold = o%resp_acc_hold diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index def99add7a..ee269e69cb 100755 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -873,6 +873,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort) !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 @@ -887,14 +888,6 @@ subroutine Growth_Derivatives( currentSite, currentCohort) write(fates_log(),*) 'No target leaf area in GrowthDerivs? Bleaf(cohort) <= 0?' call endrun(msg=errMsg(sourcefile, __LINE__)) -! 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 -! write(iulog,*) 'ED: no leaf area in gd', currentCohort%indexnumber,currentCohort%n,currentCohort%bdead, & -! currentCohort%dbh,currentCohort%balive - endif !Do we have enough carbon left over to make up the rest of the turnover demand? diff --git a/biogeophys/EDAccumulateFluxesMod.F90 b/biogeophys/EDAccumulateFluxesMod.F90 index 78563a3a2d..3310207f85 100644 --- a/biogeophys/EDAccumulateFluxesMod.F90 +++ b/biogeophys/EDAccumulateFluxesMod.F90 @@ -9,14 +9,19 @@ module EDAccumulateFluxesMod ! Rosie Fisher. March 2014. ! ! !USES: + use abortutils, only : endrun + use shr_log_mod , only : errMsg => shr_log_errMsg + use FatesConstantsMod, only : fates_huge, fates_tiny implicit none private ! public :: AccumulateFluxes_ED logical :: DEBUG = .false. ! for debugging this module - !------------------------------------------------------------------------------ + character(len=*), parameter, private :: sourcefile = & + __FILE__ + contains !------------------------------------------------------------------------------ @@ -33,6 +38,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 +74,45 @@ 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 + + ! Trap invalid values from photosynthesis and resp + ! ----------------------------------------------------------------------- + + if(ieee_is_nan(ccohort%gpp_tstep))then + write(iulog,*)'GPP NaN Trap Triggered',s,ifp,ccohort%gpp_tstep + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + if(ieee_is_nan(ccohort%resp_tstep))then + write(iulog,*)'RESP NaN Trap Triggered',s,ifp,ccohort%resp_tstep + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + if(ieee_is_nan(ccohort%npp_tstep))then + write(iulog,*)'NPP NaN Trap Triggered',s,ifp,ccohort%npp_tstep + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + if(.not.ieee_is_finite(ccohort%gpp_tstep))then + write(iulog,*)'GPP Infinite Trap Triggered',s,ifp,ccohort%gpp_tstep + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + if(.not.ieee_is_finite(ccohort%resp_tstep))then + write(iulog,*)'RESP Infinite Trap Triggered',s,ifp,ccohort%resp_tstep + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + if(.not.ieee_is_finite(ccohort%npp_tstep))then + write(iulog,*)'NPP Infinite Trap Triggered',s,ifp,ccohort%npp_tstep + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + 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 diff --git a/main/ChecksBalancesMod.F90 b/main/ChecksBalancesMod.F90 index fa813b391e..7a9e3e7467 100644 --- a/main/ChecksBalancesMod.F90 +++ b/main/ChecksBalancesMod.F90 @@ -119,14 +119,6 @@ subroutine SummarizeNetFluxes( nsites, sites, bc_in, is_beg_day ) sum(currentPatch%seed_decay) + sum(currentPatch%leaf_litter_out) + & sum(currentPatch%root_litter_out)) * & ( currentPatch%area/AREA ) * 1.e3_r8 / ( 365.0_r8*SHR_CONST_CDAY ) - print*,"TS TERMS:", sites(s)%fates_to_bgc_this_ts, & - sum(currentPatch%CWD_AG_out), & - sum(currentPatch%CWD_BG_out), & - sum(currentPatch%seed_decay), & - sum(currentPatch%leaf_litter_out), & - sum(currentPatch%root_litter_out), & - ( currentPatch%area/AREA ) * 1.e3_r8 / ( 365.0_r8*SHR_CONST_CDAY ) - ! sites(s)%tot_seed_rain_flux = sites(s)%tot_seed_rain_flux + & sum(sites(s)%seed_rain_flux) * 1.e3_r8 / ( 365.0_r8*SHR_CONST_CDAY ) @@ -214,15 +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 - print*,"ERR TERMS:",sites(s)%totfatesc, & - sites(s)%totfatesc_old, & - sites(s)%npp_timeintegrated, & - sites(s)%tot_seed_rain_flux*SHR_CONST_CDAY, & - sites(s)%fates_to_bgc_this_ts*SHR_CONST_CDAY, & - sites(s)%fire_c_to_atm*SHR_CONST_CDAY, & - sites(s)%cbal_err_fates - - sites(s)%cbal_err_bgc = sites(s)%totbgcc - & sites(s)%totbgcc_old - & (sites(s)%fates_to_bgc_last_ts*SHR_CONST_CDAY - & diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index 3df36d6b56..4368db1756 100644 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -55,4 +55,13 @@ 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) + + + + + end module FatesConstantsMod From 3ad68c9ac921c179d74c433ef5ba68fed26139e4 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 17 Jan 2017 11:49:54 -0800 Subject: [PATCH 4/6] Removal of ieee_ nan and inf checks. There was concern that many compilers have not adopted ieee_ checking functions and we would lose compiler compatibility. --- biogeophys/EDAccumulateFluxesMod.F90 | 35 ---------------------------- 1 file changed, 35 deletions(-) diff --git a/biogeophys/EDAccumulateFluxesMod.F90 b/biogeophys/EDAccumulateFluxesMod.F90 index 3310207f85..bd2437c92f 100644 --- a/biogeophys/EDAccumulateFluxesMod.F90 +++ b/biogeophys/EDAccumulateFluxesMod.F90 @@ -11,7 +11,6 @@ module EDAccumulateFluxesMod ! !USES: use abortutils, only : endrun use shr_log_mod , only : errMsg => shr_log_errMsg - use FatesConstantsMod, only : fates_huge, fates_tiny implicit none private ! @@ -79,40 +78,6 @@ subroutine AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) write(iulog,*) 'EDAccumFlux 67 ',ccohort%resp_tstep endif - ! Trap invalid values from photosynthesis and resp - ! ----------------------------------------------------------------------- - - if(ieee_is_nan(ccohort%gpp_tstep))then - write(iulog,*)'GPP NaN Trap Triggered',s,ifp,ccohort%gpp_tstep - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - if(ieee_is_nan(ccohort%resp_tstep))then - write(iulog,*)'RESP NaN Trap Triggered',s,ifp,ccohort%resp_tstep - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - if(ieee_is_nan(ccohort%npp_tstep))then - write(iulog,*)'NPP NaN Trap Triggered',s,ifp,ccohort%npp_tstep - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - if(.not.ieee_is_finite(ccohort%gpp_tstep))then - write(iulog,*)'GPP Infinite Trap Triggered',s,ifp,ccohort%gpp_tstep - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - if(.not.ieee_is_finite(ccohort%resp_tstep))then - write(iulog,*)'RESP Infinite Trap Triggered',s,ifp,ccohort%resp_tstep - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - if(.not.ieee_is_finite(ccohort%npp_tstep))then - write(iulog,*)'NPP Infinite Trap Triggered',s,ifp,ccohort%npp_tstep - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - 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 From d40d1807cd23346a5c7b91dc7ed351416a7b3798 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 19 Jan 2017 15:52:35 -0800 Subject: [PATCH 5/6] A fix was added to properly account for the npp flux to roots and sapwood in cases where there were no leaves. This case was not being triggered in typical testing, as there were no deciduous trees entering this condition. Its still not entirely clear why we are encountering trees that have no leaves yet are carbon positive, but that is not completely within the scope of this fix. This fix was tested using a multi-year f45 simulation and a special PFT set that jaholm developed that includes some temperate species, total of 5 pfts. This configuration also tripped a false positive warning message that occured while processing output boundary conditions in btran. It was assuming that transpiration uptake weighting was summing to unity after receiving conductance weighting from different pfts. I think previously we were getting unity just by luck of all pfts having the same root extinction parameters. The warning message complained and re-normalized. I simply removed the complaint, as the re-normalization was required for most cases. I also removed a pesky LHP log message by binding it in a debug. --- biogeochem/EDCanopyStructureMod.F90 | 2 +- biogeochem/EDCohortDynamicsMod.F90 | 77 +++++++++++++++-------------- biogeophys/EDBtranMod.F90 | 9 ++-- 3 files changed, 46 insertions(+), 42 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index f5419cedd7..20ae7efe54 100755 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/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 diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 6588ceeeba..63f7ef845f 100755 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -184,6 +184,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 @@ -218,68 +221,68 @@ 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)/udata%deltat - end if + new_bl = currentcohort%balive*leaf_frac + + new_br = pftcon%froot_leaf(ft) * (currentcohort%balive + currentcohort%laimemory) * leaf_frac - currentcohort%bl = currentcohort%balive*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) / udata%deltat + 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 - + max(0._r8,new_br - currentcohort%br) / udata%deltat + + currentcohort%npp_bsw = max(0._r8,new_bsw - currentcohort%bsw)/udata%deltat + 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)/udata%deltat + max(0.0_r8,new_br-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_bsw = max(0.0_r8, new_bsw-currentcohort%bsw)/udata%deltat currentcohort%npp_bdead = currentCohort%dbdeaddt end if + currentcohort%bl = 0.0_r8 + currentcohort%br = new_bl + currentcohort%bsw = new_bsw + endif if (abs(currentcohort%balive -currentcohort%bl- currentcohort%br - currentcohort%bsw)>1e-12) then diff --git a/biogeophys/EDBtranMod.F90 b/biogeophys/EDBtranMod.F90 index 8ac4a51b36..8283a4d52d 100644 --- a/biogeophys/EDBtranMod.F90 +++ b/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 From 5c4e164d1c22a1c8616cc33497f51c21947bd78c Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 19 Jan 2017 17:22:13 -0800 Subject: [PATCH 6/6] Typo fix in setting updated root biomass during growth. --- biogeochem/EDCohortDynamicsMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 63f7ef845f..012926ab09 100755 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -280,7 +280,7 @@ subroutine allocate_live_biomass(cc_p,mode) end if currentcohort%bl = 0.0_r8 - currentcohort%br = new_bl + currentcohort%br = new_br currentcohort%bsw = new_bsw endif