From 1ace5a192c264374c3adc9055f0f720c9fbcea69 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 11 Apr 2019 17:45:17 -0700 Subject: [PATCH 1/3] Moved structure reset code out of parteh and into edcohortdynamics --- biogeochem/EDCohortDynamicsMod.F90 | 92 ++++++++++++++++++++++++++++ biogeochem/FatesAllometryMod.F90 | 1 + main/EDMainMod.F90 | 17 +++++- parteh/PRTAllometricCarbonMod.F90 | 97 ++++++------------------------ 4 files changed, 127 insertions(+), 80 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index a83f303a39..97cb1cc218 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -14,6 +14,7 @@ module EDCohortDynamicsMod use FatesConstantsMod , only : itrue,ifalse use FatesConstantsMod , only : fates_unset_r8 use FatesConstantsMod , only : nearzero + use FatesConstantsMod , only : calloc_abs_error use FatesInterfaceMod , only : hlm_days_per_year use FatesInterfaceMod , only : nleafage use EDPftvarcon , only : EDPftvarcon_inst @@ -47,6 +48,10 @@ module EDCohortDynamicsMod use FatesSizeAgeTypeIndicesMod, only : sizetype_class_index use FatesAllometryMod , only : bleaf use FatesAllometryMod , only : bfineroot + use FatesAllometryMod , only : bsap_allom + use FatesAllometryMod , only : bagw_allom + use FatesAllometryMod , only : bbgw_allom + use FatesAllometryMod , only : bdead_allom use FatesAllometryMod , only : h_allom use FatesAllometryMod , only : carea_allom use FatesAllometryMod , only : ForceDBH @@ -93,6 +98,7 @@ module EDCohortDynamicsMod public :: count_cohorts public :: InitPRTCohort public :: UpdateCohortBioPhysRates + public :: EvaluateAndCorrectDBH logical, parameter :: debug = .false. ! local debug flag @@ -1706,4 +1712,90 @@ end subroutine UpdateCohortBioPhysRates ! ============================================================================ + + subroutine EvaluateAndCorrectDBH(currentCohort,delta_dbh,delta_hite) + + ! ----------------------------------------------------------------------------------- + ! If the current diameter of a plant is somehow less than what is allometrically + ! consistent with stuctural biomass (or, in the case of grasses, leaf biomass) + ! then correct (increase) the dbh to match that. + ! ----------------------------------------------------------------------------------- + + ! argument + type(ed_cohort_type),intent(inout) :: currentCohort + real(r8),intent(out) :: delta_dbh + real(r8),intent(out) :: delta_hite + + ! locals + real(r8) :: dbh + real(r8) :: canopy_trim + integer :: ipft + real(r8) :: sapw_area + real(r8) :: target_sapw_c + real(r8) :: target_agw_c + real(r8) :: target_bgw_c + real(r8) :: target_struct_c + real(r8) :: target_leaf_c + real(r8) :: struct_c + real(r8) :: hite_out + real(r8) :: leaf_c + + dbh = currentCohort%dbh + ipft = currentCohort%pft + canopy_trim = currentCohort%canopy_trim + + delta_dbh = 0._r8 + delta_hite = 0._r8 + + if( EDPftvarcon_inst%woody(ipft) == itrue) then + + struct_c = currentCohort%prt%GetState(struct_organ, all_carbon_elements) + + ! Target sapwood biomass according to allometry and trimming [kgC] + call bsap_allom(dbh,ipft,canopy_trim,sapw_area,target_sapw_c) + + ! Target total above ground biomass in woody/fibrous tissues [kgC] + call bagw_allom(dbh,ipft,target_agw_c) + + ! Target total below ground biomass in woody/fibrous tissues [kgC] + call bbgw_allom(dbh,ipft,target_bgw_c) + + ! Target total dead (structrual) biomass [kgC] + call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, ipft, target_struct_c) + + ! ------------------------------------------------------------------------------------ + ! If structure is larger than target, then we need to correct some integration errors + ! by slightly increasing dbh to match it. + ! For grasses, if leaf biomass is larger than target, then we reset dbh to match + ! ----------------------------------------------------------------------------------- + + if( (struct_c - target_struct_c ) > calloc_abs_error ) then + call ForceDBH( ipft, canopy_trim, dbh, hite_out, bdead=struct_c ) + delta_dbh = dbh - currentCohort%dbh + delta_hite = hite_out - currentCohort%hite + currentCohort%dbh = dbh + currentCohort%hite = hite_out + end if + + else + + ! This returns the sum of leaf carbon over all (age) bins + leaf_c = currentCohort%prt%GetState(leaf_organ, all_carbon_elements) + + ! Target leaf biomass according to allometry and trimming + call bleaf(dbh,ipft,canopy_trim,target_leaf_c) + + if( ( leaf_c - target_leaf_c ) > calloc_abs_error ) then + call ForceDBH( ipft, canopy_trim, dbh, hite_out, bl=leaf_c ) + delta_dbh = dbh - currentCohort%dbh + delta_hite = hite_out - currentCohort%hite + currentCohort%dbh = dbh + currentCohort%hite = hite_out + end if + + end if + return + end subroutine EvaluateAndCorrectDBH + + end module EDCohortDynamicsMod diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 3d9ac56359..42eae69cdc 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -2200,6 +2200,7 @@ end function decay_coeff_kn ! ===================================================================================== + subroutine ForceDBH( ipft, canopy_trim, d, h, bdead, bl ) ! ========================================================================= diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index b6e8c8b252..c515a5e678 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -24,6 +24,7 @@ module EDMainMod use EDCohortDynamicsMod , only : fuse_cohorts use EDCohortDynamicsMod , only : sort_cohorts use EDCohortDynamicsMod , only : count_cohorts + use EDCohortDynamicsMod , only : EvaluateAndCorrectDBH use EDPatchDynamicsMod , only : disturbance_rates use EDPatchDynamicsMod , only : fuse_patches use EDPatchDynamicsMod , only : spawn_patches @@ -270,7 +271,9 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) real(r8) :: dbh_old ! dbh of plant before daily PRT [cm] real(r8) :: hite_old ! height of plant before daily PRT [m] real(r8) :: leaf_c - + real(r8) :: delta_dbh ! correction for dbh + real(r8) :: delta_hite ! correction for hite + !----------------------------------------------------------------------- small_no = 0.0000000000_r8 ! Obviously, this is arbitrary. RF - changed to zero @@ -314,8 +317,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) ! Apply Plant Allocation and Reactive Transport ! ----------------------------------------------------------------------------- - hite_old = currentCohort%hite - dbh_old = currentCohort%dbh + ! ----------------------------------------------------------------------------- ! Identify the net carbon gain for this dynamics interval @@ -373,6 +375,15 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) currentPatch%canopy_layer_tlai, currentCohort%treelai,currentCohort%vcmax25top,7 ) + ! If the current diameter of a plant is somehow less than what is consistent + ! with what is allometrically consistent with the stuctural biomass, then + ! correct the dbh to match. + + call EvaluateAndCorrectDBH(currentCohort,delta_dbh,delta_hite) + + hite_old = currentCohort%hite + dbh_old = currentCohort%dbh + ! Conduct Growth (parteh) call currentCohort%prt%DailyPRT() call currentCohort%prt%CheckMassConservation(ft,5) diff --git a/parteh/PRTAllometricCarbonMod.F90 b/parteh/PRTAllometricCarbonMod.F90 index 3a443c9050..7e35811d3b 100644 --- a/parteh/PRTAllometricCarbonMod.F90 +++ b/parteh/PRTAllometricCarbonMod.F90 @@ -442,84 +442,27 @@ subroutine DailyPRTAllometricCarbon(this) ! ----------------------------------------------------------------------------------- ! II. Calculate target size of the biomass compartment for a given dbh. ! ----------------------------------------------------------------------------------- - - if( EDPftvarcon_inst%woody(ipft) == itrue) then - - ! Target sapwood biomass according to allometry and trimming [kgC] - call bsap_allom(dbh,ipft,canopy_trim,sapw_area,target_sapw_c) - - ! Target total above ground biomass in woody/fibrous tissues [kgC] - call bagw_allom(dbh,ipft,target_agw_c) - - ! Target total below ground biomass in woody/fibrous tissues [kgC] - call bbgw_allom(dbh,ipft,target_bgw_c) - - ! Target total dead (structrual) biomass [kgC] - call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, ipft, target_struct_c) - - ! ------------------------------------------------------------------------------------ - ! If structure is larger than target, then we need to correct some integration errors - ! by slightly increasing dbh to match it. - ! For grasses, if leaf biomass is larger than target, then we reset dbh to match - ! ----------------------------------------------------------------------------------- - - if( (struct_c - target_struct_c ) > calloc_abs_error ) then - - call ForceDBH( ipft, canopy_trim, dbh, hite_out, bdead=struct_c ) - - ! Set the structural target biomass to the current structural boimass [kgC] - target_struct_c = struct_c - - ! Target sapwood biomass according to allometry and trimming [kgC] - call bsap_allom(dbh,ipft,canopy_trim,sapw_area,target_sapw_c) - - end if - - - ! Target leaf biomass according to allometry and trimming - call bleaf(dbh,ipft,canopy_trim,target_leaf_c) - - ! Target fine-root biomass and deriv. according to allometry and trimming [kgC, kgC/cm] - call bfineroot(dbh,ipft,canopy_trim,target_fnrt_c) - - ! Target storage carbon [kgC,kgC/cm] - call bstore_allom(dbh,ipft,canopy_trim,target_store_c) - - else - - ! Target leaf biomass according to allometry and trimming - call bleaf(dbh,ipft,canopy_trim,target_leaf_c) - - - if( ( sum(leaf_c) - target_leaf_c ) > calloc_abs_error ) then - - call ForceDBH( ipft, canopy_trim, dbh, hite_out, bl=sum(leaf_c) ) - - target_leaf_c = sum(leaf_c) - - end if - - ! Target fine-root biomass and deriv. according to allometry and trimming [kgC, kgC/cm] - call bfineroot(dbh,ipft,canopy_trim,target_fnrt_c) - - ! Target sapwood biomass according to allometry and trimming [kgC] - call bsap_allom(dbh,ipft,canopy_trim,sapw_area,target_sapw_c) - - ! Target total above ground biomass in woody/fibrous tissues [kgC] - call bagw_allom(dbh,ipft,target_agw_c) - - ! Target total below ground biomass in woody/fibrous tissues [kgC] - call bbgw_allom(dbh,ipft,target_bgw_c) - - ! Target total dead (structrual) biomass and [kgC] - call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, ipft, target_struct_c) - - ! Target storage carbon [kgC] - call bstore_allom(dbh,ipft,canopy_trim,target_store_c) - - - end if + ! Target sapwood biomass according to allometry and trimming [kgC] + call bsap_allom(dbh,ipft,canopy_trim,sapw_area,target_sapw_c) + + ! Target total above ground biomass in woody/fibrous tissues [kgC] + call bagw_allom(dbh,ipft,target_agw_c) + + ! Target total below ground biomass in woody/fibrous tissues [kgC] + call bbgw_allom(dbh,ipft,target_bgw_c) + + ! Target total dead (structrual) biomass [kgC] + call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, ipft, target_struct_c) + + ! Target leaf biomass according to allometry and trimming + call bleaf(dbh,ipft,canopy_trim,target_leaf_c) + + ! Target fine-root biomass and deriv. according to allometry and trimming [kgC, kgC/cm] + call bfineroot(dbh,ipft,canopy_trim,target_fnrt_c) + + ! Target storage carbon [kgC,kgC/cm] + call bstore_allom(dbh,ipft,canopy_trim,target_store_c) ! ----------------------------------------------------------------------------------- From 95635c3329ac9e3c70354178130e7dedc50aea20 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 12 Apr 2019 11:09:22 -0700 Subject: [PATCH 2/3] Added site level input carbon tracking diagnostics to seed_rain --- biogeochem/EDPhysiologyMod.F90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index f312d6e22f..488b6ec0b6 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -892,6 +892,10 @@ subroutine seeds_in( currentSite, cp_pnt ) EDPftvarcon_inst%seed_rain(p) !KgC/m2/year currentSite%seed_rain_flux(p) = currentSite%seed_rain_flux(p) + & EDPftvarcon_inst%seed_rain(p) * currentPatch%area/AREA !KgC/m2/year + + currentSite%flux_in = currentSite%flux_in + & + EDPftvarcon_inst%seed_rain(p) * currentPatch%area * hlm_freq_day + enddo From 84ec021b02727be251059b75e5469f25e680397d Mon Sep 17 00:00:00 2001 From: ckoven Date: Mon, 22 Apr 2019 12:43:04 -0600 Subject: [PATCH 3/3] added ability to make identical cohorts for multiple PFTs from census file when PFT=0 in census file Conflicts: main/FatesInventoryInitMod.F90 --- main/FatesInventoryInitMod.F90 | 132 ++++++++++++++++++++------------- 1 file changed, 80 insertions(+), 52 deletions(-) diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index e074944ae7..4de1931b5c 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -805,6 +805,7 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & real(r8) :: b_dead real(r8) :: b_store real(r8) :: a_sapwood ! area of sapwood at reference height [m2] + integer :: i_pft, ncohorts_to_create character(len=128),parameter :: wr_fmt = & @@ -853,9 +854,9 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & call endrun(msg=errMsg(sourcefile, __LINE__)) end if - if (c_pft <= 0 ) then + if (c_pft < 0 ) then write(fates_log(), *) 'inventory pft: ',c_pft - write(fates_log(), *) 'The inventory produced a cohort with <=0 pft index' + write(fates_log(), *) 'The inventory produced a cohort with <0 pft index' call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -882,62 +883,89 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & write(fates_log(), *) 'The inventory produced a cohort with very large density /m2' call endrun(msg=errMsg(sourcefile, __LINE__)) end if - - allocate(temp_cohort) ! A temporary cohort is needed because we want to make - ! use of the allometry functions - ! Don't need to allocate leaf age classes (not used) - temp_cohort%pft = c_pft - temp_cohort%n = c_nplant * cpatch%area - temp_cohort%dbh = c_dbh - call h_allom(c_dbh,c_pft,temp_cohort%hite) - temp_cohort%canopy_trim = 1.0_r8 + if (c_pft .eq. 0 ) then + write(fates_log(), *) 'inventory pft: ',c_pft + write(fates_log(), *) 'SPECIAL CASE TRIGGERED: PFT == 0 and therefore this subroutine' + write(fates_log(), *) 'will assign a cohort with n = n_orig/numpft to every cohort in range 1 to numpft' + ncohorts_to_create = numpft + else + ncohorts_to_create = 1 + end if - ! Calculate total above-ground biomass from allometry + do i_pft = 1,ncohorts_to_create + allocate(temp_cohort) ! A temporary cohort is needed because we want to make + ! use of the allometry functions + ! Don't need to allocate leaf age classes (not used) - call bagw_allom(temp_cohort%dbh,c_pft,b_agw) - ! Calculate coarse root biomass from allometry - call bbgw_allom(temp_cohort%dbh,c_pft,b_bgw) - - ! Calculate the leaf biomass (calculates a maximum first, then applies canopy trim - ! and sla scaling factors) - call bleaf(temp_cohort%dbh,c_pft,temp_cohort%canopy_trim,b_leaf) - - ! Calculate fine root biomass - call bfineroot(temp_cohort%dbh,c_pft,temp_cohort%canopy_trim,b_fineroot) - - ! Calculate sapwood biomass - call bsap_allom(temp_cohort%dbh,c_pft,temp_cohort%canopy_trim, a_sapwood, b_sapwood) - - call bdead_allom( b_agw, b_bgw, b_sapwood, c_pft, b_dead ) - call bstore_allom(temp_cohort%dbh, c_pft, temp_cohort%canopy_trim, b_store) - - temp_cohort%laimemory = 0._r8 - cstatus = leaves_on - - if( EDPftvarcon_inst%season_decid(c_pft) == itrue .and. csite%is_cold ) then - temp_cohort%laimemory = b_leaf - b_leaf = 0._r8 - cstatus = leaves_off - endif - - if ( EDPftvarcon_inst%stress_decid(c_pft) == itrue .and. csite%is_drought ) then - temp_cohort%laimemory = b_leaf - b_leaf = 0._r8 - cstatus = leaves_off - endif - - ! Since spread is a canopy level calculation, we need to provide an initial guess here. - call create_cohort(csite, cpatch, c_pft, temp_cohort%n, temp_cohort%hite, temp_cohort%dbh, & - b_leaf, b_fineroot, b_sapwood, b_dead, b_store, & - temp_cohort%laimemory, cstatus, rstatus, temp_cohort%canopy_trim, & - 1, csite%spread, equal_leaf_aclass, bc_in) + if (c_pft .ne. 0 ) then + ! normal case: assign each cohort to its specified PFT + temp_cohort%pft = c_pft + else + ! special case, make an identical cohort for each PFT + temp_cohort%pft = i_pft + endif - - deallocate(temp_cohort) ! get rid of temporary cohort + temp_cohort%n = c_nplant * cpatch%area / real(ncohorts_to_create,r8) + temp_cohort%dbh = c_dbh + + call h_allom(c_dbh,temp_cohort%pft,temp_cohort%hite) + temp_cohort%canopy_trim = 1.0_r8 + + ! Calculate total above-ground biomass from allometry + + call bagw_allom(temp_cohort%dbh,temp_cohort%pft,b_agw) + ! Calculate coarse root biomass from allometry + call bbgw_allom(temp_cohort%dbh,temp_cohort%pft,b_bgw) + + ! Calculate the leaf biomass (calculates a maximum first, then applies canopy trim + ! and sla scaling factors) + call bleaf(temp_cohort%dbh,temp_cohort%pft,temp_cohort%canopy_trim,b_leaf) + + ! Calculate fine root biomass + call bfineroot(temp_cohort%dbh,temp_cohort%pft,temp_cohort%canopy_trim,b_fineroot) + + ! Calculate sapwood biomass + call bsap_allom(temp_cohort%dbh,temp_cohort%pft,temp_cohort%canopy_trim, a_sapwood, b_sapwood) + + call bdead_allom( b_agw, b_bgw, b_sapwood, temp_cohort%pft, b_dead ) + + call bstore_allom(temp_cohort%dbh, temp_cohort%pft, temp_cohort%canopy_trim, b_store) + + temp_cohort%laimemory = 0._r8 + cstatus = leaves_on + + if( EDPftvarcon_inst%season_decid(temp_cohort%pft) == itrue .and. csite%is_cold ) then + temp_cohort%laimemory = b_leaf + b_leaf = 0._r8 + cstatus = leaves_off + endif + + if ( EDPftvarcon_inst%stress_decid(temp_cohort%pft) == itrue .and. csite%is_drought ) then + temp_cohort%laimemory = b_leaf + b_leaf = 0._r8 + cstatus = leaves_off + endif + + ! Since spread is a canopy level calculation, we need to provide an initial guess here. + if( debug_inv) then + write(fates_log(),*) 'calling create_cohort: ', temp_cohort%pft, temp_cohort%n, & + temp_cohort%hite, temp_cohort%dbh, & + b_leaf, b_fineroot, b_sapwood, b_dead, b_store, & + temp_cohort%laimemory, cstatus, rstatus, temp_cohort%canopy_trim, & + 1, csite%spread + endif + + call create_cohort(csite, cpatch, temp_cohort%pft, temp_cohort%n, temp_cohort%hite, & + temp_cohort%dbh, b_leaf, b_fineroot, b_sapwood, b_dead, b_store, & + temp_cohort%laimemory, cstatus, rstatus, temp_cohort%canopy_trim, & + 1, csite%spread, equal_leaf_aclass, bc_in) + + deallocate(temp_cohort) ! get rid of temporary cohort + end do return - end subroutine set_inventory_edcohort_type1 + end subroutine set_inventory_edcohort_type1 end module FatesInventoryInitMod