From c4d96b58038c51a4f52cdd2e4f0d5050c0b3ab36 Mon Sep 17 00:00:00 2001 From: Marcos Longo Date: Tue, 3 May 2022 11:54:38 -0700 Subject: [PATCH] Revised the carbon allocation routine to ensure allocation to fine roots and leaves are truly proportional to demand. Description: 1. Bug fix in DailyPRTAllometricCarbon (parteh/PRTAllometricCarbonMod.F90). When allocating to different tissues, the code was subtracting allocation of tissues before calculating the amount for the next tissue, potentially under-allocating carbon to fine roots. 2. Minor code updates to simplify allocation calculations. It now uses a single allocation factor based on availability and total need, which is applied to all tissues. 3. Replaced a few if statements with select case, which simplifies adding other hypotheses in the future (and it's safer for selecting cases). This pull request addresses the bug discussed in #784 and supersedes pull request #800. Additional minor changes in former PR #800 will be added as subsequent pull requests. --- biogeochem/FatesSoilBGCFluxMod.F90 | 5 +- biogeophys/FatesPlantRespPhotosynthMod.F90 | 1 + main/EDPftvarcon.F90 | 11 +- parteh/PRTAllometricCarbonMod.F90 | 150 +++++++++++---------- parteh/PRTParamsFATESMod.F90 | 43 +++--- 5 files changed, 113 insertions(+), 97 deletions(-) diff --git a/biogeochem/FatesSoilBGCFluxMod.F90 b/biogeochem/FatesSoilBGCFluxMod.F90 index d14ce7b005..cebab25728 100644 --- a/biogeochem/FatesSoilBGCFluxMod.F90 +++ b/biogeochem/FatesSoilBGCFluxMod.F90 @@ -219,7 +219,8 @@ subroutine UnPackNutrientAquisitionBCs(sites, bc_in) end do ! We can exit if this is a c-only simulation - if(hlm_parteh_mode.eq.prt_carbon_allom_hyp) then + select case (hlm_parteh_mode) + case (prt_carbon_allom_hyp) ! These can now be zero'd do s = 1, nsites bc_in(s)%plant_nh4_uptake_flux(:,:) = 0._r8 @@ -227,7 +228,7 @@ subroutine UnPackNutrientAquisitionBCs(sites, bc_in) bc_in(s)%plant_p_uptake_flux(:,:) = 0._r8 end do return - end if + end select do s = 1, nsites diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index adffe67d85..08607ec9f2 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -500,6 +500,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) end select + ! MLO - Shouldn't these numbers be parameters too? lmr25top = 2.525e-6_r8 * (1.5_r8 ** ((25._r8 - 20._r8)/10._r8)) lmr25top = lmr25top * lnc_top / (umolC_to_kgC * g_per_kg) diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index a149132a8d..b65cad5cc1 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -1502,7 +1502,8 @@ subroutine FatesCheckParams(is_master) if(.not.is_master) return - if (hlm_parteh_mode .eq. prt_cnp_flex_allom_hyp) then + select case (hlm_parteh_mode) + case (prt_cnp_flex_allom_hyp) ! Check to see if either RD/ECA/MIC is turned on @@ -1538,14 +1539,18 @@ subroutine FatesCheckParams(is_master) end if end if - elseif (hlm_parteh_mode .ne. prt_carbon_allom_hyp) then + case (prt_carbon_allom_hyp) + ! No additional checks needed for now. + continue + + case default write(fates_log(),*) 'FATES Plant Allocation and Reactive Transport has' write(fates_log(),*) 'only 2 modules supported, allometric carbon and CNP.' write(fates_log(),*) 'fates_parteh_mode must be set to 1 or 2 in the namelist' write(fates_log(),*) 'Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) - end if + end select ! If any PFTs are specified as either prescribed N or P uptake ! then they all must be ! diff --git a/parteh/PRTAllometricCarbonMod.F90 b/parteh/PRTAllometricCarbonMod.F90 index 5bdf624502..814d3daa36 100644 --- a/parteh/PRTAllometricCarbonMod.F90 +++ b/parteh/PRTAllometricCarbonMod.F90 @@ -51,6 +51,9 @@ module PRTAllometricCarbonMod use PRTParametersMod , only : prt_params + use EDTypesMod , only : leaves_on + use EDTypesMod , only : leaves_off + implicit none private @@ -114,7 +117,7 @@ module PRTAllometricCarbonMod ! ------------------------------------------------------------------------------------- - type, public, extends(prt_vartypes) :: callom_prt_vartypes + type, public, extends(prt_vartypes) :: callom_prt_vartypes contains @@ -144,7 +147,7 @@ module PRTAllometricCarbonMod public :: InitPRTGlobalAllometricCarbon -contains + contains subroutine InitPRTGlobalAllometricCarbon() @@ -313,6 +316,9 @@ subroutine DailyPRTAllometricCarbon(this) real(r8) :: struct_below_target ! dead (structural) biomass below target amount [kgC] real(r8) :: total_below_target ! total biomass below the allometric target [kgC] + real(r8) :: allocation_factor ! allocation factor (relative to demand) to + ! reconstruct tissues + real(r8) :: flux_adj ! adjustment made to growth flux term to minimize error [kgC] real(r8) :: store_target_fraction ! ratio between storage and leaf biomass when on allometry [kgC] @@ -451,12 +457,13 @@ subroutine DailyPRTAllometricCarbon(this) call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, ipft, target_struct_c) ! Target leaf biomass according to allometry and trimming - if(leaf_status==2) then + select case (leaf_status) + case (leaves_on) call bleaf(dbh,ipft,canopy_trim,target_leaf_c) - else + case (leaves_off) target_leaf_c = 0._r8 - end if - + end select + ! Target fine-root biomass and deriv. according to allometry and trimming [kgC, kgC/cm] call bfineroot(dbh,ipft,canopy_trim,target_fnrt_c) @@ -466,43 +473,42 @@ subroutine DailyPRTAllometricCarbon(this) ! ----------------------------------------------------------------------------------- ! III. Prioritize some amount of carbon to replace leaf/root turnover - ! Make sure it isnt a negative payment, and either pay what is available + ! Make sure it isn't a negative payment, and either pay what is available ! or forcefully pay from storage. ! ----------------------------------------------------------------------------------- - if( prt_params%evergreen(ipft) ==1 ) then + if( prt_params%evergreen(ipft) == itrue ) then leaf_c_demand = max(0.0_r8, & prt_params%leaf_stor_priority(ipft)*sum(this%variables(leaf_c_id)%turnover(:))) else leaf_c_demand = 0.0_r8 end if - fnrt_c_demand = max(0.0_r8, & + fnrt_c_demand = max(0.0_r8, & prt_params%leaf_stor_priority(ipft)*this%variables(fnrt_c_id)%turnover(icd)) total_c_demand = leaf_c_demand + fnrt_c_demand - - if (total_c_demand> nearzero ) then + + if (total_c_demand > nearzero) then ! We pay this even if we don't have the carbon ! Just don't pay so much carbon that storage+carbon_balance can't pay for it + allocation_factor = max(0.0_r8,min(1.0_r8,(store_c+carbon_balance)/total_c_demand)) - leaf_c_flux = min(leaf_c_demand, & - max(0.0_r8,(store_c+carbon_balance)* & - (leaf_c_demand/total_c_demand))) - - ! Add carbon to the youngest age pool (i.e iexp_leaf = index 1) - carbon_balance = carbon_balance - leaf_c_flux - leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux - - ! If we are testing b4b, then we pay this even if we don't have the carbon - fnrt_c_flux = min(fnrt_c_demand, & - max(0.0_r8, (store_c+carbon_balance)* & - (fnrt_c_demand/total_c_demand))) + ! MLO. Edited the code to switch the order of operations. The previous code would + ! subtract leaf flux from carbon balance before estimating the fine root flux, + ! potentially allowing less fluxes to fine roots than possible. + leaf_c_flux = leaf_c_demand * allocation_factor + fnrt_c_flux = fnrt_c_demand * allocation_factor - carbon_balance = carbon_balance - fnrt_c_flux - fnrt_c = fnrt_c + fnrt_c_flux + ! Add carbon to the youngest age pool (i.e iexp_leaf = index 1) and fine roots + leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux + fnrt_c = fnrt_c + fnrt_c_flux + ! Remove fluxes from carbon balance. In case we may have drawn carbon from storage, + ! carbon_balance will become negative, in which case we will deplete carbon from + ! storage in the next step. + carbon_balance = carbon_balance - ( leaf_c_flux + fnrt_c_flux ) end if ! ----------------------------------------------------------------------------------- @@ -511,19 +517,24 @@ subroutine DailyPRTAllometricCarbon(this) ! ----------------------------------------------------------------------------------- if( carbon_balance < 0.0_r8 ) then - + + ! Store_c_flux will be negative, so store_c will be depleted store_c_flux = carbon_balance carbon_balance = carbon_balance - store_c_flux store_c = store_c + store_c_flux else - store_below_target = max(target_store_c - store_c,0.0_r8) + ! Accumulate some carbon in storage. If storage is completely depleted, aim to + ! increase storage, but not to replenish completely so we can still use some + ! carbon for growth. + store_below_target = max(0.0_r8,target_store_c - store_c) store_target_fraction = max(0.0_r8, store_c/target_store_c ) store_c_flux = min(store_below_target,carbon_balance * & max(exp(-1.*store_target_fraction**4._r8) - exp( -1.0_r8 ),0.0_r8)) + ! Move carbon from carbon balance to storage carbon_balance = carbon_balance - store_c_flux store_c = store_c + store_c_flux @@ -533,24 +544,29 @@ subroutine DailyPRTAllometricCarbon(this) ! V. If carbon is still available, prioritize some allocation to replace ! the rest of the leaf/fineroot deficit ! carbon balance is guaranteed to be >=0 beyond this point + ! MLO. Renamed demand with below target to make it consistent with the + ! definitions at the variable declaration part. ! ----------------------------------------------------------------------------------- - - leaf_c_demand = max(0.0_r8,(target_leaf_c - sum(leaf_c(1:nleafage)))) - fnrt_c_demand = max(0.0_r8,(target_fnrt_c - fnrt_c)) - total_c_demand = leaf_c_demand + fnrt_c_demand - - if( (carbon_balance > nearzero ) .and. (total_c_demand>nearzero)) then + leaf_below_target = max(0.0_r8,target_leaf_c - sum(leaf_c(1:nleafage))) + fnrt_below_target = max(0.0_r8,target_fnrt_c - fnrt_c) + + total_below_target = leaf_below_target + fnrt_below_target + + if ( (carbon_balance > nearzero) .and. (total_below_target > nearzero) ) then + ! Find fraction of carbon to be allocated to leaves and fine roots + allocation_factor = min(1.0_r8, carbon_balance / total_below_target) + + ! MLO. Edited the code to switch the order of operations. The previous code would + ! subtract leaf flux from carbon balance before estimating the fine root flux, + ! potentially allowing less fluxes to fine roots than possible. + leaf_c_flux = leaf_below_target * allocation_factor + fnrt_c_flux = fnrt_below_target * allocation_factor - leaf_c_flux = min(leaf_c_demand, & - carbon_balance*(leaf_c_demand/total_c_demand)) - carbon_balance = carbon_balance - leaf_c_flux leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux - - fnrt_c_flux = min(fnrt_c_demand, & - carbon_balance*(fnrt_c_demand/total_c_demand)) - carbon_balance = carbon_balance - fnrt_c_flux - fnrt_c = fnrt_c + fnrt_c_flux + fnrt_c = fnrt_c + fnrt_c_flux + + carbon_balance = carbon_balance - ( leaf_c_flux + fnrt_c_flux ) end if @@ -571,31 +587,22 @@ subroutine DailyPRTAllometricCarbon(this) sapw_below_target + store_below_target if ( total_below_target > nearzero ) then - - if( total_below_target > carbon_balance) then - leaf_c_flux = carbon_balance * leaf_below_target/total_below_target - fnrt_c_flux = carbon_balance * fnrt_below_target/total_below_target - sapw_c_flux = carbon_balance * sapw_below_target/total_below_target - store_c_flux = carbon_balance * store_below_target/total_below_target - else - leaf_c_flux = leaf_below_target - fnrt_c_flux = fnrt_below_target - sapw_c_flux = sapw_below_target - store_c_flux = store_below_target - end if - - carbon_balance = carbon_balance - leaf_c_flux - leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux - - carbon_balance = carbon_balance - fnrt_c_flux - fnrt_c = fnrt_c + fnrt_c_flux - - carbon_balance = carbon_balance - sapw_c_flux - sapw_c = sapw_c + sapw_c_flux - - carbon_balance = carbon_balance - store_c_flux - store_c = store_c + store_c_flux - + ! Find allocation factor based on available carbon and total demand to meet target. + allocation_factor = min(1.0_r8, carbon_balance / total_below_target) + + ! Find fluxes to individual pools + leaf_c_flux = leaf_below_target * allocation_factor + fnrt_c_flux = fnrt_below_target * allocation_factor + sapw_c_flux = sapw_below_target * allocation_factor + store_c_flux = store_below_target * allocation_factor + + leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux + fnrt_c = fnrt_c + fnrt_c_flux + sapw_c = sapw_c + sapw_c_flux + store_c = store_c + store_c_flux + + carbon_balance = carbon_balance - & + ( leaf_c_flux + fnrt_c_flux + sapw_c_flux + store_c_flux ) end if end if @@ -690,11 +697,12 @@ subroutine DailyPRTAllometricCarbon(this) c_pool(dbh_id) = dbh ! Only grow leaves if we are in a "leaf-on" status - if(leaf_status==2) then - c_mask(leaf_c_id) = grow_leaf - else - c_mask(leaf_c_id) = .false. - end if + select case (leaf_status) + case (leaves_on) + c_mask(leaf_c_id) = grow_leaf + case default + c_mask(leaf_c_id) = .false. + end select c_mask(fnrt_c_id) = grow_fnrt c_mask(sapw_c_id) = grow_sapw c_mask(store_c_id) = grow_store diff --git a/parteh/PRTParamsFATESMod.F90 b/parteh/PRTParamsFATESMod.F90 index 208ff848fb..c017b4e138 100644 --- a/parteh/PRTParamsFATESMod.F90 +++ b/parteh/PRTParamsFATESMod.F90 @@ -1014,12 +1014,12 @@ subroutine PRTCheckParams(is_master) ! Check to make sure the organ ids are valid if this is the ! cnp_flex_allom_hypothesis - if ((hlm_parteh_mode .eq. prt_cnp_flex_allom_hyp) .or. & - (hlm_parteh_mode .eq. prt_carbon_allom_hyp) ) then + select case (hlm_parteh_mode) + case (prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp) do io = 1,norgans if(prt_params%organ_id(io) == repro_organ) then - write(fates_log(),*) 'with flexible cnp or c-only alloc hypothesese' + write(fates_log(),*) 'with flexible cnp or c-only alloc hypotheses' write(fates_log(),*) 'reproductive tissues are a special case' write(fates_log(),*) 'and therefore should not be included in' write(fates_log(),*) 'the parameter file organ list' @@ -1027,16 +1027,16 @@ subroutine PRTCheckParams(is_master) write(fates_log(),*) 'Aborting' end if if(prt_params%organ_id(io) == store_organ) then - write(fates_log(),*) 'with flexible cnp or c-only alloc hypothesese' + write(fates_log(),*) 'with flexible cnp or c-only alloc hypotheses' write(fates_log(),*) 'storage is a special case' write(fates_log(),*) 'and therefore should not be included in' write(fates_log(),*) 'the parameter file organ list' write(fates_log(),*) 'fates_prt_organ_id: ',prt_params%organ_id(:) write(fates_log(),*) 'Aborting' end if - + end do - end if + end select pftloop: do ipft = 1,npft @@ -1128,9 +1128,8 @@ subroutine PRTCheckParams(is_master) ! should not be re-translocating mass upon turnover. ! Note to advanced users. Feel free to remove these checks... ! ------------------------------------------------------------------- - - if ((hlm_parteh_mode .eq. prt_cnp_flex_allom_hyp) .or. & - (hlm_parteh_mode .eq. prt_carbon_allom_hyp) ) then + select case (hlm_parteh_mode) + case (prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp) do i = 1,norgans io = prt_params%organ_id(i) @@ -1165,10 +1164,11 @@ subroutine PRTCheckParams(is_master) end if end do - end if - - if (hlm_parteh_mode .eq. prt_cnp_flex_allom_hyp) then - + end select + + select case (hlm_parteh_mode) + case (prt_cnp_flex_allom_hyp) + ! Make sure nutrient storage fractions are positive if( prt_params%nitr_store_ratio(ipft) < 0._r8 ) then write(fates_log(),*) 'With parteh allometric CNP hypothesis' @@ -1235,9 +1235,9 @@ subroutine PRTCheckParams(is_master) end if end do - - end if - + + end select + ! Growth respiration ! if (parteh_mode .eq. prt_carbon_allom_hyp) then @@ -1258,8 +1258,8 @@ subroutine PRTCheckParams(is_master) ! end if ! end if - if ((hlm_parteh_mode .eq. prt_cnp_flex_allom_hyp) .or. & - (hlm_parteh_mode .eq. prt_carbon_allom_hyp) ) then + select case (hlm_parteh_mode) + case (prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp) ! The first nitrogen stoichiometry is used in all cases if ( (any(prt_params%nitr_stoich_p1(ipft,:) < 0.0_r8)) .or. & (any(prt_params%nitr_stoich_p1(ipft,:) >= 1.0_r8))) then @@ -1269,9 +1269,10 @@ subroutine PRTCheckParams(is_master) write(fates_log(),*) ' Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) end if - end if + end select - if(hlm_parteh_mode .eq. prt_cnp_flex_allom_hyp) then + select case (hlm_parteh_mode) + case (prt_cnp_flex_allom_hyp) do i = 1,norgans if ( (prt_params%nitr_stoich_p1(ipft,i) < 0._r8) .or. & @@ -1309,7 +1310,7 @@ subroutine PRTCheckParams(is_master) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - end if + end select ! Check turnover time-scales