Revised the carbon allocation routine to ensure allocation to fine ro…

…ots 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 NGEET#784 and supersedes pull request NGEET#800.
Additional minor changes in former PR NGEET#800 will be added as subsequent pull requests.

biogeochem/FatesSoilBGCFluxMod.F90

@@ -219,15 +219,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
 

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)
 

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 !

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