Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Carbon allocation bug fix #863

Merged
merged 3 commits into from
Jun 15, 2022
Merged

Carbon allocation bug fix #863

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
c4d96b5
Revised the carbon allocation routine to ensure allocation to fine ro…
mpaiao May 3, 2022
687f826
Merge branch 'master' into mpaiao-pr-carbon-alloc
glemieux Jun 2, 2022
4ee18db
Merge remote-tracking branch 'ngeet/master' into mpaiao-pr-carbon-alloc
glemieux Jun 13, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions biogeochem/FatesSoilBGCFluxMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -220,15 +220,16 @@ 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
bc_in(s)%plant_no3_uptake_flux(:,:) = 0._r8
bc_in(s)%plant_p_uptake_flux(:,:) = 0._r8
end do
return
end if
end select

do s = 1, nsites

Expand Down
1 change: 1 addition & 0 deletions biogeophys/FatesPlantRespPhotosynthMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -505,6 +505,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)

Expand Down
11 changes: 8 additions & 3 deletions main/EDPftvarcon.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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 !
Expand Down
150 changes: 79 additions & 71 deletions parteh/PRTAllometricCarbonMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,9 @@ module PRTAllometricCarbonMod

use PRTParametersMod , only : prt_params

use EDTypesMod , only : leaves_on
use EDTypesMod , only : leaves_off

implicit none
private

Expand Down Expand Up @@ -114,7 +117,7 @@ module PRTAllometricCarbonMod
! -------------------------------------------------------------------------------------


type, public, extends(prt_vartypes) :: callom_prt_vartypes
type, public, extends(prt_vartypes) :: callom_prt_vartypes

contains

Expand Down Expand Up @@ -144,7 +147,7 @@ module PRTAllometricCarbonMod
public :: InitPRTGlobalAllometricCarbon


contains
contains


subroutine InitPRTGlobalAllometricCarbon()
Expand Down Expand Up @@ -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]

Expand Down Expand Up @@ -455,12 +461,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)

Expand All @@ -470,43 +477,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

! -----------------------------------------------------------------------------------
Expand All @@ -515,19 +521,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

Expand All @@ -537,24 +548,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

Expand All @@ -575,31 +591,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

Expand Down Expand Up @@ -685,11 +692,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(struct_c_id) = grow_struct
Expand Down
Loading