From ad53cdea7e19c3e2eb9dacefecba29c96406c371 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 28 Jul 2022 22:03:16 -0400 Subject: [PATCH 01/10] Removed laimemory, sapwmemory and structmemory --- biogeochem/EDCohortDynamicsMod.F90 | 28 +----- biogeochem/EDPhysiologyMod.F90 | 138 +++++++++++------------------ main/EDInitMod.F90 | 12 +-- main/EDPftvarcon.F90 | 9 ++ main/EDTypesMod.F90 | 6 -- main/FatesInventoryInitMod.F90 | 21 ++--- main/FatesRestartInterfaceMod.F90 | 31 ------- 7 files changed, 67 insertions(+), 178 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 4948f68129..a7d203da5d 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -147,8 +147,7 @@ module EDCohortDynamicsMod subroutine create_cohort(currentSite, patchptr, pft, nn, hite, coage, dbh, & - prt, leafmemory, sapwmemory, structmemory, & - status, recruitstatus,ctrim, carea, clayer, spread, bc_in) + prt, status, recruitstatus,ctrim, carea, clayer, spread, bc_in) ! ! !DESCRIPTION: ! create new cohort @@ -181,12 +180,6 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, coage, dbh, & real(r8), intent(in) :: dbh ! dbh: cm class(prt_vartypes),target :: prt ! The allocated PARTEH ! object - real(r8), intent(in) :: leafmemory ! target leaf biomass- set from - ! previous year: kGC per indiv - real(r8), intent(in) :: sapwmemory ! target sapwood biomass- set from - ! previous year: kGC per indiv - real(r8), intent(in) :: structmemory ! target structural biomass- set from - ! previous year: kGC per indiv real(r8), intent(in) :: ctrim ! What is the fraction of the maximum ! leaf biomass that we are targeting? real(r8), intent(in) :: spread ! The community assembly effects how @@ -237,9 +230,6 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, coage, dbh, & new_cohort%canopy_trim = ctrim new_cohort%canopy_layer = clayer new_cohort%canopy_layer_yesterday = real(clayer, r8) - new_cohort%leafmemory = leafmemory - new_cohort%sapwmemory = sapwmemory - new_cohort%structmemory = structmemory ! This sets things like vcmax25top, that depend on the ! leaf age fractions (which are defined by PARTEH) @@ -544,9 +534,6 @@ subroutine nan_cohort(cc_p) currentCohort%dbh = nan ! 'diameter at breast height' in cm currentCohort%coage = nan ! age of the cohort in years currentCohort%hite = nan ! height: meters - currentCohort%leafmemory = nan ! target leaf biomass- set from previous year: kGC per indiv - currentCohort%sapwmemory = nan ! target sapwood biomass- set from previous year: kGC per indiv - currentCohort%structmemory = nan ! target structural biomass- set from previous year: kGC per indiv currentCohort%lai = nan ! leaf area index of cohort m2/m2 currentCohort%sai = nan ! stem area index of cohort m2/m2 currentCohort%g_sb_laweight = nan ! Total leaf conductance of cohort (stomata+blayer) weighted by leaf-area [m/s]*[m2] @@ -1198,7 +1185,6 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) write(fates_log(),*) 'Cohort I, Cohort II' write(fates_log(),*) 'n:',currentCohort%n,nextc%n write(fates_log(),*) 'isnew:',currentCohort%isnew,nextc%isnew - write(fates_log(),*) 'leafmemory:',currentCohort%leafmemory,nextc%leafmemory write(fates_log(),*) 'hite:',currentCohort%hite,nextc%hite write(fates_log(),*) 'coage:',currentCohort%coage,nextc%coage write(fates_log(),*) 'dbh:',currentCohort%dbh,nextc%dbh @@ -1236,15 +1222,6 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) ! ----------------------------------------------------------------- call UpdateCohortBioPhysRates(currentCohort) - currentCohort%leafmemory = (currentCohort%n*currentCohort%leafmemory & - + nextc%n*nextc%leafmemory)/newn - - currentCohort%sapwmemory = (currentCohort%n*currentCohort%sapwmemory & - + nextc%n*nextc%sapwmemory)/newn - - currentCohort%structmemory = (currentCohort%n*currentCohort%structmemory & - + nextc%n*nextc%structmemory)/newn - currentCohort%canopy_trim = (currentCohort%n*currentCohort%canopy_trim & + nextc%n*nextc%canopy_trim)/newn @@ -1810,9 +1787,6 @@ subroutine copy_cohort( currentCohort,copyc ) n%dbh = o%dbh n%coage = o%coage n%hite = o%hite - n%leafmemory = o%leafmemory - n%sapwmemory = o%sapwmemory - n%structmemory = o%structmemory n%lai = o%lai n%sai = o%sai n%g_sb_laweight = o%g_sb_laweight diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 9a43af3226..d3815241fc 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -425,8 +425,6 @@ subroutine trim_canopy( currentSite ) real(r8) :: initial_trim ! Initial trim real(r8) :: optimum_trim ! Optimum trim value - real(r8) :: initial_leafmem ! Initial leafmemory - real(r8) :: optimum_leafmem ! Optimum leafmemory !---------------------------------------------------------------------- @@ -446,15 +444,13 @@ subroutine trim_canopy( currentSite ) currentCohort => currentPatch%tallest do while (associated(currentCohort)) - ! Save off the incoming trim and leafmemory + ! Save off the incoming trim initial_trim = currentCohort%canopy_trim - initial_leafmem = currentCohort%leafmemory ! Add debug diagnstic output to determine which cohort if (debug) then write(fates_log(),*) 'Current cohort:', icohort write(fates_log(),*) 'Starting canopy trim:', initial_trim - write(fates_log(),*) 'Starting leafmemory:', currentCohort%leafmemory endif trimmed = .false. @@ -596,10 +592,6 @@ subroutine trim_canopy( currentSite ) if (currentCohort%hite > EDPftvarcon_inst%hgt_min(ipft)) then currentCohort%canopy_trim = currentCohort%canopy_trim - & EDPftvarcon_inst%trim_inc(ipft) - if (prt_params%evergreen(ipft) /= 1)then - currentCohort%leafmemory = currentCohort%leafmemory * & - (1.0_r8 - EDPftvarcon_inst%trim_inc(ipft)) - endif trimmed = .true. @@ -641,17 +633,11 @@ subroutine trim_canopy( currentSite ) ! optimum_trim = (nnu_clai_b(1,1) / cumulative_lai_cohort) * initial_trim - optimum_leafmem = (nnu_clai_b(1,1) / cumulative_lai_cohort) * initial_leafmem ! Determine if the optimum trim value makes sense. The smallest cohorts tend to have unrealistic fits. if (optimum_trim > 0. .and. optimum_trim < 1.) then currentCohort%canopy_trim = optimum_trim - ! If the cohort pft is not evergreen we reduce the leafmemory as well - if (prt_params%evergreen(ipft) /= 1) then - currentCohort%leafmemory = optimum_leafmem - endif - trimmed = .true. endif @@ -1077,7 +1063,11 @@ subroutine phenology_leafonoff(currentSite) real(r8) :: struct_c ! structural wood carbon [kg] real(r8) :: store_c ! storage carbon [kg] real(r8) :: store_c_transfer_frac ! Fraction of storage carbon used to flush leaves - real(r8) :: totalmemory ! total memory of carbon [kg] + real(r8) :: deficit_c ! Amount of C needed to get flushing pools "on-allometry" + real(r8) :: target_leaf_c + real(r8) :: target_sapw_c + real(r8) :: target_bgw_c, target_agw_c, target_struct_c + real(r8) :: sapw_area integer :: ipft real(r8), parameter :: leaf_drop_fraction = 1.0_r8 real(r8), parameter :: carbon_store_buffer = 0.10_r8 @@ -1096,10 +1086,10 @@ subroutine phenology_leafonoff(currentSite) if(debug) call currentCohort%prt%CheckMassConservation(ipft,0) - store_c = currentCohort%prt%GetState(store_organ, all_carbon_elements) - leaf_c = currentCohort%prt%GetState(leaf_organ, all_carbon_elements) - sapw_c = currentCohort%prt%GetState(sapw_organ, all_carbon_elements) - struct_c = currentCohort%prt%GetState(struct_organ, all_carbon_elements) + store_c = currentCohort%prt%GetState(store_organ, carbon12_element) + leaf_c = currentCohort%prt%GetState(leaf_organ, carbon12_element) + sapw_c = currentCohort%prt%GetState(sapw_organ, carbon12_element) + struct_c = currentCohort%prt%GetState(struct_organ, carbon12_element) stem_drop_fraction = EDPftvarcon_inst%phen_stem_drop_fraction(ipft) @@ -1113,18 +1103,24 @@ subroutine phenology_leafonoff(currentSite) currentCohort%status_coh = leaves_on ! Leaves are on, so change status to ! stop flow of carbon out of bstore. + call bleaf(currentCohort%dbh,currentCohort%pft,currentCohort%canopy_trim,target_leaf_c) + call bsap_allom(currentCohort%dbh,currentCohort%pft,currentCohort%canopy_trim,sapw_area,target_sapw_c) + call bbgw_allom(currentCohort%dbh,currentCohort%pft,target_bgw_c) + call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, currentCohort%pft, target_struct_c) + + if(prt_params%woody(ipft) == itrue)then + deficit_c = target_leaf_c + else + deficit_c = target_leaf_c + (target_sapw_c-sapw_c) + (target_struct_c-struct_c) + end if + if(store_c>nearzero) then - ! flush either the amount required from the leafmemory, or -most- of the storage pool + + ! flush either the amount to get to the target, or -most- of the storage pool ! RF: added a criterion to stop the entire store pool emptying and triggering termination mortality ! n.b. this might not be necessary if we adopted a more gradual approach to leaf flushing... - store_c_transfer_frac = min((EDPftvarcon_inst%phenflush_fraction(ipft)* & - currentCohort%leafmemory)/store_c,(1.0_r8-carbon_store_buffer)) - - if(prt_params%woody(ipft).ne.itrue)then - totalmemory=currentCohort%leafmemory+currentCohort%sapwmemory+currentCohort%structmemory - store_c_transfer_frac = min((EDPftvarcon_inst%phenflush_fraction(ipft)* & - totalmemory)/store_c, (1.0_r8-carbon_store_buffer)) - endif + store_c_transfer_frac = min((EDPftvarcon_inst%phenflush_fraction(ipft)*deficit_c)/store_c, & + (1.0_r8-carbon_store_buffer)) else store_c_transfer_frac = 0.0_r8 @@ -1135,21 +1131,21 @@ subroutine phenology_leafonoff(currentSite) if(prt_params%woody(ipft) == itrue) then call PRTPhenologyFlush(currentCohort%prt, ipft, leaf_organ, store_c_transfer_frac) - currentCohort%leafmemory = 0.0_r8 else - ! Check that the stem drop fraction is set to non-zero amount otherwise flush all carbon store to leaves + ! Check that the stem drop fraction is set to non-zero amount + ! otherwise flush all carbon store to leaves if (stem_drop_fraction .gt. 0.0_r8) then call PRTPhenologyFlush(currentCohort%prt, ipft, leaf_organ, & - store_c_transfer_frac*currentCohort%leafmemory/totalmemory) + store_c_transfer_frac*target_leaf_c/deficit_c) call PRTPhenologyFlush(currentCohort%prt, ipft, sapw_organ, & - store_c_transfer_frac*currentCohort%sapwmemory/totalmemory) + store_c_transfer_frac*(target_sapw_c-sapw_c)/deficit_c) call PRTPhenologyFlush(currentCohort%prt, ipft, struct_organ, & - store_c_transfer_frac*currentCohort%structmemory/totalmemory) + store_c_transfer_frac*(target_struct_c-struct_c)/deficit_c) else @@ -1158,10 +1154,6 @@ subroutine phenology_leafonoff(currentSite) end if - currentCohort%leafmemory = 0.0_r8 - currentCohort%structmemory = 0.0_r8 - currentCohort%sapwmemory = 0.0_r8 - endif endif !pft phenology endif ! growing season @@ -1179,11 +1171,6 @@ subroutine phenology_leafonoff(currentSite) ! This sets the cohort to the "leaves off" flag currentCohort%status_coh = leaves_off - ! Remember what the leaf mass was for next year - ! the same amount back on in the spring... - - currentCohort%leafmemory = leaf_c - ! Drop Leaves (this routine will update the leaf state variables, ! for carbon and any other element that are prognostic. It will ! also track the turnover masses that will be sent to litter later on) @@ -1193,10 +1180,6 @@ subroutine phenology_leafonoff(currentSite) if(prt_params%woody(ipft).ne.itrue)then - currentCohort%sapwmemory = sapw_c * stem_drop_fraction - - currentCohort%structmemory = struct_c * stem_drop_fraction - call PRTDeciduousTurnover(currentCohort%prt,ipft, & sapw_organ, stem_drop_fraction) @@ -1226,19 +1209,22 @@ subroutine phenology_leafonoff(currentSite) currentCohort%status_coh = leaves_on ! Leaves are on, so change status to ! stop flow of carbon out of bstore. - if(store_c>nearzero) then - - store_c_transfer_frac = & - min((EDPftvarcon_inst%phenflush_fraction(ipft)*currentCohort%leafmemory)/store_c, & - (1.0_r8-carbon_store_buffer)) - - if(prt_params%woody(ipft).ne.itrue)then - - totalmemory=currentCohort%leafmemory+currentCohort%sapwmemory+currentCohort%structmemory - store_c_transfer_frac = min(EDPftvarcon_inst%phenflush_fraction(ipft)*totalmemory/store_c, & - (1.0_r8-carbon_store_buffer)) + call bleaf(currentCohort%dbh,currentCohort%pft,currentCohort%canopy_trim,target_leaf_c) + call bsap_allom(currentCohort%dbh,currentCohort%pft,currentCohort%canopy_trim,sapw_area,target_sapw_c) + call bbgw_allom(currentCohort%dbh,currentCohort%pft,target_bgw_c) + call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, currentCohort%pft, target_struct_c) - endif + if(prt_params%woody(ipft) == itrue)then + deficit_c = target_leaf_c + else + deficit_c = target_leaf_c + (target_sapw_c-sapw_c) + (target_struct_c-struct_c) + end if + + if(store_c>nearzero) then + + store_c_transfer_frac = & + min((EDPftvarcon_inst%phenflush_fraction(ipft)*deficit_c)/store_c, & + (1.0_r8-carbon_store_buffer)) else store_c_transfer_frac = 0.0_r8 @@ -1251,22 +1237,20 @@ subroutine phenology_leafonoff(currentSite) call PRTPhenologyFlush(currentCohort%prt, ipft, & leaf_organ, store_c_transfer_frac) - currentCohort%leafmemory = 0.0_r8 - else ! Check that the stem drop fraction is set to non-zero amount otherwise flush all carbon store to leaves if (stem_drop_fraction .gt. 0.0_r8) then call PRTPhenologyFlush(currentCohort%prt, ipft, leaf_organ, & - store_c_transfer_frac*currentCohort%leafmemory/totalmemory) + store_c_transfer_frac*target_leaf_c/deficit_c) call PRTPhenologyFlush(currentCohort%prt, ipft, sapw_organ, & - store_c_transfer_frac*currentCohort%sapwmemory/totalmemory) + store_c_transfer_frac*(target_sapw_c-sapw_c)/deficit_c) call PRTPhenologyFlush(currentCohort%prt, ipft, struct_organ, & - store_c_transfer_frac*currentCohort%structmemory/totalmemory) - + store_c_transfer_frac*(target_struct_c-struct_c)/deficit_c) + else call PRTPhenologyFlush(currentCohort%prt, ipft, leaf_organ, & @@ -1274,10 +1258,6 @@ subroutine phenology_leafonoff(currentSite) end if - currentCohort%leafmemory = 0.0_r8 - currentCohort%structmemory = 0.0_r8 - currentCohort%sapwmemory = 0.0_r8 - endif ! woody plant check endif !currentCohort status again? endif !currentSite status @@ -1291,17 +1271,11 @@ subroutine phenology_leafonoff(currentSite) ! This sets the cohort to the "leaves off" flag currentCohort%status_coh = leaves_off - ! Remember what the leaf mass was for next year - currentCohort%leafmemory = leaf_c - call PRTDeciduousTurnover(currentCohort%prt,ipft, & leaf_organ, leaf_drop_fraction) if(prt_params%woody(ipft).ne.itrue)then - currentCohort%sapwmemory = sapw_c * stem_drop_fraction - currentCohort%structmemory = struct_c * stem_drop_fraction - call PRTDeciduousTurnover(currentCohort%prt,ipft, & sapw_organ, stem_drop_fraction) @@ -1881,22 +1855,15 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) ! Default assumption is that leaves are on cohortstatus = leaves_on - temp_cohort%leafmemory = 0.0_r8 - temp_cohort%sapwmemory = 0.0_r8 - temp_cohort%structmemory = 0.0_r8 - ! But if the plant is seasonally (cold) deciduous, and the site status is flagged ! as "cold", then set the cohort's status to leaves_off, and remember the leaf biomass if ((prt_params%season_decid(ft) == itrue) .and. & (any(currentSite%cstatus == [phen_cstat_nevercold,phen_cstat_iscold]))) then - temp_cohort%leafmemory = c_leaf c_leaf = 0.0_r8 ! If plant is not woody then set sapwood and structural biomass as well if (prt_params%woody(ft).ne.itrue) then - temp_cohort%sapwmemory = c_sapw * stem_drop_fraction - temp_cohort%structmemory = c_struct * stem_drop_fraction c_sapw = (1.0_r8 - stem_drop_fraction) * c_sapw c_struct = (1.0_r8 - stem_drop_fraction) * c_struct endif @@ -1908,13 +1875,10 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) ! biomass if ((prt_params%stress_decid(ft) == itrue) .and. & (any(currentSite%dstatus == [phen_dstat_timeoff,phen_dstat_moistoff]))) then - temp_cohort%leafmemory = c_leaf c_leaf = 0.0_r8 ! If plant is not woody then set sapwood and structural biomass as well if(prt_params%woody(ft).ne.itrue)then - temp_cohort%sapwmemory = c_sapw * stem_drop_fraction - temp_cohort%structmemory = c_struct * stem_drop_fraction c_sapw = (1.0_r8 - stem_drop_fraction) * c_sapw c_struct = (1.0_r8 - stem_drop_fraction) * c_struct endif @@ -1988,7 +1952,7 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) endif ! Only bother allocating a new cohort if there is a reasonable amount of it - any_recruits: if (temp_cohort%n > min_n_safemath )then + any_recruits: if (temp_cohort%n > min_n_safemath )then ! ----------------------------------------------------------------------------- ! PART II. @@ -2091,10 +2055,10 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) ! ----------------------------------------------------------------------------------- call prt%CheckInitialConditions() + ! This initializes the cohort call create_cohort(currentSite,currentPatch, temp_cohort%pft, temp_cohort%n, & temp_cohort%hite, temp_cohort%coage, temp_cohort%dbh, prt, & - temp_cohort%leafmemory, temp_cohort%sapwmemory, temp_cohort%structmemory, & cohortstatus, recruitstatus, & temp_cohort%canopy_trim,temp_cohort%c_area, & currentPatch%NCL_p, currentSite%spread, bc_in) diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 451834a9ec..688c07010e 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -794,9 +794,6 @@ subroutine init_cohorts( site_in, patch_in, bc_in) call bstore_allom(temp_cohort%dbh, pft, temp_cohort%canopy_trim, c_store) - temp_cohort%leafmemory = 0._r8 - temp_cohort%sapwmemory = 0._r8 - temp_cohort%structmemory = 0._r8 cstatus = leaves_on stem_drop_fraction = EDPftvarcon_inst%phen_stem_drop_fraction(temp_cohort%pft) @@ -805,9 +802,6 @@ subroutine init_cohorts( site_in, patch_in, bc_in) if( prt_params%season_decid(pft) == itrue .and. & any(site_in%cstatus == [phen_cstat_nevercold,phen_cstat_iscold])) then - temp_cohort%leafmemory = c_leaf - temp_cohort%sapwmemory = c_sapw * stem_drop_fraction - temp_cohort%structmemory = c_struct * stem_drop_fraction c_leaf = 0._r8 c_sapw = (1.0_r8-stem_drop_fraction) * c_sapw c_struct = (1.0_r8-stem_drop_fraction) * c_struct @@ -816,9 +810,6 @@ subroutine init_cohorts( site_in, patch_in, bc_in) if ( prt_params%stress_decid(pft) == itrue .and. & any(site_in%dstatus == [phen_dstat_timeoff,phen_dstat_moistoff])) then - temp_cohort%leafmemory = c_leaf - temp_cohort%sapwmemory = c_sapw * stem_drop_fraction - temp_cohort%structmemory = c_struct * stem_drop_fraction c_leaf = 0._r8 c_sapw = (1.0_r8-stem_drop_fraction) * c_sapw c_struct = (1.0_r8-stem_drop_fraction) * c_struct @@ -900,8 +891,7 @@ subroutine init_cohorts( site_in, patch_in, bc_in) call prt_obj%CheckInitialConditions() call create_cohort(site_in, patch_in, pft, temp_cohort%n, temp_cohort%hite, & - temp_cohort%coage, temp_cohort%dbh, prt_obj, temp_cohort%leafmemory, & - temp_cohort%sapwmemory, temp_cohort%structmemory, cstatus, rstatus, & + temp_cohort%coage, temp_cohort%dbh, prt_obj, cstatus, rstatus, & temp_cohort%canopy_trim, temp_cohort%c_area,1, site_in%spread, bc_in) diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index 7d6d4f641b..e381ab8bc7 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -1690,6 +1690,15 @@ subroutine FatesCheckParams(is_master) end if end if + if( (prt_params%woody(ipft) == itrue) .and. & + (EDPftvarcon_inst%phen_stem_drop_fraction(ipft) > nearzero ) ) then + write(fates_log(),*) ' Non-zero stem-drop fractions are not allowed for woody plants' + write(fates_log(),*) ' PFT#: ',ipft + write(fates_log(),*) ' part_params%woody:',prt_params%woody(ipft) + write(fates_log(),*) ' phen_stem_drop_fraction: ', EDPFtvarcon_inst%phen_stem_drop_fraction(ipft) + write(fates_log(),*) ' Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if ! Check if freezing tolerance is within reasonable bounds ! ---------------------------------------------------------------------------------- diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 39b7c9cbf2..edf9976218 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -214,9 +214,6 @@ module EDTypesMod real(r8) :: coage ! cohort age in years real(r8) :: hite ! height: meters integer :: indexnumber ! unique number for each cohort. (within clump?) - real(r8) :: leafmemory ! target leaf biomass- set from previous year: kGC per indiv - real(r8) :: sapwmemory ! target sapwood biomass- set from previous year: kGC per indiv - real(r8) :: structmemory ! target structural biomass- set from previous year: kGC per indiv integer :: canopy_layer ! canopy status of cohort (1 = canopy, 2 = understorey, etc.) real(r8) :: canopy_layer_yesterday ! recent canopy status of cohort ! (1 = canopy, 2 = understorey, etc.) @@ -1037,9 +1034,6 @@ subroutine dump_cohort(ccohort) write(fates_log(),*) 'co%dbh = ', ccohort%dbh write(fates_log(),*) 'co%hite = ', ccohort%hite write(fates_log(),*) 'co%coage = ', ccohort%coage - write(fates_log(),*) 'co%leafmemory = ', ccohort%leafmemory - write(fates_log(),*) 'co%sapwmemory = ', ccohort%sapwmemory - write(fates_log(),*) 'co%structmemory = ', ccohort%structmemory write(fates_log(),*) 'leaf carbon = ', ccohort%prt%GetState(leaf_organ,all_carbon_elements) write(fates_log(),*) 'fineroot carbon = ', ccohort%prt%GetState(fnrt_organ,all_carbon_elements) diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index e4b5fe4e2b..a76ae9df4f 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -1041,32 +1041,22 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & call bstore_allom(temp_cohort%dbh, temp_cohort%pft, temp_cohort%canopy_trim, c_store) - temp_cohort%leafmemory = 0._r8 - temp_cohort%sapwmemory = 0._r8 - temp_cohort%structmemory = 0._r8 cstatus = leaves_on - stem_drop_fraction = EDPftvarcon_inst%phen_stem_drop_fraction(temp_cohort%pft) if( prt_params%season_decid(temp_cohort%pft) == itrue .and. & any(csite%cstatus == [phen_cstat_nevercold,phen_cstat_iscold])) then - temp_cohort%leafmemory = c_leaf - temp_cohort%sapwmemory = c_sapw * stem_drop_fraction - temp_cohort%structmemory = c_struct * stem_drop_fraction c_leaf = 0._r8 - c_sapw = (1._r8 - stem_drop_fraction) * c_sapw - c_struct = (1._r8 - stem_drop_fraction) * c_struct + c_sapw = (1._r8 - stem_drop_fraction) * c_sapw + c_struct = (1._r8 - stem_drop_fraction) * c_struct cstatus = leaves_off endif if ( prt_params%stress_decid(temp_cohort%pft) == itrue .and. & any(csite%dstatus == [phen_dstat_timeoff,phen_dstat_moistoff])) then - temp_cohort%leafmemory = c_leaf - temp_cohort%sapwmemory = c_sapw * stem_drop_fraction - temp_cohort%structmemory = c_struct * stem_drop_fraction c_leaf = 0._r8 - c_sapw = (1._r8 - stem_drop_fraction) * c_sapw - c_struct = (1._r8 - stem_drop_fraction) * c_struct + c_sapw = (1._r8 - stem_drop_fraction) * c_sapw + c_struct = (1._r8 - stem_drop_fraction) * c_struct cstatus = leaves_off endif @@ -1161,8 +1151,7 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & call create_cohort(csite, cpatch, temp_cohort%pft, temp_cohort%n, temp_cohort%hite, & temp_cohort%coage, temp_cohort%dbh, & - prt_obj, temp_cohort%leafmemory,temp_cohort%sapwmemory, temp_cohort%structmemory, & - cstatus, rstatus, temp_cohort%canopy_trim,temp_cohort%c_area, & + prt_obj, cstatus, rstatus, temp_cohort%canopy_trim,temp_cohort%c_area, & 1, csite%spread, bc_in) deallocate(temp_cohort) ! get rid of temporary cohort diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index c61c87fd4b..838441be8b 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -100,9 +100,6 @@ module FatesRestartInterfaceMod integer :: ir_coage_co integer :: ir_g_sb_laweight_co integer :: ir_height_co - integer :: ir_leafmemory_co - integer :: ir_sapwmemory_co - integer :: ir_structmemory_co integer :: ir_nplant_co integer :: ir_gpp_acc_co integer :: ir_npp_acc_co @@ -714,21 +711,6 @@ subroutine define_restart_vars(this, initialize_variables) long_name='ed cohort - plant height', units='m', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_height_co ) - call this%set_restart_var(vname='fates_leafmemory', vtype=cohort_r8, & - long_name='ed cohort - target leaf biomass set from prev year', & - units='kgC/indiv', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_leafmemory_co ) - - call this%set_restart_var(vname='fates_sapwmemory', vtype=cohort_r8, & - long_name='ed cohort - target sapwood biomass set from prev year', & - units='kgC/indiv', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_sapwmemory_co ) - - call this%set_restart_var(vname='fates_structmemory', vtype=cohort_r8, & - long_name='ed cohort - target structural biomass set from prev year', & - units='kgC/indiv', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_structmemory_co ) - call this%set_restart_var(vname='fates_nplant', vtype=cohort_r8, & long_name='ed cohort - number of plants in the cohort', & units='/patch', flushval = flushzero, & @@ -1773,9 +1755,6 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_coage_co => this%rvars(ir_coage_co)%r81d, & rio_g_sb_laweight_co => this%rvars(ir_g_sb_laweight_co)%r81d, & rio_height_co => this%rvars(ir_height_co)%r81d, & - rio_leafmemory_co => this%rvars(ir_leafmemory_co)%r81d, & - rio_sapwmemory_co => this%rvars(ir_sapwmemory_co)%r81d, & - rio_structmemory_co => this%rvars(ir_structmemory_co)%r81d, & rio_nplant_co => this%rvars(ir_nplant_co)%r81d, & rio_gpp_acc_co => this%rvars(ir_gpp_acc_co)%r81d, & rio_npp_acc_co => this%rvars(ir_npp_acc_co)%r81d, & @@ -2002,11 +1981,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_dbh_co(io_idx_co) = ccohort%dbh rio_coage_co(io_idx_co) = ccohort%coage rio_height_co(io_idx_co) = ccohort%hite - rio_leafmemory_co(io_idx_co) = ccohort%leafmemory - rio_sapwmemory_co(io_idx_co) = ccohort%sapwmemory - rio_structmemory_co(io_idx_co) = ccohort%structmemory rio_g_sb_laweight_co(io_idx_co)= ccohort%g_sb_laweight - rio_nplant_co(io_idx_co) = ccohort%n rio_gpp_acc_co(io_idx_co) = ccohort%gpp_acc rio_npp_acc_co(io_idx_co) = ccohort%npp_acc @@ -2612,9 +2587,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_coage_co => this%rvars(ir_coage_co)%r81d, & rio_g_sb_laweight_co => this%rvars(ir_g_sb_laweight_co)%r81d, & rio_height_co => this%rvars(ir_height_co)%r81d, & - rio_leafmemory_co => this%rvars(ir_leafmemory_co)%r81d, & - rio_sapwmemory_co => this%rvars(ir_sapwmemory_co)%r81d, & - rio_structmemory_co => this%rvars(ir_structmemory_co)%r81d, & rio_nplant_co => this%rvars(ir_nplant_co)%r81d, & rio_gpp_acc_co => this%rvars(ir_gpp_acc_co)%r81d, & rio_npp_acc_co => this%rvars(ir_npp_acc_co)%r81d, & @@ -2814,9 +2786,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ccohort%coage = rio_coage_co(io_idx_co) ccohort%g_sb_laweight= rio_g_sb_laweight_co(io_idx_co) ccohort%hite = rio_height_co(io_idx_co) - ccohort%leafmemory = rio_leafmemory_co(io_idx_co) - ccohort%sapwmemory = rio_sapwmemory_co(io_idx_co) - ccohort%structmemory = rio_structmemory_co(io_idx_co) ccohort%n = rio_nplant_co(io_idx_co) ccohort%gpp_acc = rio_gpp_acc_co(io_idx_co) ccohort%npp_acc = rio_npp_acc_co(io_idx_co) From d364757b99c31bbcb7d194696e9b6a50ad0d6148 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 28 Jul 2022 22:18:15 -0400 Subject: [PATCH 02/10] removed %lai and %sai where were redundant with treelai and canopy_area --- biogeochem/EDCanopyStructureMod.F90 | 436 +++++++++++----------------- biogeochem/EDCohortDynamicsMod.F90 | 6 - main/EDTypesMod.F90 | 3 - 3 files changed, 172 insertions(+), 273 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index f195f59f0a..6fb7529bf1 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1449,8 +1449,6 @@ subroutine leaf_area_profile( currentSite ) ! ! currentCohort%treelai ! LAI per unit crown area (m2/m2) ! currentCohort%treesai ! SAI per unit crown area (m2/m2) - ! currentCohort%lai ! LAI per unit canopy area (m2/m2) - ! currentCohort%sai ! SAI per unit canopy area (m2/m2) ! currentCohort%NV ! The number of discrete vegetation ! ! layers needed to describe this crown ! @@ -1477,7 +1475,7 @@ subroutine leaf_area_profile( currentSite ) ! !USES: use EDtypesMod , only : area, dinc_vai, dlower_vai, hitemax, n_hite_bins - + ! ! !ARGUMENTS type(ed_site_type) , intent(inout) :: currentSite @@ -1497,7 +1495,6 @@ subroutine leaf_area_profile( currentSite ) real(r8) :: layer_bottom_hite ! notional bottom height of this canopy layer (m) integer :: smooth_leaf_distribution ! is the leaf distribution this option (1) or not (0) real(r8) :: frac_canopy(N_HITE_BINS) ! amount of canopy in each height class - real(r8) :: patch_lai ! LAI summed over the patch in m2/m2 of canopy area real(r8) :: minh(N_HITE_BINS) ! minimum height in height class (m) real(r8) :: maxh(N_HITE_BINS) ! maximum height in height class (m) real(r8) :: dh ! vertical detph of height class (m) @@ -1525,7 +1522,6 @@ subroutine leaf_area_profile( currentSite ) currentPatch%canopy_layer_tlai(:) = 0._r8 currentPatch%ncan(:,:) = 0 currentPatch%nrad(:,:) = 0 - patch_lai = 0._r8 currentPatch%tlai_profile(:,:,:) = 0._r8 currentPatch%tsai_profile(:,:,:) = 0._r8 currentPatch%elai_profile(:,:,:) = 0._r8 @@ -1541,307 +1537,223 @@ subroutine leaf_area_profile( currentSite ) if (currentPatch%total_canopy_area > nearzero ) then - call UpdatePatchLAI(currentPatch, patch_lai) - - if(smooth_leaf_distribution == 1)then - - ! ----------------------------------------------------------------------------- - ! we are going to ignore the concept of canopy layers, and put all of the leaf - ! area into height banded bins. using the same domains as we had before, except - ! that CL always = 1 - ! ----------------------------------------------------------------------------- - - ! this is a crude way of dividing up the bins. Should it be a function of actual maximum height? - dh = 1.0_r8*(HITEMAX/N_HITE_BINS) - do iv = 1,N_HITE_BINS - if (iv == 1) then - minh(iv) = 0.0_r8 - maxh(iv) = dh - else - minh(iv) = (iv-1)*dh - maxh(iv) = (iv)*dh - endif - enddo - - currentCohort => currentPatch%shortest - do while(associated(currentCohort)) - ft = currentCohort%pft - min_chite = currentCohort%hite - currentCohort%hite * prt_params%crown_depth_frac(ft) - max_chite = currentCohort%hite - do iv = 1,N_HITE_BINS - frac_canopy(iv) = 0.0_r8 - ! this layer is in the middle of the canopy - if(max_chite > maxh(iv).and.min_chite < minh(iv))then - frac_canopy(iv)= min(1.0_r8,dh / (currentCohort%hite*prt_params%crown_depth_frac(ft) )) - ! this is the layer with the bottom of the canopy in it. - elseif(min_chite < maxh(iv).and.min_chite > minh(iv).and.max_chite > maxh(iv))then - frac_canopy(iv) = (maxh(iv) -min_chite ) / (currentCohort%hite*prt_params%crown_depth_frac(ft) ) - ! this is the layer with the top of the canopy in it. - elseif(max_chite > minh(iv).and.max_chite < maxh(iv).and.min_chite < minh(iv))then - frac_canopy(iv) = (max_chite - minh(iv)) / (currentCohort%hite*prt_params%crown_depth_frac(ft)) - elseif(max_chite < maxh(iv).and.min_chite > minh(iv))then !the whole cohort is within this layer. - frac_canopy(iv) = 1.0_r8 - endif - - ! no m2 of leaf per m2 of ground in each height class - currentPatch%tlai_profile(1,ft,iv) = currentPatch%tlai_profile(1,ft,iv) + frac_canopy(iv) * & - currentCohort%lai - currentPatch%tsai_profile(1,ft,iv) = currentPatch%tsai_profile(1,ft,iv) + frac_canopy(iv) * & - currentCohort%sai + call UpdatePatchLAI(currentPatch) - !snow burial - if(currentSite%snow_depth > maxh(iv))then - fraction_exposed = 0._r8 - endif - if(currentSite%snow_depth < minh(iv))then - fraction_exposed = 1._r8 - endif - if(currentSite%snow_depth >= minh(iv) .and. currentSite%snow_depth <= maxh(iv)) then !only partly hidden... - fraction_exposed = 1._r8 - max(0._r8,(min(1.0_r8,(currentSite%snow_depth-minh(iv))/dh))) - endif - - currentPatch%elai_profile(1,ft,iv) = currentPatch%tlai_profile(1,ft,iv) * fraction_exposed - currentPatch%esai_profile(1,ft,iv) = currentPatch%tsai_profile(1,ft,iv) * fraction_exposed + ! ----------------------------------------------------------------------------- + ! Standard canopy layering model. + ! Go through all cohorts and add their leaf area + ! and canopy area to the accumulators. + ! ----------------------------------------------------------------------------- - enddo ! (iv) hite bins + currentCohort => currentPatch%shortest + do while(associated(currentCohort)) + ft = currentCohort%pft + cl = currentCohort%canopy_layer - currentCohort => currentCohort%taller + ! ---------------------------------------------------------------- + ! How much of each tree is stem area index? Assuming that there is + ! This may indeed be zero if there is a sensecent grass + ! ---------------------------------------------------------------- - enddo !currentCohort + if( (currentCohort%treelai+currentCohort%treesai) > 0._r8)then + fleaf = currentCohort%treelai / (currentCohort%treelai + currentCohort%treesai) + else + fleaf = 0._r8 + endif - ! ----------------------------------------------------------------------------- - ! Perform a leaf area conservation check on the LAI profile - lai = 0.0_r8 - do ft = 1,numpft - lai = lai+ sum(currentPatch%tlai_profile(1,ft,:)) - enddo + currentPatch%nrad(cl,ft) = currentPatch%ncan(cl,ft) - if(lai > patch_lai)then - write(fates_log(), *) 'FATES: problem with lai assignments' + if (currentPatch%nrad(cl,ft) > nlevleaf ) then + write(fates_log(), *) 'Number of radiative leaf layers is larger' + write(fates_log(), *) ' than the maximum allowed.' + write(fates_log(), *) ' cl: ',cl + write(fates_log(), *) ' ft: ',ft + write(fates_log(), *) ' nlevleaf: ',nlevleaf + write(fates_log(), *) ' currentPatch%nrad(cl,ft): ', currentPatch%nrad(cl,ft) call endrun(msg=errMsg(sourcefile, __LINE__)) - endif + end if - else ! smooth leaf distribution + ! -------------------------------------------------------------------------- + ! Whole layers. Make a weighted average of the leaf area in each layer + ! before dividing it by the total area. Fill up layer for whole layers. + ! -------------------------------------------------------------------------- - ! ----------------------------------------------------------------------------- - ! Standard canopy layering model. - ! Go through all cohorts and add their leaf area - ! and canopy area to the accumulators. - ! ----------------------------------------------------------------------------- + do iv = 1,currentCohort%NV + ! This loop builds the arrays that define the effective (not snow covered) + ! and total (includes snow covered) area indices for leaves and stems + ! We calculate the absolute elevation of each layer to help determine if the layer + ! is obscured by snow. - currentCohort => currentPatch%shortest - do while(associated(currentCohort)) - ft = currentCohort%pft - cl = currentCohort%canopy_layer + layer_top_hite = currentCohort%hite - & + ( real(iv-1,r8)/currentCohort%NV * currentCohort%hite * & + prt_params%crown_depth_frac(currentCohort%pft) ) - ! ---------------------------------------------------------------- - ! How much of each tree is stem area index? Assuming that there is - ! This may indeed be zero if there is a sensecent grass - ! ---------------------------------------------------------------- + layer_bottom_hite = currentCohort%hite - & + ( real(iv,r8)/currentCohort%NV * currentCohort%hite * & + prt_params%crown_depth_frac(currentCohort%pft) ) - if( (currentCohort%treelai+currentCohort%treesai) > 0._r8)then - fleaf = currentCohort%lai / (currentCohort%lai + currentCohort%sai) - else - fleaf = 0._r8 + fraction_exposed = 1.0_r8 + if(currentSite%snow_depth > layer_top_hite)then + fraction_exposed = 0._r8 + endif + if(currentSite%snow_depth < layer_bottom_hite)then + fraction_exposed = 1._r8 + endif + if(currentSite%snow_depth >= layer_bottom_hite .and. & + currentSite%snow_depth <= layer_top_hite) then !only partly hidden... + fraction_exposed = 1._r8 - max(0._r8,(min(1.0_r8,(currentSite%snow_depth -layer_bottom_hite)/ & + (layer_top_hite-layer_bottom_hite )))) endif - currentPatch%nrad(cl,ft) = currentPatch%ncan(cl,ft) + if(iv==currentCohort%NV) then + remainder = (currentCohort%treelai + currentCohort%treesai) - & + (dlower_vai(iv) - dinc_vai(iv)) + if(remainder > dinc_vai(iv) )then + write(fates_log(), *)'ED: issue with remainder', & + currentCohort%treelai,currentCohort%treesai,dinc_vai(iv), & + currentCohort%NV,remainder - if (currentPatch%nrad(cl,ft) > nlevleaf ) then - write(fates_log(), *) 'Number of radiative leaf layers is larger' - write(fates_log(), *) ' than the maximum allowed.' - write(fates_log(), *) ' cl: ',cl - write(fates_log(), *) ' ft: ',ft - write(fates_log(), *) ' nlevleaf: ',nlevleaf - write(fates_log(), *) ' currentPatch%nrad(cl,ft): ', currentPatch%nrad(cl,ft) - call endrun(msg=errMsg(sourcefile, __LINE__)) + call endrun(msg=errMsg(sourcefile, __LINE__)) + endif + else + remainder = dinc_vai(iv) end if + currentPatch%tlai_profile(cl,ft,iv) = currentPatch%tlai_profile(cl,ft,iv) + & + remainder * fleaf * currentCohort%c_area/currentPatch%total_canopy_area - ! -------------------------------------------------------------------------- - ! Whole layers. Make a weighted average of the leaf area in each layer - ! before dividing it by the total area. Fill up layer for whole layers. - ! -------------------------------------------------------------------------- + currentPatch%elai_profile(cl,ft,iv) = currentPatch%elai_profile(cl,ft,iv) + & + remainder * fleaf * currentCohort%c_area/currentPatch%total_canopy_area * & + fraction_exposed - do iv = 1,currentCohort%NV + currentPatch%tsai_profile(cl,ft,iv) = currentPatch%tsai_profile(cl,ft,iv) + & + remainder * (1._r8 - fleaf) * currentCohort%c_area/currentPatch%total_canopy_area - ! This loop builds the arrays that define the effective (not snow covered) - ! and total (includes snow covered) area indices for leaves and stems - ! We calculate the absolute elevation of each layer to help determine if the layer - ! is obscured by snow. + currentPatch%esai_profile(cl,ft,iv) = currentPatch%esai_profile(cl,ft,iv) + & + remainder * (1._r8 - fleaf) * currentCohort%c_area/currentPatch%total_canopy_area * & + fraction_exposed - layer_top_hite = currentCohort%hite - & - ( real(iv-1,r8)/currentCohort%NV * currentCohort%hite * & - prt_params%crown_depth_frac(currentCohort%pft) ) + currentPatch%canopy_area_profile(cl,ft,iv) = currentPatch%canopy_area_profile(cl,ft,iv) + & + currentCohort%c_area/currentPatch%total_canopy_area - layer_bottom_hite = currentCohort%hite - & - ( real(iv,r8)/currentCohort%NV * currentCohort%hite * & - prt_params%crown_depth_frac(currentCohort%pft) ) + currentPatch%layer_height_profile(cl,ft,iv) = currentPatch%layer_height_profile(cl,ft,iv) + & + (remainder * fleaf * currentCohort%c_area/currentPatch%total_canopy_area * & + (layer_top_hite+layer_bottom_hite)/2.0_r8) !average height of layer. - fraction_exposed = 1.0_r8 - if(currentSite%snow_depth > layer_top_hite)then - fraction_exposed = 0._r8 - endif - if(currentSite%snow_depth < layer_bottom_hite)then - fraction_exposed = 1._r8 - endif - if(currentSite%snow_depth >= layer_bottom_hite .and. & - currentSite%snow_depth <= layer_top_hite) then !only partly hidden... - fraction_exposed = 1._r8 - max(0._r8,(min(1.0_r8,(currentSite%snow_depth -layer_bottom_hite)/ & - (layer_top_hite-layer_bottom_hite )))) - endif - - if(iv==currentCohort%NV) then - remainder = (currentCohort%treelai + currentCohort%treesai) - & - (dlower_vai(iv) - dinc_vai(iv)) - if(remainder > dinc_vai(iv) )then - write(fates_log(), *)'ED: issue with remainder', & - currentCohort%treelai,currentCohort%treesai,dinc_vai(iv), & - currentCohort%NV,remainder - - call endrun(msg=errMsg(sourcefile, __LINE__)) - endif - else - remainder = dinc_vai(iv) - end if + end do - currentPatch%tlai_profile(cl,ft,iv) = currentPatch%tlai_profile(cl,ft,iv) + & - remainder * fleaf * currentCohort%c_area/currentPatch%total_canopy_area + currentCohort => currentCohort%taller - currentPatch%elai_profile(cl,ft,iv) = currentPatch%elai_profile(cl,ft,iv) + & - remainder * fleaf * currentCohort%c_area/currentPatch%total_canopy_area * & - fraction_exposed + enddo !cohort - currentPatch%tsai_profile(cl,ft,iv) = currentPatch%tsai_profile(cl,ft,iv) + & - remainder * (1._r8 - fleaf) * currentCohort%c_area/currentPatch%total_canopy_area + ! -------------------------------------------------------------------------- - currentPatch%esai_profile(cl,ft,iv) = currentPatch%esai_profile(cl,ft,iv) + & - remainder * (1._r8 - fleaf) * currentCohort%c_area/currentPatch%total_canopy_area * & - fraction_exposed + ! If there is an upper-story, the top canopy layer + ! should have a value of exactly 1.0 in its top leaf layer + ! -------------------------------------------------------------------------- - currentPatch%canopy_area_profile(cl,ft,iv) = currentPatch%canopy_area_profile(cl,ft,iv) + & + if ( (currentPatch%NCL_p > 1) .and. & + (sum(currentPatch%canopy_area_profile(1,:,1)) < 0.9999 )) then + write(fates_log(), *) 'FATES: canopy_area_profile was less than 1 at the canopy top' + write(fates_log(), *) 'cl: ',1 + write(fates_log(), *) 'iv: ',1 + write(fates_log(), *) 'sum(cpatch%canopy_area_profile(1,:,1)): ', & + sum(currentPatch%canopy_area_profile(1,:,1)) + currentCohort => currentPatch%shortest + do while(associated(currentCohort)) + if(currentCohort%canopy_layer==1)then + write(fates_log(), *) 'FATES: cohorts',currentCohort%dbh,currentCohort%c_area, & + currentPatch%total_canopy_area,currentPatch%area + write(fates_log(), *) 'ED: fracarea', currentCohort%pft, & currentCohort%c_area/currentPatch%total_canopy_area - - currentPatch%layer_height_profile(cl,ft,iv) = currentPatch%layer_height_profile(cl,ft,iv) + & - (remainder * fleaf * currentCohort%c_area/currentPatch%total_canopy_area * & - (layer_top_hite+layer_bottom_hite)/2.0_r8) !average height of layer. - - end do - + endif currentCohort => currentCohort%taller + enddo !currentCohort + call endrun(msg=errMsg(sourcefile, __LINE__)) - enddo !cohort - - ! -------------------------------------------------------------------------- + end if - ! If there is an upper-story, the top canopy layer - ! should have a value of exactly 1.0 in its top leaf layer - ! -------------------------------------------------------------------------- - if ( (currentPatch%NCL_p > 1) .and. & - (sum(currentPatch%canopy_area_profile(1,:,1)) < 0.9999 )) then - write(fates_log(), *) 'FATES: canopy_area_profile was less than 1 at the canopy top' - write(fates_log(), *) 'cl: ',1 - write(fates_log(), *) 'iv: ',1 - write(fates_log(), *) 'sum(cpatch%canopy_area_profile(1,:,1)): ', & - sum(currentPatch%canopy_area_profile(1,:,1)) - currentCohort => currentPatch%shortest - do while(associated(currentCohort)) - if(currentCohort%canopy_layer==1)then - write(fates_log(), *) 'FATES: cohorts',currentCohort%dbh,currentCohort%c_area, & - currentPatch%total_canopy_area,currentPatch%area - write(fates_log(), *) 'ED: fracarea', currentCohort%pft, & - currentCohort%c_area/currentPatch%total_canopy_area - endif - currentCohort => currentCohort%taller - enddo !currentCohort - call endrun(msg=errMsg(sourcefile, __LINE__)) + ! -------------------------------------------------------------------------- + ! In the following loop we are now normalizing the effective and + ! total area profiles to convert from units of leaf/stem area per vegetated + ! canopy area, into leaf/stem area per area of their own radiative column + ! which is typically the footprint of all cohorts contained in the canopy + ! layer x pft bins. + ! Also perform some checks on area normalization. + ! Check the area of each leaf layer, across pfts. + ! It should never be larger than 1 or less than 0. + ! -------------------------------------------------------------------------- - end if + do cl = 1,currentPatch%NCL_p + do iv = 1,currentPatch%ncan(cl,ft) + if( debug .and. sum(currentPatch%canopy_area_profile(cl,:,iv)) > 1.0001_r8 ) then - ! -------------------------------------------------------------------------- - ! In the following loop we are now normalizing the effective and - ! total area profiles to convert from units of leaf/stem area per vegetated - ! canopy area, into leaf/stem area per area of their own radiative column - ! which is typically the footprint of all cohorts contained in the canopy - ! layer x pft bins. - ! Also perform some checks on area normalization. - ! Check the area of each leaf layer, across pfts. - ! It should never be larger than 1 or less than 0. - ! -------------------------------------------------------------------------- + write(fates_log(), *) 'FATES: A canopy_area_profile exceeded 1.0' + write(fates_log(), *) 'cl: ',cl + write(fates_log(), *) 'iv: ',iv + write(fates_log(), *) 'sum(cpatch%canopy_area_profile(cl,:,iv)): ', & + sum(currentPatch%canopy_area_profile(cl,:,iv)) + currentCohort => currentPatch%shortest + do while(associated(currentCohort)) + if(currentCohort%canopy_layer==cl)then + write(fates_log(), *) 'FATES: cohorts in layer cl = ',cl, & + currentCohort%dbh,currentCohort%c_area, & + currentPatch%total_canopy_area,currentPatch%area + write(fates_log(), *) 'ED: fracarea', currentCohort%pft, & + currentCohort%c_area/currentPatch%total_canopy_area + endif + currentCohort => currentCohort%taller + enddo !currentCohort + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end do - do cl = 1,currentPatch%NCL_p + do ft = 1,numpft do iv = 1,currentPatch%ncan(cl,ft) - if( debug .and. sum(currentPatch%canopy_area_profile(cl,:,iv)) > 1.0001_r8 ) then - - write(fates_log(), *) 'FATES: A canopy_area_profile exceeded 1.0' - write(fates_log(), *) 'cl: ',cl - write(fates_log(), *) 'iv: ',iv - write(fates_log(), *) 'sum(cpatch%canopy_area_profile(cl,:,iv)): ', & - sum(currentPatch%canopy_area_profile(cl,:,iv)) - currentCohort => currentPatch%shortest - do while(associated(currentCohort)) - if(currentCohort%canopy_layer==cl)then - write(fates_log(), *) 'FATES: cohorts in layer cl = ',cl, & - currentCohort%dbh,currentCohort%c_area, & - currentPatch%total_canopy_area,currentPatch%area - write(fates_log(), *) 'ED: fracarea', currentCohort%pft, & - currentCohort%c_area/currentPatch%total_canopy_area - endif - currentCohort => currentCohort%taller - enddo !currentCohort - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - end do + if( currentPatch%canopy_area_profile(cl,ft,iv) > nearzero ) then - do ft = 1,numpft - do iv = 1,currentPatch%ncan(cl,ft) + currentPatch%tlai_profile(cl,ft,iv) = currentPatch%tlai_profile(cl,ft,iv) / & + currentPatch%canopy_area_profile(cl,ft,iv) - if( currentPatch%canopy_area_profile(cl,ft,iv) > nearzero ) then + currentPatch%tsai_profile(cl,ft,iv) = currentPatch%tsai_profile(cl,ft,iv) / & + currentPatch%canopy_area_profile(cl,ft,iv) - currentPatch%tlai_profile(cl,ft,iv) = currentPatch%tlai_profile(cl,ft,iv) / & - currentPatch%canopy_area_profile(cl,ft,iv) + currentPatch%elai_profile(cl,ft,iv) = currentPatch%elai_profile(cl,ft,iv) / & + currentPatch%canopy_area_profile(cl,ft,iv) - currentPatch%tsai_profile(cl,ft,iv) = currentPatch%tsai_profile(cl,ft,iv) / & - currentPatch%canopy_area_profile(cl,ft,iv) - - currentPatch%elai_profile(cl,ft,iv) = currentPatch%elai_profile(cl,ft,iv) / & - currentPatch%canopy_area_profile(cl,ft,iv) - - currentPatch%esai_profile(cl,ft,iv) = currentPatch%esai_profile(cl,ft,iv) / & - currentPatch%canopy_area_profile(cl,ft,iv) - end if - - if(currentPatch%tlai_profile(cl,ft,iv)>nearzero )then - currentPatch%layer_height_profile(cl,ft,iv) = currentPatch%layer_height_profile(cl,ft,iv) & - /currentPatch%tlai_profile(cl,ft,iv) - end if + currentPatch%esai_profile(cl,ft,iv) = currentPatch%esai_profile(cl,ft,iv) / & + currentPatch%canopy_area_profile(cl,ft,iv) + end if - enddo + if(currentPatch%tlai_profile(cl,ft,iv)>nearzero )then + currentPatch%layer_height_profile(cl,ft,iv) = currentPatch%layer_height_profile(cl,ft,iv) & + /currentPatch%tlai_profile(cl,ft,iv) + end if enddo - enddo - ! -------------------------------------------------------------------------- - ! Set the mask that identifies which PFT x can-layer combinations have - ! scattering elements in them. - ! -------------------------------------------------------------------------- + enddo + enddo - do cl = 1,currentPatch%NCL_p - do ft = 1,numpft - do iv = 1, currentPatch%nrad(cl,ft) - if(currentPatch%canopy_area_profile(cl,ft,iv) > 0._r8)then - currentPatch%canopy_mask(cl,ft) = 1 - endif - end do !iv - enddo !ft - enddo ! loop over cl + ! -------------------------------------------------------------------------- + ! Set the mask that identifies which PFT x can-layer combinations have + ! scattering elements in them. + ! -------------------------------------------------------------------------- - endif !leaf distribution + do cl = 1,currentPatch%NCL_p + do ft = 1,numpft + do iv = 1, currentPatch%nrad(cl,ft) + if(currentPatch%canopy_area_profile(cl,ft,iv) > 0._r8)then + currentPatch%canopy_mask(cl,ft) = 1 + endif + end do !iv + enddo !ft + enddo ! loop over cl end if @@ -2180,7 +2092,7 @@ end subroutine CanopyLayerArea ! =============================================================================================== - subroutine UpdatePatchLAI(currentPatch, patch_lai) + subroutine UpdatePatchLAI(currentPatch) ! -------------------------------------------------------------------------------------------- ! This subroutine works through the current patch cohorts and updates the canopy_layer_tlai @@ -2192,7 +2104,6 @@ subroutine UpdatePatchLAI(currentPatch, patch_lai) ! Arguments type(ed_patch_type),intent(inout), target :: currentPatch - real(r8), intent(inout) :: patch_lai ! Local Variables type(ed_cohort_type), pointer :: currentCohort @@ -2219,11 +2130,10 @@ subroutine UpdatePatchLAI(currentPatch, patch_lai) ! Update the number of number of vegetation layers currentPatch%ncan(cl,ft) = max(currentPatch%ncan(cl,ft),currentCohort%NV) - ! Update the patch canopy layer tlai - currentPatch%canopy_layer_tlai(cl) = currentPatch%canopy_layer_tlai(cl) + currentCohort%lai + ! Update the patch canopy layer tlai (LAI per canopy area) + currentPatch%canopy_layer_tlai(cl) = currentPatch%canopy_layer_tlai(cl) + & + currentCohort%treelai *currentCohort%c_area/currentPatch%total_canopy_area - ! Calculate the total patch lai - patch_lai = patch_lai + currentCohort%lai end if currentCohort => currentCohort%shorter @@ -2263,9 +2173,7 @@ subroutine UpdateCohortLAI(currentCohort, canopy_layer_tlai, patcharea) currentCohort%vcmax25top,4) end if - ! Update the cohort lai and sai - currentCohort%lai = currentCohort%treelai *currentCohort%c_area/patcharea - currentCohort%sai = currentCohort%treesai *currentCohort%c_area/patcharea + ! Number of actual vegetation layers in this cohort's crown currentCohort%nv = count((currentCohort%treelai+currentCohort%treesai) .gt. dlower_vai(:)) + 1 diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index a7d203da5d..2d796a9048 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -277,8 +277,6 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, coage, dbh, & patchptr%canopy_layer_tlai, new_cohort%treelai,new_cohort%vcmax25top,2 ) end if - new_cohort%lai = new_cohort%treelai * new_cohort%c_area/patchptr%area - ! Put cohort at the right place in the linked list storebigcohort => patchptr%tallest @@ -534,8 +532,6 @@ subroutine nan_cohort(cc_p) currentCohort%dbh = nan ! 'diameter at breast height' in cm currentCohort%coage = nan ! age of the cohort in years currentCohort%hite = nan ! height: meters - currentCohort%lai = nan ! leaf area index of cohort m2/m2 - currentCohort%sai = nan ! stem area index of cohort m2/m2 currentCohort%g_sb_laweight = nan ! Total leaf conductance of cohort (stomata+blayer) weighted by leaf-area [m/s]*[m2] currentCohort%canopy_trim = nan ! What is the fraction of the maximum leaf biomass that we are targeting? :- currentCohort%leaf_cost = nan ! How much does it cost to maintain leaves: kgC/m2/year-1 @@ -1787,8 +1783,6 @@ subroutine copy_cohort( currentCohort,copyc ) n%dbh = o%dbh n%coage = o%coage n%hite = o%hite - n%lai = o%lai - n%sai = o%sai n%g_sb_laweight = o%g_sb_laweight n%leaf_cost = o%leaf_cost n%canopy_layer = o%canopy_layer diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index edf9976218..dfdfcdf43b 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -1041,9 +1041,6 @@ subroutine dump_cohort(ccohort) write(fates_log(),*) 'structural (dead) carbon = ', ccohort%prt%GetState(struct_organ,all_carbon_elements) write(fates_log(),*) 'storage carbon = ', ccohort%prt%GetState(store_organ,all_carbon_elements) write(fates_log(),*) 'reproductive carbon = ', ccohort%prt%GetState(repro_organ,all_carbon_elements) - - write(fates_log(),*) 'co%lai = ', ccohort%lai - write(fates_log(),*) 'co%sai = ', ccohort%sai write(fates_log(),*) 'co%g_sb_laweight = ', ccohort%g_sb_laweight write(fates_log(),*) 'co%leaf_cost = ', ccohort%leaf_cost write(fates_log(),*) 'co%canopy_layer = ', ccohort%canopy_layer From 3c5090d2f8343971f64f652e9305154af982e256 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 28 Jul 2022 22:35:49 -0400 Subject: [PATCH 03/10] Cleaned up logic on flushing for grasses without stemdrop --- biogeochem/EDPhysiologyMod.F90 | 100 +++++++++++++++------------------ 1 file changed, 45 insertions(+), 55 deletions(-) diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index d3815241fc..f3072cd4eb 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -1108,10 +1108,12 @@ subroutine phenology_leafonoff(currentSite) call bbgw_allom(currentCohort%dbh,currentCohort%pft,target_bgw_c) call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, currentCohort%pft, target_struct_c) - if(prt_params%woody(ipft) == itrue)then - deficit_c = target_leaf_c - else + if (stem_drop_fraction .gt. 0.0_r8) then + ! Note, this is only true for some grasses, woody plants don't + ! have a stem drop fraction deficit_c = target_leaf_c + (target_sapw_c-sapw_c) + (target_struct_c-struct_c) + else + deficit_c = target_leaf_c end if if(store_c>nearzero) then @@ -1128,33 +1130,27 @@ subroutine phenology_leafonoff(currentSite) ! This call will request that storage carbon will be transferred to ! leaf tissues. It is specified as a fraction of the available storage - if(prt_params%woody(ipft) == itrue) then - - call PRTPhenologyFlush(currentCohort%prt, ipft, leaf_organ, store_c_transfer_frac) - + ! Check that the stem drop fraction is set to non-zero amount + ! otherwise flush all carbon store to leaves + if (stem_drop_fraction .gt. 0.0_r8) then + + call PRTPhenologyFlush(currentCohort%prt, ipft, leaf_organ, & + store_c_transfer_frac*target_leaf_c/deficit_c) + + call PRTPhenologyFlush(currentCohort%prt, ipft, sapw_organ, & + store_c_transfer_frac*(target_sapw_c-sapw_c)/deficit_c) + + call PRTPhenologyFlush(currentCohort%prt, ipft, struct_organ, & + store_c_transfer_frac*(target_struct_c-struct_c)/deficit_c) + else + + call PRTPhenologyFlush(currentCohort%prt, ipft, leaf_organ, & + store_c_transfer_frac) + + end if - ! Check that the stem drop fraction is set to non-zero amount - ! otherwise flush all carbon store to leaves - if (stem_drop_fraction .gt. 0.0_r8) then - - call PRTPhenologyFlush(currentCohort%prt, ipft, leaf_organ, & - store_c_transfer_frac*target_leaf_c/deficit_c) - - call PRTPhenologyFlush(currentCohort%prt, ipft, sapw_organ, & - store_c_transfer_frac*(target_sapw_c-sapw_c)/deficit_c) - - call PRTPhenologyFlush(currentCohort%prt, ipft, struct_organ, & - store_c_transfer_frac*(target_struct_c-struct_c)/deficit_c) - - else - - call PRTPhenologyFlush(currentCohort%prt, ipft, leaf_organ, & - store_c_transfer_frac) - - end if - - endif + endif endif !pft phenology endif ! growing season @@ -1214,10 +1210,12 @@ subroutine phenology_leafonoff(currentSite) call bbgw_allom(currentCohort%dbh,currentCohort%pft,target_bgw_c) call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, currentCohort%pft, target_struct_c) - if(prt_params%woody(ipft) == itrue)then - deficit_c = target_leaf_c - else + if (stem_drop_fraction .gt. 0.0_r8) then + ! Note, this is only true for some grasses, woody plants don't + ! have a stem drop fraction deficit_c = target_leaf_c + (target_sapw_c-sapw_c) + (target_struct_c-struct_c) + else + deficit_c = target_leaf_c end if if(store_c>nearzero) then @@ -1232,31 +1230,23 @@ subroutine phenology_leafonoff(currentSite) ! This call will request that storage carbon will be transferred to ! leaf tissues. It is specified as a fraction of the available storage - if(prt_params%woody(ipft) == itrue) then - - call PRTPhenologyFlush(currentCohort%prt, ipft, & - leaf_organ, store_c_transfer_frac) - + if (stem_drop_fraction .gt. 0.0_r8) then + + call PRTPhenologyFlush(currentCohort%prt, ipft, leaf_organ, & + store_c_transfer_frac*target_leaf_c/deficit_c) + + call PRTPhenologyFlush(currentCohort%prt, ipft, sapw_organ, & + store_c_transfer_frac*(target_sapw_c-sapw_c)/deficit_c) + + call PRTPhenologyFlush(currentCohort%prt, ipft, struct_organ, & + store_c_transfer_frac*(target_struct_c-struct_c)/deficit_c) + else - - ! Check that the stem drop fraction is set to non-zero amount otherwise flush all carbon store to leaves - if (stem_drop_fraction .gt. 0.0_r8) then - - call PRTPhenologyFlush(currentCohort%prt, ipft, leaf_organ, & - store_c_transfer_frac*target_leaf_c/deficit_c) - - call PRTPhenologyFlush(currentCohort%prt, ipft, sapw_organ, & - store_c_transfer_frac*(target_sapw_c-sapw_c)/deficit_c) - - call PRTPhenologyFlush(currentCohort%prt, ipft, struct_organ, & - store_c_transfer_frac*(target_struct_c-struct_c)/deficit_c) - - else - - call PRTPhenologyFlush(currentCohort%prt, ipft, leaf_organ, & - store_c_transfer_frac) - - end if + + call PRTPhenologyFlush(currentCohort%prt, ipft, leaf_organ, & + store_c_transfer_frac) + + end if endif ! woody plant check endif !currentCohort status again? From a37749b814bcb6da824be224ca37fb600082dfd5 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 29 Jul 2022 10:20:09 -0400 Subject: [PATCH 04/10] Added call to bagw during leaf flushing --- biogeochem/EDPhysiologyMod.F90 | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index f3072cd4eb..a47d1fe433 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -1066,7 +1066,7 @@ subroutine phenology_leafonoff(currentSite) real(r8) :: deficit_c ! Amount of C needed to get flushing pools "on-allometry" real(r8) :: target_leaf_c real(r8) :: target_sapw_c - real(r8) :: target_bgw_c, target_agw_c, target_struct_c + real(r8) :: target_agw_c, target_bgw_c, target_agw_c, target_struct_c real(r8) :: sapw_area integer :: ipft real(r8), parameter :: leaf_drop_fraction = 1.0_r8 @@ -1104,9 +1104,12 @@ subroutine phenology_leafonoff(currentSite) ! stop flow of carbon out of bstore. call bleaf(currentCohort%dbh,currentCohort%pft,currentCohort%canopy_trim,target_leaf_c) - call bsap_allom(currentCohort%dbh,currentCohort%pft,currentCohort%canopy_trim,sapw_area,target_sapw_c) + call bsap_allom(currentCohort%dbh,currentCohort%pft, & + currentCohort%canopy_trim,sapw_area,target_sapw_c) + call bagw_allom(currentCohort%dbh,currentCohort%pft,target_agw_c) call bbgw_allom(currentCohort%dbh,currentCohort%pft,target_bgw_c) - call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, currentCohort%pft, target_struct_c) + call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, & + currentCohort%pft, target_struct_c) if (stem_drop_fraction .gt. 0.0_r8) then ! Note, this is only true for some grasses, woody plants don't @@ -1205,10 +1208,14 @@ subroutine phenology_leafonoff(currentSite) currentCohort%status_coh = leaves_on ! Leaves are on, so change status to ! stop flow of carbon out of bstore. - call bleaf(currentCohort%dbh,currentCohort%pft,currentCohort%canopy_trim,target_leaf_c) - call bsap_allom(currentCohort%dbh,currentCohort%pft,currentCohort%canopy_trim,sapw_area,target_sapw_c) + call bleaf(currentCohort%dbh,currentCohort%pft,& + currentCohort%canopy_trim,target_leaf_c) + call bsap_allom(currentCohort%dbh,currentCohort%pft, & + currentCohort%canopy_trim,sapw_area,target_sapw_c) + call bagw_allom(currentCohort%dbh,currentCohort%pft,target_agw_c) call bbgw_allom(currentCohort%dbh,currentCohort%pft,target_bgw_c) - call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, currentCohort%pft, target_struct_c) + call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, & + currentCohort%pft, target_struct_c) if (stem_drop_fraction .gt. 0.0_r8) then ! Note, this is only true for some grasses, woody plants don't From d4b7913e9eb7a5b571bd5a54fa9afa0a63cd6fa9 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Sun, 31 Jul 2022 10:08:46 -0600 Subject: [PATCH 05/10] Removed extra end-if closures I added while updating phenology logic --- biogeochem/EDPhysiologyMod.F90 | 28 +++++++++++++--------------- 1 file changed, 13 insertions(+), 15 deletions(-) diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index a47d1fe433..4ecd61c340 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -1066,7 +1066,7 @@ subroutine phenology_leafonoff(currentSite) real(r8) :: deficit_c ! Amount of C needed to get flushing pools "on-allometry" real(r8) :: target_leaf_c real(r8) :: target_sapw_c - real(r8) :: target_agw_c, target_bgw_c, target_agw_c, target_struct_c + real(r8) :: target_agw_c, target_bgw_c, target_struct_c real(r8) :: sapw_area integer :: ipft real(r8), parameter :: leaf_drop_fraction = 1.0_r8 @@ -1097,9 +1097,9 @@ subroutine phenology_leafonoff(currentSite) ! The site level flags signify that it is no-longer too cold ! for leaves. Time to signal flushing - if (prt_params%season_decid(ipft) == itrue)then - if ( currentSite%cstatus == phen_cstat_notcold )then ! we have just moved to leaves being on . - if (currentCohort%status_coh == leaves_off)then ! Are the leaves currently off? + if_colddec: if (prt_params%season_decid(ipft) == itrue)then + if_notcold: if ( currentSite%cstatus == phen_cstat_notcold )then ! we have just moved to leaves being on . + if_leaves_off: if (currentCohort%status_coh == leaves_off)then ! Are the leaves currently off? currentCohort%status_coh = leaves_on ! Leaves are on, so change status to ! stop flow of carbon out of bstore. @@ -1153,16 +1153,15 @@ subroutine phenology_leafonoff(currentSite) end if - endif - endif !pft phenology - endif ! growing season + endif if_leaves_off + endif if_notcold !COLD LEAF OFF - if (currentSite%cstatus == phen_cstat_nevercold .or. & + if_cold: if (currentSite%cstatus == phen_cstat_nevercold .or. & currentSite%cstatus == phen_cstat_iscold) then ! past leaf drop day? Leaves still on tree? - - if (currentCohort%status_coh == leaves_on) then ! leaves have not dropped - + + if_leaves_on: if (currentCohort%status_coh == leaves_on) then ! leaves have not dropped + ! leaf off occur on individuals bigger than specific size for grass if (currentCohort%dbh > EDPftvarcon_inst%phen_cold_size_threshold(ipft) & .or. prt_params%woody(ipft)==itrue) then @@ -1187,9 +1186,9 @@ subroutine phenology_leafonoff(currentSite) endif ! woody plant check endif ! individual dbh size check - endif !leaf status - endif !currentSite status - endif !season_decid + endif if_leaves_on !leaf status + endif if_cold !currentSite status + endif if_colddec !season_decid ! DROUGHT LEAF ON ! Site level flag indicates it is no longer in drought condition @@ -1255,7 +1254,6 @@ subroutine phenology_leafonoff(currentSite) end if - endif ! woody plant check endif !currentCohort status again? endif !currentSite status From 73596b7a35ea65d6595fcc46bd14eab085d0b8a8 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 4 Aug 2022 12:36:02 -0400 Subject: [PATCH 06/10] Added nearzero to if logic --- biogeochem/EDCanopyStructureMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 6fb7529bf1..4066aec6b6 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1555,7 +1555,7 @@ subroutine leaf_area_profile( currentSite ) ! This may indeed be zero if there is a sensecent grass ! ---------------------------------------------------------------- - if( (currentCohort%treelai+currentCohort%treesai) > 0._r8)then + if( (currentCohort%treelai+currentCohort%treesai) > nearzero)then fleaf = currentCohort%treelai / (currentCohort%treelai + currentCohort%treesai) else fleaf = 0._r8 From 82aa37007a9ceba5de3e7dab5be888a35fb2c6da Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 22 Aug 2022 14:06:14 -0400 Subject: [PATCH 07/10] Removed declaration of cohort structure lai and sai --- main/EDTypesMod.F90 | 3 --- 1 file changed, 3 deletions(-) diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index dfdfcdf43b..16ceda6115 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -218,9 +218,6 @@ module EDTypesMod real(r8) :: canopy_layer_yesterday ! recent canopy status of cohort ! (1 = canopy, 2 = understorey, etc.) ! real to be conservative during fusion - - real(r8) :: lai ! leaf area index of cohort: m2 leaf area of entire cohort per m2 of canopy area of a patch - real(r8) :: sai ! stem area index of cohort: m2 leaf area of entire cohort per m2 of canopy area of a patch real(r8) :: g_sb_laweight ! Total conductance (stomata+boundary layer) of the cohort, weighted by its leaf area [m/s]*[m2] real(r8) :: canopy_trim ! What is the fraction of the maximum leaf biomass that we are targeting? :- real(r8) :: leaf_cost ! How much does it cost to maintain leaves: kgC/m2/year-1 From 8dbc9f6a16ea981287bcbab71ad59d10a4892e14 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 22 Aug 2022 14:14:40 -0400 Subject: [PATCH 08/10] Removed unused smooth_leaf switch --- biogeochem/EDCanopyStructureMod.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 4066aec6b6..b3175202ba 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1493,7 +1493,6 @@ subroutine leaf_area_profile( currentSite ) real(r8) :: fraction_exposed ! how much of this layer is not covered by snow? real(r8) :: layer_top_hite ! notional top height of this canopy layer (m) real(r8) :: layer_bottom_hite ! notional bottom height of this canopy layer (m) - integer :: smooth_leaf_distribution ! is the leaf distribution this option (1) or not (0) real(r8) :: frac_canopy(N_HITE_BINS) ! amount of canopy in each height class real(r8) :: minh(N_HITE_BINS) ! minimum height in height class (m) real(r8) :: maxh(N_HITE_BINS) ! maximum height in height class (m) @@ -1504,7 +1503,6 @@ subroutine leaf_area_profile( currentSite ) !---------------------------------------------------------------------- - smooth_leaf_distribution = 0 ! Here we are trying to generate a profile of leaf area, indexed by 'z' and by pft ! We assume that each point in the canopy recieved the light attenuated by the average From bf99a0b6ba0901411875ce566c8a69d2b507c180 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 13 Sep 2022 10:07:37 -0600 Subject: [PATCH 09/10] use lai instead of treelai for b4b --- biogeochem/EDCanopyStructureMod.F90 | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index b3175202ba..3305d2d857 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1499,7 +1499,8 @@ subroutine leaf_area_profile( currentSite ) real(r8) :: dh ! vertical detph of height class (m) real(r8) :: min_chite ! bottom of cohort canopy (m) real(r8) :: max_chite ! top of cohort canopy (m) - real(r8) :: lai ! summed lai for checking m2 m-2 + real(r8) :: lai ! we use this to preserve b4b right now (m2/m2) + real(r8) :: sai ! "" !---------------------------------------------------------------------- @@ -1552,9 +1553,11 @@ subroutine leaf_area_profile( currentSite ) ! How much of each tree is stem area index? Assuming that there is ! This may indeed be zero if there is a sensecent grass ! ---------------------------------------------------------------- - + lai = currentCohort%treelai * currentCohort%c_area/currentPatch%total_canopy_area + sai = currentCohort%treesai * currentCohort%c_area/currentPatch%total_canopy_area if( (currentCohort%treelai+currentCohort%treesai) > nearzero)then - fleaf = currentCohort%treelai / (currentCohort%treelai + currentCohort%treesai) + !fleaf = currentCohort%treelai / (currentCohort%treelai + currentCohort%treesai) + fleaf = lai / (lai+sai) else fleaf = 0._r8 endif From cdb8fcae592db0ff6ccb557550cbe46ea5790714 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 13 Sep 2022 16:21:54 -0600 Subject: [PATCH 10/10] minor comment updates on lai usage during fleaf --- biogeochem/EDCanopyStructureMod.F90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 3305d2d857..0c508edce3 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1499,8 +1499,8 @@ subroutine leaf_area_profile( currentSite ) real(r8) :: dh ! vertical detph of height class (m) real(r8) :: min_chite ! bottom of cohort canopy (m) real(r8) :: max_chite ! top of cohort canopy (m) - real(r8) :: lai ! we use this to preserve b4b right now (m2/m2) - real(r8) :: sai ! "" + real(r8) :: lai ! leaf area per canopy area + real(r8) :: sai ! stem area per canopy area !---------------------------------------------------------------------- @@ -1556,7 +1556,9 @@ subroutine leaf_area_profile( currentSite ) lai = currentCohort%treelai * currentCohort%c_area/currentPatch%total_canopy_area sai = currentCohort%treesai * currentCohort%c_area/currentPatch%total_canopy_area if( (currentCohort%treelai+currentCohort%treesai) > nearzero)then - !fleaf = currentCohort%treelai / (currentCohort%treelai + currentCohort%treesai) + + ! See issue: https://github.com/NGEET/fates/issues/899 + ! fleaf = currentCohort%treelai / (currentCohort%treelai + currentCohort%treesai) fleaf = lai / (lai+sai) else fleaf = 0._r8