From 64d68b66d519dd3b24d150b77003a4cd8b9de153 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 9 Aug 2021 14:21:01 -0400 Subject: [PATCH 01/12] Changes that migrate the nlevleaf arrays to nlevleafmem, to reduce memory and computation --- biogeochem/EDCanopyStructureMod.F90 | 42 +-- biogeochem/EDCohortDynamicsMod.F90 | 149 +++++++++-- biogeochem/EDPhysiologyMod.F90 | 281 +++++++++++---------- biogeochem/FatesAllometryMod.F90 | 52 ++++ biogeophys/EDAccumulateFluxesMod.F90 | 14 +- biogeophys/FatesPlantRespPhotosynthMod.F90 | 79 +++--- main/EDTypesMod.F90 | 18 +- main/FatesHistoryInterfaceMod.F90 | 26 +- main/FatesRestartInterfaceMod.F90 | 26 +- 9 files changed, 457 insertions(+), 230 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 6aee6de8de..262d3238e9 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -14,12 +14,14 @@ module EDCanopyStructureMod use EDPftvarcon , only : EDPftvarcon_inst use PRTParametersMod , only : prt_params use FatesAllometryMod , only : carea_allom + use FatesAllometryMod , only : h2d_allom use EDCohortDynamicsMod , only : copy_cohort, terminate_cohorts, fuse_cohorts use EDCohortDynamicsMod , only : InitPRTObject use EDCohortDynamicsMod , only : InitPRTBoundaryConditions use EDCohortDynamicsMod , only : SendCohortToLitter use FatesAllometryMod , only : tree_lai use FatesAllometryMod , only : tree_sai + use FatesAllometryMod , only : GetNLevVeg use EDtypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use EDTypesMod , only : nclmax use EDTypesMod , only : nlevleaf @@ -1423,8 +1425,10 @@ subroutine leaf_area_profile( currentSite ) ! 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 + ! currentCohort%nveg_ ! The number of discrete vegetation ! ! layers needed to describe this crown + ! ! _max: on-allometry + ! ! _act: actual (not necessarily on allometry) ! ! The following patch level diagnostics are updated here: ! @@ -1479,9 +1483,13 @@ subroutine leaf_area_profile( currentSite ) real(r8) :: leaf_c ! leaf carbon [kg] !---------------------------------------------------------------------- - - - + real(r8) :: test_trim + real(r8) :: test_n + real(r8) :: test_vcmax25 + integer :: test_cl + real(r8) :: dbh + integer :: idbh + smooth_leaf_distribution = 0 ! Here we are trying to generate a profile of leaf area, indexed by 'z' and by pft @@ -1541,15 +1549,17 @@ subroutine leaf_area_profile( currentSite ) currentCohort%lai = currentCohort%treelai *currentCohort%c_area/currentPatch%total_canopy_area currentCohort%sai = currentCohort%treesai *currentCohort%c_area/currentPatch%total_canopy_area - ! Number of actual vegetation layers in this cohort's crown - currentCohort%nv = ceiling((currentCohort%treelai+currentCohort%treesai)/dinc_ed) - - currentPatch%ncan(cl,ft) = max(currentPatch%ncan(cl,ft),currentCohort%NV) - patch_lai = patch_lai + currentCohort%lai currentPatch%canopy_layer_tlai(cl) = currentPatch%canopy_layer_tlai(cl) + currentCohort%lai + + call GetNLevVeg(currentCohort%dbh, leaf_c, currentSite%spread, currentCohort%pft, & + currentCohort%canopy_trim, currentCohort%canopy_layer, & + currentPatch%canopy_layer_tlai, currentCohort%vcmax25top, & + currentCohort%nveg_act,currentCohort%nveg_max) + currentPatch%ncan(cl,ft) = max(currentPatch%ncan(cl,ft),currentCohort%nveg_act) + currentCohort => currentCohort%shorter enddo !currentCohort @@ -1681,7 +1691,7 @@ subroutine leaf_area_profile( currentSite ) ! before dividing it by the total area. Fill up layer for whole layers. ! -------------------------------------------------------------------------- - do iv = 1,currentCohort%NV + do iv = 1,currentCohort%nveg_act ! This loop builds the arrays that define the effective (not snow covered) ! and total (includes snow covered) area indices for leaves and stems @@ -1689,11 +1699,11 @@ subroutine leaf_area_profile( currentSite ) ! is obscured by snow. layer_top_hite = currentCohort%hite - & - ( real(iv-1,r8)/currentCohort%NV * currentCohort%hite * & + ( real(iv-1,r8)/currentCohort%nveg_act * currentCohort%hite * & EDPftvarcon_inst%crown(currentCohort%pft) ) layer_bottom_hite = currentCohort%hite - & - ( real(iv,r8)/currentCohort%NV * currentCohort%hite * & + ( real(iv,r8)/currentCohort%nveg_act * currentCohort%hite * & EDPftvarcon_inst%crown(currentCohort%pft) ) fraction_exposed = 1.0_r8 @@ -1709,13 +1719,13 @@ subroutine leaf_area_profile( currentSite ) (layer_top_hite-layer_bottom_hite )))) endif - if(iv==currentCohort%NV) then + if(iv==currentCohort%nveg_act) then remainder = (currentCohort%treelai + currentCohort%treesai) - & - (dinc_ed*real(currentCohort%nv-1,r8)) + (dinc_ed*real(currentCohort%nveg_act-1,r8)) if(remainder > dinc_ed )then write(fates_log(), *)'ED: issue with remainder', & currentCohort%treelai,currentCohort%treesai,dinc_ed, & - currentCohort%NV,remainder + currentCohort%nveg_act,remainder call endrun(msg=errMsg(sourcefile, __LINE__)) endif else @@ -1859,7 +1869,7 @@ subroutine leaf_area_profile( currentSite ) endif !leaf distribution end if - + currentPatch => currentPatch%younger enddo !patch diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index b6714ee3e9..5c529887c5 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -34,6 +34,7 @@ module EDCohortDynamicsMod use EDTypesMod , only : min_npm2, min_nppatch use EDTypesMod , only : min_n_safemath use EDTypesMod , only : nlevleaf + use EDTypesMod , only : nlevleafmem use PRTGenericMod , only : max_nleafage use EDTypesMod , only : ican_upper use EDTypesMod , only : site_fluxdiags_type @@ -65,7 +66,8 @@ module EDCohortDynamicsMod use FatesAllometryMod , only : carea_allom use FatesAllometryMod , only : ForceDBH use FatesAllometryMod , only : tree_lai, tree_sai - use FatesAllometryMod , only : set_root_fraction + use FatesAllometryMod , only : set_root_fraction + use FatesAllometryMod , only : GetNLevVeg use PRTGenericMod, only : prt_carbon_allom_hyp use PRTGenericMod, only : prt_cnp_flex_allom_hyp use PRTGenericMod, only : prt_vartypes @@ -300,8 +302,6 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, coage, dbh, & call InitPRTBoundaryConditions(new_cohort) - - ! Recuits do not have mortality rates, nor have they moved any ! carbon when they are created. They will bias our statistics ! until they have experienced a full day. We need a newly recruited flag. @@ -309,6 +309,12 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, coage, dbh, & ! growth, disturbance and mortality. new_cohort%isnew = .true. + ! New cohorts are not trimmable because they have not experienced + ! time to generate an annual carbon balance. This value will be + ! set to true to all existing cohorts at the end of the next trimming + ! attempt + new_cohort%is_trimmable = .false. + if( hlm_use_planthydro.eq.itrue ) then nlevrhiz = currentSite%si_hydr%nlevrhiz @@ -513,7 +519,8 @@ subroutine nan_cohort(cc_p) currentCohort%indexnumber = fates_unset_int ! unique number for each cohort. (within clump?) currentCohort%canopy_layer = fates_unset_int ! canopy status of cohort (1 = canopy, 2 = understorey, etc.) currentCohort%canopy_layer_yesterday = nan ! recent canopy status of cohort (1 = canopy, 2 = understorey, etc.) - currentCohort%NV = fates_unset_int ! Number of leaf layers: - + currentCohort%nveg_act = fates_unset_int + currentCohort%nveg_max = fates_unset_int currentCohort%status_coh = fates_unset_int ! growth status of plant (2 = leaves on , 1 = leaves off) currentCohort%size_class = fates_unset_int ! size class index currentCohort%size_class_lasttimestep = fates_unset_int ! size class index @@ -629,7 +636,6 @@ subroutine zero_cohort(cc_p) currentCohort => cc_p - currentCohort%NV = 0 currentCohort%status_coh = 0 currentCohort%rdark = 0._r8 currentCohort%resp_m = 0._r8 @@ -646,8 +652,6 @@ subroutine zero_cohort(cc_p) currentcohort%gpp_tstep = 0._r8 currentcohort%resp_tstep = 0._r8 currentcohort%resp_acc_hold = 0._r8 - - currentcohort%year_net_uptake(:) = 999._r8 ! this needs to be 999, or trimming of new cohorts will break. currentcohort%ts_net_uptake(:) = 0._r8 currentcohort%fraction_crown_burned = 0._r8 currentCohort%size_class = 1 @@ -1149,6 +1153,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) write(fates_log(),*) 'canopy_trim:',currentCohort%canopy_trim,nextc%canopy_trim write(fates_log(),*) 'canopy_layer_yesterday:', & currentCohort%canopy_layer_yesterday,nextc%canopy_layer_yesterday + write(fates_log(),*) 'is_trimmable:',currentCohort%is_trimmable,nextc%is_trimmable do i=1, nlevleaf write(fates_log(),*) 'leaf level: ',i,'year_net_uptake', & currentCohort%year_net_uptake(i),nextc%year_net_uptake(i) @@ -1183,8 +1188,10 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) 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 + + call FuseVegLayerMem(currentCohort,nextc,currentPatch,currentSite%spread) + + ! c13disc_acc calculation; weighted mean by GPP if ((currentCohort%n * currentCohort%gpp_acc + nextc%n * nextc%gpp_acc) .eq. 0.0_r8) then @@ -1434,15 +1441,10 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) currentCohort%ddbhdt = (currentCohort%n*currentCohort%ddbhdt + & nextc%n*nextc%ddbhdt)/newn - do i=1, nlevleaf - if (currentCohort%year_net_uptake(i) == 999._r8 .or. nextc%year_net_uptake(i) == 999._r8) then - currentCohort%year_net_uptake(i) = & - min(nextc%year_net_uptake(i),currentCohort%year_net_uptake(i)) - else - currentCohort%year_net_uptake(i) = (currentCohort%n*currentCohort%year_net_uptake(i) + & - nextc%n*nextc%year_net_uptake(i))/newn - endif - enddo + + + + end if !(currentCohort%isnew) @@ -1740,6 +1742,10 @@ subroutine insert_cohort(pcc, ptall, pshort, tnull, snull, storebigcohort, store end subroutine insert_cohort + + + + !-------------------------------------------------------------------------------------! subroutine copy_cohort( currentCohort,copyc ) ! @@ -1776,7 +1782,8 @@ subroutine copy_cohort( currentCohort,copyc ) n%leaf_cost = o%leaf_cost n%canopy_layer = o%canopy_layer n%canopy_layer_yesterday = o%canopy_layer_yesterday - n%nv = o%nv + n%nveg_act = o%nveg_act + n%nveg_max = o%nveg_max n%status_coh = o%status_coh n%canopy_trim = o%canopy_trim n%excl_weight = o%excl_weight @@ -1893,6 +1900,110 @@ subroutine copy_cohort( currentCohort,copyc ) end subroutine copy_cohort + ! ===================================================================================== + + subroutine FuseVegLayerMem(ccohort,ncohort,cpatch,site_spread) + + ! This routine handles the fusion of cohorts' vegetation + ! layer carbon balance memory, and its trimming value + ! This is tricky becuase: + ! + ! 1) two alike cohorts may or may not have existed since + ! the last time we started tracking net uptake + ! in which case we either revert to the one that did + ! or the typical weighted average + ! 2) two alike cohorts may not have had the same trimming + ! factor, and thus have been tracking carbon + ! balance over two different LAI depths + + type(ed_cohort_type) :: ccohort ! current (recipient) + type(ed_cohort_type) :: ncohort ! next (donor) + type(ed_patch_type) :: cpatch ! current patch (both cohorts are on this patch) + real(r8),intent(in) :: site_spread + + real(r8) :: joint_net_uptake(nlevleaf) + real(r8) :: joint_wgt(nlevleaf) + real(r8) :: leaf_c + integer :: iv ! veg layer index + integer :: ivm ! veg layer memory index + integer :: iv0,iv1 ! lowest and highest veg layer index matching memory + + joint_net_uptake(:) = 0._r8 + joint_wgt(:) = 0._r8 + + if(ccohort%is_trimmable .and. ncohort%is_trimmable) then + + ccohort%canopy_trim = (ccohort%n*ccohort%canopy_trim & + + ncohort%n*ncohort%canopy_trim)/(ccohort%n+ncohort%n) + + do ivm = 1, min(nlevleafmem,ccohort%nveg_max) + iv = ccohort%nveg_max - ivm + 1 + joint_net_uptake(iv) = ccohort%n*ccohort%year_net_uptake(ivm) + joint_wgt(iv) = ccohort%n + end do + do ivm = 1, min(nlevleafmem,ncohort%nveg_max) + iv = ncohort%nveg_max - ivm + 1 + joint_net_uptake(iv) = joint_net_uptake(iv)+ncohort%n*ncohort%year_net_uptake(ivm) + joint_wgt(iv) = joint_wgt(iv)+ncohort%n + end do + + iv0 = min(max(ccohort%nveg_max-nlevleafmem+1,1),max(ncohort%nveg_max-nlevleafmem+1,1)) + iv1 = max(ccohort%nveg_max,ncohort%nveg_max) + + if(debug)then + if(joint_wgt(iv0) currentSite%youngest_patch do while(associated(currentPatch)) - + ! Add debug diagnstic output to determine which patch if (debug) then write(fates_log(),*) 'Current patch:', ipatch write(fates_log(),*) 'Current patch cohorts:', currentPatch%countcohorts endif - + icohort = 1 - + currentCohort => currentPatch%tallest do while (associated(currentCohort)) - ! Save off the incoming trim and laimemory - initial_trim = currentCohort%canopy_trim - initial_laimem = currentCohort%laimemory + ! Save off the incoming trim and laimemory + initial_trim = currentCohort%canopy_trim + initial_laimem = currentCohort%laimemory ! 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 laimemory:', currentCohort%laimemory - endif + write(fates_log(),*) 'Current cohort:', icohort + write(fates_log(),*) 'Starting canopy trim:', initial_trim + write(fates_log(),*) 'Starting laimemory:', currentCohort%laimemory + endif trimmed = .false. ipft = currentCohort%pft - call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread,currentCohort%pft,currentCohort%c_area) - leaf_c = currentCohort%prt%GetState(leaf_organ, all_carbon_elements) + ! Identify current canopy layer (cl) + cl = currentCohort%canopy_layer + + leaf_c = currentCohort%prt%GetState(leaf_organ,carbon12_element) currentCohort%treelai = tree_lai(leaf_c, currentCohort%pft, currentCohort%c_area, & - currentCohort%n, currentCohort%canopy_layer, & - currentPatch%canopy_layer_tlai,currentCohort%vcmax25top ) + currentCohort%n, currentCohort%canopy_layer, & + currentPatch%canopy_layer_tlai,currentCohort%vcmax25top ) currentCohort%treesai = tree_sai(currentCohort%pft, currentCohort%dbh, currentCohort%canopy_trim, & - currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, & - currentPatch%canopy_layer_tlai, currentCohort%treelai, & - currentCohort%vcmax25top,0 ) - - currentCohort%nv = ceiling((currentCohort%treelai+currentCohort%treesai)/dinc_ed) - - if (currentCohort%nv > nlevleaf)then - write(fates_log(),*) 'nv > nlevleaf',currentCohort%nv, & - currentCohort%treelai,currentCohort%treesai, & - currentCohort%c_area,currentCohort%n,leaf_c + currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, & + currentPatch%canopy_layer_tlai, currentCohort%treelai , & + currentCohort%vcmax25top,5) + + call GetNLevVeg(currentCohort%dbh, leaf_c, currentSite%spread, currentCohort%pft, & + currentCohort%canopy_trim, currentCohort%canopy_layer, & + currentPatch%canopy_layer_tlai, currentCohort%vcmax25top, & + currentCohort%nveg_act,currentCohort%nveg_max) + + if (currentCohort%nveg_max > nlevleaf)then + write(fates_log(),*) 'nveg_max > nlevleaf',currentCohort%nveg_max, & + currentCohort%treelai,currentCohort%treesai, & + currentCohort%c_area,currentCohort%n,leaf_c call endrun(msg=errMsg(sourcefile, __LINE__)) endif @@ -486,48 +497,55 @@ subroutine trim_canopy( currentSite ) bfr_per_bleaf = tar_bfr/tar_bl endif - ! Identify current canopy layer (cl) - cl = currentCohort%canopy_layer - ! PFT-level maximum SLA value, even if under a thick canopy (same units as slatop) sla_max = prt_params%slamax(ipft) ! Initialize nnu_clai_a nnu_clai_a(:,:) = 0._r8 nnu_clai_b(:,:) = 0._r8 - - !Leaf cost vs netuptake for each leaf layer. - do z = 1, currentCohort%nv - ! Calculate the cumulative total vegetation area index (no snow occlusion, stems and leaves) + ! We dont trimm, if it is not trimmable - leaf_inc = dinc_ed * & - currentCohort%treelai/(currentCohort%treelai+currentCohort%treesai) - - ! Now calculate the cumulative top-down lai of the current layer's midpoint within the current cohort - lai_layers_above = leaf_inc * (z-1) - lai_current = min(leaf_inc, currentCohort%treelai - lai_layers_above) - cumulative_lai_cohort = lai_layers_above + 0.5*lai_current + if_trimmable: if(currentCohort%is_trimmable)then - ! Now add in the lai above the current cohort for calculating the sla leaf level - lai_canopy_above = sum(currentPatch%canopy_layer_tlai(1:cl-1)) - cumulative_lai = lai_canopy_above + cumulative_lai_cohort + ! Leaf cost vs netuptake for each leaf layer. + ! Lets loop over the last few leaf layers, only the number of layers + ! necessary for evaluating the trend + + leafmem_loop: do zm = 1, min(nll,currentCohort%nveg_max) + + ! The actual leaf layer index (going from top to bottom) + z = currentCohort%nveg_max - zm + 1 + + ! Calculate the cumulative total vegetation area index (no snow occlusion, stems and leaves) + + leaf_inc = dinc_ed * & + currentCohort%treelai/(currentCohort%treelai+currentCohort%treesai) + + ! Now calculate the cumulative top-down lai of the current layer's midpoint within the current cohort + lai_layers_above = leaf_inc * (z-1) + lai_current = min(leaf_inc, currentCohort%treelai - lai_layers_above) + cumulative_lai_cohort = lai_layers_above + 0.5*lai_current + + ! Now add in the lai above the current cohort for calculating the sla leaf level + lai_canopy_above = sum(currentPatch%canopy_layer_tlai(1:cl-1)) + cumulative_lai = lai_canopy_above + cumulative_lai_cohort - ! There was activity this year in this leaf layer. This should only occur for bottom most leaf layer - if (currentCohort%year_net_uptake(z) /= 999._r8)then - ! Calculate sla_levleaf following the sla profile with overlying leaf area ! Scale for leaf nitrogen profile kn = decay_coeff_kn(ipft,currentCohort%vcmax25top) ! Nscaler value at leaf level z nscaler_levleaf = exp(-kn * cumulative_lai) + ! Sla value at leaf level z after nitrogen profile scaling (m2/gC) + ! (RGK) THIS SHOULD BE USING THE SLAMAX AS WELL< NOT JUST A CAP sla_levleaf = prt_params%slatop(ipft)/nscaler_levleaf + ! if(sla_levleaf > sla_max)then sla_levleaf = sla_max end if - + !Leaf Cost kgC/m2/year-1 !decidous costs. if (prt_params%season_decid(ipft) == itrue .or. & @@ -544,15 +562,16 @@ subroutine trim_canopy( currentSite ) bfr_per_bleaf / prt_params%root_long(ipft) endif + currentCohort%leaf_cost = currentCohort%leaf_cost * & - (prt_params%grperc(ipft) + 1._r8) + (prt_params%grperc(ipft) + 1._r8) else !evergreen costs ! Leaf cost at leaf level z accounting for sla profile currentCohort%leaf_cost = 1.0_r8/(sla_levleaf* & sum(prt_params%leaf_long(ipft,:))*1000.0_r8) !convert from sla in m2g-1 to m2kg-1 - - + + if ( int(prt_params%allom_fmode(ipft)) .eq. 1 ) then ! if using trimmed leaf for fine root biomass allometry, add the cost of the root increment ! to the leaf increment; otherwise do not. @@ -561,45 +580,39 @@ subroutine trim_canopy( currentSite ) bfr_per_bleaf / prt_params%root_long(ipft) endif currentCohort%leaf_cost = currentCohort%leaf_cost * & - (prt_params%grperc(ipft) + 1._r8) + (prt_params%grperc(ipft) + 1._r8) endif ! Construct the arrays for a least square fit of the net_net_uptake versus the cumulative lai ! if at least nll leaf layers are present in the current cohort and only for the bottom nll ! leaf layers. - if (currentCohort%nv > nll .and. currentCohort%nv - z < nll) then - - ! Build the A matrix for the LHS of the linear system. A = [n sum(x); sum(x) sum(x^2)] - ! where n = nll and x = yearly_net_uptake-leafcost - nnu_clai_a(1,1) = nnu_clai_a(1,1) + 1 ! Increment for each layer used - nnu_clai_a(1,2) = nnu_clai_a(1,2) + currentCohort%year_net_uptake(z) - currentCohort%leaf_cost - nnu_clai_a(2,1) = nnu_clai_a(1,2) - nnu_clai_a(2,2) = nnu_clai_a(2,2) + (currentCohort%year_net_uptake(z) - currentCohort%leaf_cost)**2 - - ! Build the B matrix for the RHS of the linear system. B = [sum(y); sum(x*y)] - ! where x = yearly_net_uptake-leafcost and y = cumulative_lai_cohort - nnu_clai_b(1,1) = nnu_clai_b(1,1) + cumulative_lai_cohort - nnu_clai_b(2,1) = nnu_clai_b(2,1) + (cumulative_lai_cohort * & - (currentCohort%year_net_uptake(z) - currentCohort%leaf_cost)) - end if + + ! Build the A matrix for the LHS of the linear system. A = [n sum(x); sum(x) sum(x^2)] + ! where n = nll and x = yearly_net_uptake-leafcost + nnu_clai_a(1,1) = nnu_clai_a(1,1) + 1 ! Increment for each layer used + nnu_clai_a(1,2) = nnu_clai_a(1,2) + currentCohort%year_net_uptake(zm) - currentCohort%leaf_cost + nnu_clai_a(2,1) = nnu_clai_a(1,2) + nnu_clai_a(2,2) = nnu_clai_a(2,2) + (currentCohort%year_net_uptake(zm) - currentCohort%leaf_cost)**2 + + ! Build the B matrix for the RHS of the linear system. B = [sum(y); sum(x*y)] + ! where x = yearly_net_uptake-leafcost and y = cumulative_lai_cohort + nnu_clai_b(1,1) = nnu_clai_b(1,1) + cumulative_lai_cohort + nnu_clai_b(2,1) = nnu_clai_b(2,1) + (cumulative_lai_cohort * & + (currentCohort%year_net_uptake(zm) - currentCohort%leaf_cost)) + ! Check leaf cost against the yearly net uptake for that cohort leaf layer - if (currentCohort%year_net_uptake(z) < currentCohort%leaf_cost) then + if (currentCohort%year_net_uptake(zm) < currentCohort%leaf_cost) then ! Make sure the cohort trim fraction is great than the pft trim limit if (currentCohort%canopy_trim > EDPftvarcon_inst%trim_limit(ipft)) then - ! if ( debug ) then - ! write(fates_log(),*) 'trimming leaves', & - ! currentCohort%canopy_trim,currentCohort%leaf_cost - ! endif - ! keep trimming until none of the canopy is in negative carbon balance. if (currentCohort%hite > EDPftvarcon_inst%hgt_min(ipft)) then currentCohort%canopy_trim = currentCohort%canopy_trim - & - EDPftvarcon_inst%trim_inc(ipft) + EDPftvarcon_inst%trim_inc(ipft) if (prt_params%evergreen(ipft) /= 1)then currentCohort%laimemory = currentCohort%laimemory * & - (1.0_r8 - EDPftvarcon_inst%trim_inc(ipft)) + (1.0_r8 - EDPftvarcon_inst%trim_inc(ipft)) endif trimmed = .true. @@ -607,74 +620,78 @@ subroutine trim_canopy( currentSite ) endif ! hite check endif ! trim limit check endif ! net uptake check - endif ! leaf activity check - enddo ! z, leaf layer loop + + enddo leafmem_loop + end if if_trimmable ! Compute the optimal cumulative lai based on the cohort net-net uptake profile if at least 2 leaf layers - if (nnu_clai_a(1,1) > 1) then - - ! Compute the optimum size of the work array - lwork = -1 ! Ask sgels to compute optimal number of entries for work - call dgels(trans, m, n, nrhs, nnu_clai_a, lda, nnu_clai_b, ldb, work, lwork, info) - lwork = int(work(1)) ! Pick the optimum. TBD, can work(1) come back with greater than work size? - - ! if (debug) then - ! write(fates_log(),*) 'LLSF lwork output (info, lwork):', info, lwork - ! endif - - ! Compute the minimum of 2-norm of of the least squares fit to solve for X - ! Note that dgels returns the solution by overwriting the nnu_clai_b array. - ! The result has the form: X = [b; m] - ! where b = y-intercept (i.e. the cohort lai that has zero yearly net-net uptake) - ! and m is the slope of the linear fit - call dgels(trans, m, n, nrhs, nnu_clai_a, lda, nnu_clai_b, ldb, work, lwork, info) - - if (info < 0) then - write(fates_log(),*) 'LLSF optimium LAI calculation returned illegal value' - call endrun(msg=errMsg(sourcefile, __LINE__)) - endif - - if (debug) then - write(fates_log(),*) 'LLSF optimium LAI (intercept,slope):', nnu_clai_b - write(fates_log(),*) 'LLSF optimium LAI:', nnu_clai_b(1,1) - write(fates_log(),*) 'LLSF optimium LAI info:', info - write(fates_log(),*) 'LAI fraction (optimum_lai/cumulative_lai):', nnu_clai_b(1,1) / cumulative_lai_cohort - endif - - ! Calculate the optimum trim based on the initial canopy trim value - if (cumulative_lai_cohort > 0._r8) then ! Sometime cumulative_lai comes in at 0.0? - - ! - optimum_trim = (nnu_clai_b(1,1) / cumulative_lai_cohort) * initial_trim - optimum_laimem = (nnu_clai_b(1,1) / cumulative_lai_cohort) * initial_laimem - - ! 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 laimemory as well - if (prt_params%evergreen(ipft) /= 1) then - currentCohort%laimemory = optimum_laimem - endif + if_optimumtrim: if (nnu_clai_a(1,1) > 1) then + + ! Compute the optimum size of the work array + lwork = -1 ! Ask sgels to compute optimal number of entries for work + call dgels(trans, m, n, nrhs, nnu_clai_a, lda, nnu_clai_b, ldb, work, lwork, info) + lwork = int(work(1)) ! Pick the optimum. TBD, can work(1) come back with greater than work size? + + ! if (debug) then + ! write(fates_log(),*) 'LLSF lwork output (info, lwork):', info, lwork + ! endif + + ! Compute the minimum of 2-norm of of the least squares fit to solve for X + ! Note that dgels returns the solution by overwriting the nnu_clai_b array. + ! The result has the form: X = [b; m] + ! where b = y-intercept (i.e. the cohort lai that has zero yearly net-net uptake) + ! and m is the slope of the linear fit + call dgels(trans, m, n, nrhs, nnu_clai_a, lda, nnu_clai_b, ldb, work, lwork, info) + + if (info < 0) then + write(fates_log(),*) 'LLSF optimium LAI calculation returned illegal value' + call endrun(msg=errMsg(sourcefile, __LINE__)) + endif + + if (debug) then + write(fates_log(),*) 'LLSF optimium LAI (intercept,slope):', nnu_clai_b + write(fates_log(),*) 'LLSF optimium LAI:', nnu_clai_b(1,1) + write(fates_log(),*) 'LLSF optimium LAI info:', info + write(fates_log(),*) 'LAI fraction (optimum_lai/cumulative_lai):', nnu_clai_b(1,1) / cumulative_lai_cohort + endif + + ! Calculate the optimum trim based on the initial canopy trim value + if (cumulative_lai_cohort > 0._r8) then ! Sometime cumulative_lai comes in at 0.0? + + ! + optimum_trim = (nnu_clai_b(1,1) / cumulative_lai_cohort) * initial_trim + optimum_laimem = (nnu_clai_b(1,1) / cumulative_lai_cohort) * initial_laimem + + ! 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 laimemory as well + if (prt_params%evergreen(ipft) /= 1) then + currentCohort%laimemory = optimum_laimem + endif - trimmed = .true. + trimmed = .true. - endif - endif - endif + endif + endif + endif if_optimumtrim - ! Reset activity for the cohort for the start of the next year - currentCohort%year_net_uptake(:) = 999.0_r8 + ! Set cohort annual uptake to zero, this also + ! initializes a cohort so that the mean should proceed. Cohorts + ! that are not zero'd here, are not allowed to trim (because + ! they need full year) + currentCohort%year_net_uptake(:) = 0.0_r8 ! Add to trim fraction if cohort not trimmed at all if ( (.not.trimmed) .and.currentCohort%canopy_trim < 1.0_r8)then currentCohort%canopy_trim = currentCohort%canopy_trim + EDPftvarcon_inst%trim_inc(ipft) - endif + endif if ( debug ) then write(fates_log(),*) 'trimming:',currentCohort%canopy_trim endif - + ! currentCohort%canopy_trim = 1.0_r8 !FIX(RF,032414) this turns off ctrim for now. currentCohort => currentCohort%shorter icohort = icohort + 1 diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 8e27faae22..3c0824afbf 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -102,6 +102,7 @@ module FatesAllometryMod implicit none private + public :: GetNLevVeg ! Number of vegetation layers for the plant public :: h2d_allom ! Generic height to diameter wrapper public :: h_allom ! Generic diameter to height wrapper public :: bagw_allom ! Generic AGWB (above grnd. woody bio) wrapper @@ -758,6 +759,57 @@ real(r8) function tree_sai( pft, dbh, canopy_trim, c_area, nplant, cl, & return end function tree_sai + + + subroutine GetNLevVeg(dbh, leaf_c, site_spread, ipft, canopy_trim, & + canopy_layer,canopy_layer_tlai, vcmax25top, & + nveg_act, nveg_max) + + ! Determine the maximum number of leaf (vegetation) layers + ! for a cohort. This is based off of allometry, and assuming the plant + ! has the maximum leaf for its size and trimming. + + real(r8), intent(inout) :: dbh + real(r8), intent(in) :: leaf_c + real(r8), intent(in) :: site_spread + integer, intent(in) :: ipft + real(r8), intent(in) :: canopy_trim + integer, intent(in) :: canopy_layer + real(r8), intent(in) :: canopy_layer_tlai(:) + real(r8), intent(in) :: vcmax25top + integer, intent(out) :: nveg_act + integer, intent(out) :: nveg_max + + integer :: nv + real(r8) :: c_area + real(r8) :: leaf_c_target + real(r8) :: max_lai, max_sai + real(r8), parameter :: nplant_dum = 1.0_r8 ! Number of plants falls out (dummy) + + call carea_allom(dbh,nplant_dum,site_spread,ipft,c_area) + + max_lai = tree_lai(leaf_c, ipft, c_area, & + nplant_dum, canopy_layer, canopy_layer_tlai, vcmax25top ) + + max_sai = tree_sai(ipft, dbh, canopy_trim, & + c_area, nplant_dum, canopy_layer, canopy_layer_tlai, max_lai, & + vcmax25top, 0) + + nveg_act = ceiling((max_lai+max_sai)/dinc_ed) + + call bleaf(dbh,ipft,canopy_trim,leaf_c_target) + + max_lai = tree_lai(leaf_c_target, ipft, c_area, & + nplant_dum, canopy_layer, canopy_layer_tlai, vcmax25top ) + + max_sai = tree_sai(ipft, dbh, canopy_trim, & + c_area, nplant_dum, canopy_layer, canopy_layer_tlai, max_lai, & + vcmax25top, 0) + + nveg_max = ceiling((max_lai+max_sai)/dinc_ed) + + return + end subroutine GetNLevVeg ! ============================================================================ ! Generic sapwood biomass interface diff --git a/biogeophys/EDAccumulateFluxesMod.F90 b/biogeophys/EDAccumulateFluxesMod.F90 index 4d873cca85..80afeba6c3 100644 --- a/biogeophys/EDAccumulateFluxesMod.F90 +++ b/biogeophys/EDAccumulateFluxesMod.F90 @@ -13,6 +13,7 @@ module EDAccumulateFluxesMod use FatesGlobals, only : fates_log use shr_log_mod , only : errMsg => shr_log_errMsg use FatesConstantsMod , only : r8 => fates_r8 + use EDTypesMod, only : nlevleafmem implicit none @@ -93,13 +94,12 @@ subroutine AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) (ccohort%c13disc_clm * ccohort%gpp_tstep)) / & (ccohort%gpp_acc + ccohort%gpp_tstep) endif - - do iv=1,ccohort%nv - if(ccohort%year_net_uptake(iv) == 999._r8)then ! note that there were leaves in this layer this year. - ccohort%year_net_uptake(iv) = 0._r8 - end if - ccohort%year_net_uptake(iv) = ccohort%year_net_uptake(iv) + ccohort%ts_net_uptake(iv) - enddo + + if( ccohort%is_trimmable ) then + do iv=1,min(ccohort%nveg_max,nlevleafmem) + ccohort%year_net_uptake(iv) = ccohort%year_net_uptake(iv) + ccohort%ts_net_uptake(iv) + end do + end if ccohort => ccohort%taller enddo ! while(associated(ccohort)) diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 349f6473bf..25e32eaeba 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -32,6 +32,7 @@ module FATESPlantRespPhotosynthMod use FatesInterfaceTypesMod, only : nleafage use EDTypesMod, only : maxpft use EDTypesMod, only : nlevleaf + use EDTypesMod, only : nlevleafmem use EDTypesMod, only : nclmax use PRTGenericMod, only : max_nleafage use EDTypesMod, only : do_fates_salinity @@ -214,7 +215,8 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! above the leaf layer of interest real(r8) :: lai_current ! the LAI in the current leaf layer real(r8) :: cumulative_lai ! the cumulative LAI, top down, to the leaf layer of interest - + integer :: ivm, iv, nvm, nva ! Loop indices and bounds for cohort vegetation layers + real(r8), allocatable :: rootfr_ft(:,:) ! Root fractions per depth and PFT ! ----------------------------------------------------------------------------------- @@ -224,7 +226,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) !real(r8) :: psncanopy_pa ! patch sunlit leaf photosynthesis (umol CO2 /m**2/ s) !real(r8) :: lmrcanopy_pa ! patch sunlit leaf maintenance respiration rate (umol CO2/m**2/s) - integer :: cl,s,iv,j,ps,ft,ifp ! indices + integer :: cl,s,j,ps,ft,ifp ! indices integer :: nv ! number of leaf layers integer :: NCL_p ! number of canopy layers in patch integer :: iage ! loop counter for leaf age classes @@ -377,7 +379,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) if(currentPatch%canopy_mask(cl,ft) == 1)then ! Loop over leaf-layers - do iv = 1,currentCohort%nv + do iv = 1,currentCohort%nveg_act ! ------------------------------------------------------------ ! If we are doing plant hydro-dynamics (or any run-type @@ -560,13 +562,14 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! (where the rates become per cohort, ie /individual). Most likely ! a sum over layers. ! --------------------------------------------------------------- - nv = currentCohort%nv + nva = currentCohort%nveg_act + nvm = currentCohort%nveg_max call ScaleLeafLayerFluxToCohort(nv, & !in - currentPatch%psn_z(cl,ft,1:nv), & !in - lmr_z(1:nv,ft,cl), & !in - rs_z(1:nv,ft,cl), & !in - currentPatch%elai_profile(cl,ft,1:nv), & !in - c13disc_z(cl, ft, 1:nv), & !in + currentPatch%psn_z(cl,ft,1:nva), & !in + lmr_z(1:nva,ft,cl), & !in + rs_z(1:nva,ft,cl), & !in + currentPatch%elai_profile(cl,ft,1:nva),& !in + c13disc_z(cl, ft, 1:nva), & !in currentCohort%c_area, & !in currentCohort%n, & !in bc_in(s)%rb_pa(ifp), & !in @@ -578,8 +581,20 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) cohort_eleaf_area) !out ! Net Uptake does not need to be scaled, just transfer directly - currentCohort%ts_net_uptake(1:nv) = anet_av_z(1:nv,ft,cl) * umolC_to_kgC - + ! only track the last few bins + ! Also, here we track net uptake in the vegetation memory + ! array. This has a reverse indexing. 1 refers to the lowest possible + ! bin in the cohort's vegetation layering scheme, which is defined by + ! allometry. The cohort might not actually have leaves down to its + ! maximum depth because of a weak carbon balance or due to + ! being in a deciduous transition. + if(nva>0) then + do ivm = nvm-nva+1, min(nlevleafmem,nva) + iv = nvm - ivm + 1 + currentCohort%ts_net_uptake(ivm) = anet_av_z(iv,ft,cl) * umolC_to_kgC + end do + end if + else ! In this case, the cohort had no leaves, @@ -1245,12 +1260,12 @@ end subroutine LeafLayerPhotosynthesis ! ===================================================================================== - subroutine ScaleLeafLayerFluxToCohort(nv, & ! in currentCohort%nv - psn_llz, & ! in %psn_z(1:currentCohort%nv,ft,cl) - lmr_llz, & ! in lmr_z(1:currentCohort%nv,ft,cl) - rs_llz, & ! in rs_z(1:currentCohort%nv,ft,cl) - elai_llz, & ! in %elai_profile(cl,ft,1:currentCohort%nv) - c13disc_llz, & ! in c13disc_z(cl, ft, 1:currentCohort%nv) + subroutine ScaleLeafLayerFluxToCohort(nveg_act, & ! in nveg_act + psn_llz, & ! in %psn_z(nveg_act,ft,cl) + lmr_llz, & ! in lmr_z(nveg_act,ft,cl) + rs_llz, & ! in rs_z(nveg_act,ft,cl) + elai_llz, & ! in %elai_profile(cl,ft,nveg_act) + c13disc_llz, & ! in c13disc_z(cl, ft, 1:nveg_act) c_area, & ! in currentCohort%c_area nplant, & ! in currentCohort%n rb, & ! in bc_in(s)%rb_pa(ifp) @@ -1272,12 +1287,12 @@ subroutine ScaleLeafLayerFluxToCohort(nv, & ! in currentCohort%nv use FatesConstantsMod, only : umolC_to_kgC ! Arguments - integer, intent(in) :: nv ! number of active leaf layers - real(r8), intent(in) :: psn_llz(nv) ! layer photosynthesis rate (GPP) [umolC/m2leaf/s] - real(r8), intent(in) :: lmr_llz(nv) ! layer dark respiration rate [umolC/m2leaf/s] - real(r8), intent(in) :: rs_llz(nv) ! leaf layer stomatal resistance [s/m] - real(r8), intent(in) :: elai_llz(nv) ! exposed LAI per layer [m2 leaf/ m2 pft footprint] - real(r8), intent(in) :: c13disc_llz(nv) ! leaf layer c13 discrimination, weighted mean + integer, intent(in) :: nveg_act ! number of active leaf layers + real(r8), intent(in) :: psn_llz(nveg_act) ! layer photosynthesis rate (GPP) [umolC/m2leaf/s] + real(r8), intent(in) :: lmr_llz(nveg_act) ! layer dark respiration rate [umolC/m2leaf/s] + real(r8), intent(in) :: rs_llz(nveg_act) ! leaf layer stomatal resistance [s/m] + real(r8), intent(in) :: elai_llz(nveg_act) ! exposed LAI per layer [m2 leaf/ m2 pft footprint] + real(r8), intent(in) :: c13disc_llz(nveg_act) ! leaf layer c13 discrimination, weighted mean real(r8), intent(in) :: c_area ! crown area m2/m2 real(r8), intent(in) :: nplant ! indiv/m2 real(r8), intent(in) :: rb ! leaf boundary layer resistance (s/m) @@ -1303,7 +1318,7 @@ subroutine ScaleLeafLayerFluxToCohort(nv, & ! in currentCohort%nv gpp = 0.0_r8 rdark = 0.0_r8 - do il = 1, nv ! Loop over the leaf layers this cohort participates in + do il = 1, nveg_act ! Loop over the leaf layers this cohort participates in ! Cohort's total effective leaf area in this layer [m2] @@ -1334,13 +1349,13 @@ subroutine ScaleLeafLayerFluxToCohort(nv, & ! in currentCohort%nv - if (nv > 1) then + if (nveg_act > 1) then ! cohort%c13disc_clm as weighted mean of d13c flux at all related leave layers - sum_weight = sum(psn_llz(1:nv-1) * elai_llz(1:nv-1)) + sum_weight = sum(psn_llz(1:nveg_act-1) * elai_llz(1:nveg_act-1)) if (sum_weight .eq. 0.0_r8) then c13disc_clm = 0.0 else - c13disc_clm = sum(c13disc_llz(1:nv-1) * psn_llz(1:nv-1) * elai_llz(1:nv-1)) / sum_weight + c13disc_clm = sum(c13disc_llz(1:nveg_act-1) * psn_llz(1:nveg_act-1) * elai_llz(1:nveg_act-1)) / sum_weight end if end if @@ -1363,11 +1378,11 @@ subroutine ScaleLeafLayerFluxToCohort(nv, & ! in currentCohort%nv if ( debug ) then write(fates_log(),*) 'EDPhoto 816 ', gpp - write(fates_log(),*) 'EDPhoto 817 ', psn_llz(1:nv) - write(fates_log(),*) 'EDPhoto 820 ', nv - write(fates_log(),*) 'EDPhoto 821 ', elai_llz(1:nv) + write(fates_log(),*) 'EDPhoto 817 ', psn_llz(1:nveg_act) + write(fates_log(),*) 'EDPhoto 820 ', nveg_act + write(fates_log(),*) 'EDPhoto 821 ', elai_llz(1:nveg_act) write(fates_log(),*) 'EDPhoto 843 ', rdark - write(fates_log(),*) 'EDPhoto 873 ', nv + write(fates_log(),*) 'EDPhoto 873 ', nveg_act write(fates_log(),*) 'EDPhoto 874 ', cohort_eleaf_area endif @@ -1602,7 +1617,7 @@ subroutine UpdateCanopyNCanNRadPresent(currentPatch) currentPatch%ncan(currentCohort%canopy_layer,currentCohort%pft) = & max(currentPatch%ncan(currentCohort%canopy_layer,currentCohort%pft), & - currentCohort%NV) + currentCohort%nveg_act) currentCohort => currentCohort%shorter diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 5da7babc54..e6fe12ca13 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -36,6 +36,7 @@ module EDTypesMod ! to understory layers (all layers that ! are not the top canopy layer) + integer, parameter, public :: nlevleafmem = 4 ! number of layers for veg memory used for trimming integer, parameter, public :: nlevleaf = 30 ! number of leaf layers in canopy layer integer, parameter, public :: maxpft = 16 ! maximum number of PFTs allowed ! the parameter file may determine that fewer @@ -231,13 +232,20 @@ module EDTypesMod real(r8) :: leaf_cost ! How much does it cost to maintain leaves: kgC/m2/year-1 real(r8) :: excl_weight ! How much of this cohort is demoted each year, as a proportion of all cohorts:- real(r8) :: prom_weight ! How much of this cohort is promoted each year, as a proportion of all cohorts:- - integer :: nv ! Number of leaf layers: - + integer :: nveg_act ! Number of vegetation layers (actual) at any point in time + integer :: nveg_max ! Maximum number of vegetation layers based on allometry integer :: status_coh ! growth status of plant (2 = leaves on , 1 = leaves off) real(r8) :: c_area ! areal extent of canopy (m2) real(r8) :: treelai ! lai of an individual within cohort leaf area (m2) / crown area (m2) real(r8) :: treesai ! stem area index of an indiv. within cohort: stem area (m2) / crown area (m2) logical :: isnew ! flag to signify a new cohort, new cohorts have not experienced ! npp or mortality and should therefore not be fused or averaged + logical :: is_trimmable ! This is only true if a cohort has been in existance since the last trimming attempt. + ! this is necessary because we require a whole year's worth of carbon balance + ! information to make the decision. Newly recruited cohorts (which are not trimmable), that + ! are fused into a trimmable cohort, will inherit the trimmable flag and the annual + ! carbon balance information. All cohorts should be initialized as is_trimmable = .false. + integer :: size_class ! An index that indicates which diameter size bin the cohort currently resides in ! this is used for history output. We maintain this in the main cohort memory ! because we don't want to continually re-calculate the cohort's position when @@ -322,8 +330,8 @@ module EDTypesMod - real(r8) :: ts_net_uptake(nlevleaf) ! Net uptake of leaf layers: kgC/m2/timestep - real(r8) :: year_net_uptake(nlevleaf) ! Net uptake of leaf layers: kgC/m2/year + real(r8) :: ts_net_uptake(nlevleafmem) ! Net uptake of leaf layers: kgC/m2/timestep + real(r8) :: year_net_uptake(nlevleafmem) ! Net uptake of leaf layers: kgC/m2/year ! RESPIRATION COMPONENTS real(r8) :: rdark ! Dark respiration: kgC/indiv/s @@ -1042,7 +1050,8 @@ subroutine dump_cohort(ccohort) write(fates_log(),*) 'co%leaf_cost = ', ccohort%leaf_cost write(fates_log(),*) 'co%canopy_layer = ', ccohort%canopy_layer write(fates_log(),*) 'co%canopy_layer_yesterday = ', ccohort%canopy_layer_yesterday - write(fates_log(),*) 'co%nv = ', ccohort%nv + write(fates_log(),*) 'co%nveg_act = ', ccohort%nveg_act + write(fates_log(),*) 'co%nveg_max = ', ccohort%nveg_max write(fates_log(),*) 'co%status_coh = ', ccohort%status_coh write(fates_log(),*) 'co%canopy_trim = ', ccohort%canopy_trim write(fates_log(),*) 'co%excl_weight = ', ccohort%excl_weight @@ -1079,6 +1088,7 @@ subroutine dump_cohort(ccohort) write(fates_log(),*) 'co%frmort = ', ccohort%frmort write(fates_log(),*) 'co%asmort = ', ccohort%asmort write(fates_log(),*) 'co%isnew = ', ccohort%isnew + write(fates_log(),*) 'co%is_trimmable = ', ccohort%is_trimmable write(fates_log(),*) 'co%dndt = ', ccohort%dndt write(fates_log(),*) 'co%dhdt = ', ccohort%dhdt write(fates_log(),*) 'co%ddbhdt = ', ccohort%ddbhdt diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 97f3342b43..a61acebc7f 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -2779,7 +2779,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! hio_crownarea_si_can(io_si, ican) = hio_crownarea_si_can(io_si, ican) + ccohort%c_area / AREA ! - do ileaf=1,ccohort%nv + do ileaf=1,ccohort%nveg_act cnlf_indx = ileaf + (ican-1) * nlevleaf hio_crownarea_si_cnlf(io_si, cnlf_indx) = hio_crownarea_si_cnlf(io_si, cnlf_indx) + & ccohort%c_area / AREA @@ -3669,13 +3669,13 @@ subroutine update_history_hifrq(this,nc,nsites,sites,bc_in,dt_tstep) end associate endif - !!! canopy leaf carbon balance - ican = ccohort%canopy_layer - do ileaf=1,ccohort%nv - cnlf_indx = ileaf + (ican-1) * nlevleaf - hio_ts_net_uptake_si_cnlf(io_si, cnlf_indx) = hio_ts_net_uptake_si_cnlf(io_si, cnlf_indx) + & - ccohort%ts_net_uptake(ileaf) * g_per_kg * per_dt_tstep * ccohort%c_area / AREA - end do +! !!! canopy leaf carbon balance +! ican = ccohort%canopy_layer +! do ileaf=1,ccohort%nveg_act +! cnlf_indx = ileaf + (ican-1) * nlevleaf +! hio_ts_net_uptake_si_cnlf(io_si, cnlf_indx) = hio_ts_net_uptake_si_cnlf(io_si, cnlf_indx) + & +! ccohort%ts_net_uptake(ileaf) * g_per_kg * per_dt_tstep * ccohort%c_area / AREA +! end do ccohort => ccohort%taller enddo ! cohort loop @@ -5156,11 +5156,11 @@ subroutine define_history_vars(this, initialize_variables) ivar=ivar, initialize=initialize_variables, index = ih_fabi_sha_top_si_can ) !!! canopy-resolved fluxes and structure - call this%set_history_var(vname='NET_C_UPTAKE_CNLF', units='gC/m2/s', & - long='net carbon uptake by each canopy and leaf layer per unit ground area (i.e. divide by CROWNAREA_CNLF to make per leaf area)', & - use_default='inactive', & - avgflag='A', vtype=site_cnlf_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=2, & - ivar=ivar, initialize=initialize_variables, index = ih_ts_net_uptake_si_cnlf ) + !!call this%set_history_var(vname='NET_C_UPTAKE_CNLF', units='gC/m2/s', & + !! long='net carbon uptake by each canopy and leaf layer per unit ground area (i.e. divide by CROWNAREA_CNLF to make per leaf area)', & + !! use_default='inactive', & + !! avgflag='A', vtype=site_cnlf_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=2, & + !! ivar=ivar, initialize=initialize_variables, index = ih_ts_net_uptake_si_cnlf ) call this%set_history_var(vname='CROWNAREA_CNLF', units='m2/m2', & long='total crown area that is occupied by leaves in each canopy and leaf layer', & diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 7ae00ed0b2..cfc74672a7 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -143,6 +143,7 @@ module FatesRestartInterfaceMod integer :: ir_pft_co integer :: ir_status_co integer :: ir_isnew_co + integer :: ir_istrimmable_co ! Litter integer :: ir_agcwd_litt @@ -219,9 +220,7 @@ module FatesRestartInterfaceMod integer, parameter, public :: fates_restart_num_dimensions = 2 !(cohort,column) integer, parameter, public :: fates_restart_num_dim_kinds = 4 !(cohort-int,cohort-r8,site-int,site-r8) - ! integer constants for storing logical data - integer, parameter, public :: old_cohort = 0 - integer, parameter, public :: new_cohort = 1 + real(r8), parameter, public :: flushinvalid = -9999.0 real(r8), parameter, public :: flushzero = 0.0 @@ -864,6 +863,11 @@ subroutine define_restart_vars(this, initialize_variables) units='0/1', flushval = flushone, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_isnew_co ) + call this%set_restart_var(vname='fates_istrimmable', vtype=cohort_int, & + long_name='ed cohort - binary flag specifying if the cohort has persisted since the last trim attempt', & + units='0/1', flushval = flushone, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_istrimmable_co ) + call this%set_restart_var(vname='fates_gsblaweight',vtype=cohort_r8, & long_name='ed cohort - leaf-area weighted total stomatal+blayer conductance', & units='[m/s]*[m2]', flushval = flushzero, & @@ -1655,6 +1659,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_pft_co => this%rvars(ir_pft_co)%int1d, & rio_status_co => this%rvars(ir_status_co)%int1d, & rio_isnew_co => this%rvars(ir_isnew_co)%int1d, & + rio_istrimmable_co => this%rvars(ir_istrimmable_co)%int1d, & rio_gnd_alb_dif_pasb => this%rvars(ir_gnd_alb_dif_pasb)%r81d, & rio_gnd_alb_dir_pasb => this%rvars(ir_gnd_alb_dir_pasb)%r81d, & rio_spread_si => this%rvars(ir_spread_si)%r81d, & @@ -1893,9 +1898,14 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_pft_co(io_idx_co) = ccohort%pft rio_status_co(io_idx_co) = ccohort%status_coh if ( ccohort%isnew ) then - rio_isnew_co(io_idx_co) = new_cohort + rio_isnew_co(io_idx_co) = itrue else - rio_isnew_co(io_idx_co) = old_cohort + rio_isnew_co(io_idx_co) = ifalse + endif + if ( ccohort%is_trimmable ) then + rio_istrimmable_co(io_idx_co) = itrue + else + rio_istrimmable_co(io_idx_co) = ifalse endif if ( debug ) then @@ -2439,6 +2449,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_pft_co => this%rvars(ir_pft_co)%int1d, & rio_status_co => this%rvars(ir_status_co)%int1d, & rio_isnew_co => this%rvars(ir_isnew_co)%int1d, & + rio_istrimmable_co => this%rvars(ir_istrimmable_co)%int1d, & rio_gnd_alb_dif_pasb => this%rvars(ir_gnd_alb_dif_pasb)%r81d, & rio_gnd_alb_dir_pasb => this%rvars(ir_gnd_alb_dir_pasb)%r81d, & rio_spread_si => this%rvars(ir_spread_si)%r81d, & @@ -2645,8 +2656,9 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ccohort%resp_tstep = rio_resp_tstep_co(io_idx_co) ccohort%pft = rio_pft_co(io_idx_co) ccohort%status_coh = rio_status_co(io_idx_co) - ccohort%isnew = ( rio_isnew_co(io_idx_co) .eq. new_cohort ) - + ccohort%isnew = ( rio_isnew_co(io_idx_co) .eq. itrue ) + ccohort%is_trimmable = ( rio_istrimmable_co(io_idx_co) .eq. itrue ) + call UpdateCohortBioPhysRates(ccohort) From 6a978156b31cb103863ab0f8234a68e28e402e61 Mon Sep 17 00:00:00 2001 From: Gregory Lemieux Date: Mon, 12 Sep 2022 11:43:21 -0700 Subject: [PATCH 02/12] add year_net_uptake to restart --- main/FatesRestartInterfaceMod.F90 | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 7de6fcb5aa..be75c6b424 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -133,6 +133,7 @@ module FatesRestartInterfaceMod integer :: ir_daily_p_demand_co integer :: ir_daily_n_need_co integer :: ir_daily_p_need_co + integer :: ir_year_net_up_co !Logging integer :: ir_lmort_direct_co @@ -839,6 +840,11 @@ subroutine define_restart_vars(this, initialize_variables) units='kgN/plant/day', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_daily_n_need_co ) + call this%set_restart_var(vname='fates_year_net_up', vtype=cohort_r8, & + long_name='fates cohort- yearly net uptake', & + units='kgC/m2/year', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_year_net_up_co ) + call this%set_restart_var(vname='fates_frmort', vtype=cohort_r8, & long_name='ed cohort - freezing mortality rate', & units='/year', flushval = flushzero, & @@ -1803,6 +1809,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_daily_p_demand_co => this%rvars(ir_daily_p_demand_co)%r81d, & rio_daily_n_need_co => this%rvars(ir_daily_n_need_co)%r81d, & rio_daily_p_need_co => this%rvars(ir_daily_p_need_co)%r81d, & + rio_year_net_up_co => this%rvars(ir_year_net_up_co)%r81d, & rio_smort_co => this%rvars(ir_smort_co)%r81d, & rio_asmort_co => this%rvars(ir_asmort_co)%r81d, & rio_frmort_co => this%rvars(ir_frmort_co)%r81d, & @@ -1985,6 +1992,10 @@ subroutine set_restart_vectors(this,nc,nsites,sites) end do end do + do i = 1,nlevleaf + this%rvars(ir_year_net_up_co)%r81d(io_idx_co) = ccohort%year_net_uptake(i) + end do + if(hlm_use_planthydro==itrue)then @@ -2643,6 +2654,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_daily_p_demand_co => this%rvars(ir_daily_p_demand_co)%r81d, & rio_daily_n_need_co => this%rvars(ir_daily_n_need_co)%r81d, & rio_daily_p_need_co => this%rvars(ir_daily_p_need_co)%r81d, & + rio_year_net_up_co => this%rvars(ir_year_net_up_co)%r81d, & rio_smort_co => this%rvars(ir_smort_co)%r81d, & rio_asmort_co => this%rvars(ir_asmort_co)%r81d, & rio_frmort_co => this%rvars(ir_frmort_co)%r81d, & @@ -2893,6 +2905,10 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ccohort%treesai = this%rvars(ir_treesai_co)%r81d(io_idx_co) end if + do i = 1,nlevleaf + ccohort%year_net_uptake(i) = this%rvars(ir_year_net_up_co)%r81d(io_idx_co) + end do + io_idx_co = io_idx_co + 1 ccohort => ccohort%taller From 2bdb0df45236a3cbaee6367550e9dac2f250a694 Mon Sep 17 00:00:00 2001 From: Gregory Lemieux Date: Mon, 12 Sep 2022 15:31:29 -0600 Subject: [PATCH 03/12] add nlevleaf to restart --- main/FatesRestartInterfaceMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index be75c6b424..1b296ea155 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -36,7 +36,7 @@ module FatesRestartInterfaceMod use FatesLitterMod, only : litter_type use FatesLitterMod, only : ncwd use FatesLitterMod, only : ndcmpy - use EDTypesMod, only : nfsc + use EDTypesMod, only : nfsc, nlevleaf use PRTGenericMod, only : prt_global use PRTGenericMod, only : num_elements use FatesRunningMeanMod, only : rmean_type From 58da0ee42e7938244008ae07eb458be5581905fa Mon Sep 17 00:00:00 2001 From: Gregory Lemieux Date: Tue, 13 Sep 2022 17:06:37 -0700 Subject: [PATCH 04/12] add leaf layer index for each cohort --- main/FatesRestartInterfaceMod.F90 | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 1b296ea155..f66ad00598 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -1733,6 +1733,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) integer :: io_idx_si_pft ! each site-pft index integer :: io_idx_si_vtmem ! indices for veg-temp memory at site integer :: io_idx_pa_ncl ! each canopy layer within each patch + integer :: io_idx_co_nll ! each leaf layer within each cohort ! Some counters (for checking mostly) integer :: totalcohorts ! total cohort count on this thread (diagnostic) @@ -1877,6 +1878,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) io_idx_si_wmem = io_idx_co_1st io_idx_si_vtmem = io_idx_co_1st io_idx_pa_ncl = io_idx_co_1st + io_idx_co_nll = io_idx_co_1st ! Hydraulics counters lyr = hydraulic layer, shell = rhizosphere shell io_idx_si_lyr_shell = io_idx_co_1st @@ -1993,7 +1995,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) end do do i = 1,nlevleaf - this%rvars(ir_year_net_up_co)%r81d(io_idx_co) = ccohort%year_net_uptake(i) + this%rvars(ir_year_net_up_co)%r81d(io_idx_co_nll) = ccohort%year_net_uptake(i) + io_idx_co_nll = io_idx_co_nll + 1 end do @@ -2586,6 +2589,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) integer :: io_idx_si_cwd integer :: io_idx_si_pft integer :: io_idx_pa_ncl ! each canopy layer within each patch + integer :: io_idx_co_nll ! each leaf layer within each cohort ! Some counters (for checking mostly) integer :: totalcohorts ! total cohort count on this thread (diagnostic) @@ -2711,6 +2715,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) io_idx_si_wmem = io_idx_co_1st io_idx_si_vtmem = io_idx_co_1st io_idx_pa_ncl = io_idx_co_1st + io_idx_co_nll = io_idx_co_1st ! Hydraulics counters lyr = hydraulic layer, shell = rhizosphere shell io_idx_si_lyr_shell = io_idx_co_1st @@ -2906,7 +2911,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) end if do i = 1,nlevleaf - ccohort%year_net_uptake(i) = this%rvars(ir_year_net_up_co)%r81d(io_idx_co) + ccohort%year_net_uptake(i) = this%rvars(ir_year_net_up_co)%r81d(io_idx_co_nll) + io_idx_co_nll = io_idx_co_nll + 1 end do io_idx_co = io_idx_co + 1 From 25efd59bb87a5cfb538e3b8261b1e63bb96e5fbf Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 22 Nov 2022 23:36:37 -0500 Subject: [PATCH 05/12] updates to leaf memory arrays --- biogeochem/EDCanopyStructureMod.F90 | 14 ++++---- biogeochem/EDCohortDynamicsMod.F90 | 28 +++++++++++----- biogeochem/EDPhysiologyMod.F90 | 14 +++++--- biogeochem/FatesAllometryMod.F90 | 33 +++++++++---------- biogeophys/EDAccumulateFluxesMod.F90 | 9 +++-- biogeophys/FatesPlantRespPhotosynthMod.F90 | 22 ++++++++++--- main/EDInitMod.F90 | 38 +++++++++++++++++++++- main/EDMainMod.F90 | 4 +-- main/FatesConstantsMod.F90 | 10 ++++++ main/FatesHistoryInterfaceMod.F90 | 6 ++-- 10 files changed, 126 insertions(+), 52 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index ee6f3052ee..c47d3ceb5d 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1623,7 +1623,7 @@ subroutine leaf_area_profile( currentSite ) if(remainder > dinc_vai(iv) )then write(fates_log(), *)'ED: issue with remainder', & currentCohort%treelai,currentCohort%treesai,dinc_vai(iv), & - currentCohort%nveg_act,remainder + dlower_vai(iv),currentCohort%nveg_act,remainder call endrun(msg=errMsg(sourcefile, __LINE__)) endif @@ -2185,14 +2185,12 @@ subroutine UpdateCohortLAI(currentCohort, canopy_layer_tlai, total_canopy_area) canopy_layer_tlai, currentCohort%treelai , & currentCohort%vcmax25top,4) end if - - ! Number of actual vegetation layers in this cohort's crown - call GetNLevVeg(currentCohort%dbh, leaf_c, currentSite%spread, currentCohort%pft, & - currentCohort%canopy_trim, currentCohort%canopy_layer, & - currentPatch%canopy_layer_tlai, currentCohort%vcmax25top, & - currentCohort%nveg_act,currentCohort%nveg_max) - + ! Number of actual vegetation layers in this cohort's crown + call GetNLevVeg(currentCohort%dbh, leaf_c, currentCohort%c_area, currentCohort%pft, currentCohort%n, & + currentCohort%crowndamage, currentCohort%canopy_trim, currentCohort%canopy_layer, & + canopy_layer_tlai, currentCohort%vcmax25top, & + currentCohort%nveg_act, currentCohort%nveg_max) end subroutine UpdateCohortLAI diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index e20e1003e9..a28d88d17a 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -1218,7 +1218,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) write(fates_log(),*) 'canopy_layer_yesterday:', & currentCohort%canopy_layer_yesterday,nextc%canopy_layer_yesterday write(fates_log(),*) 'is_trimmable:',currentCohort%is_trimmable,nextc%is_trimmable - do i=1, nlevleaf + do i=1, nlevleafmem write(fates_log(),*) 'leaf level: ',i,'year_net_uptake', & currentCohort%year_net_uptake(i),nextc%year_net_uptake(i) end do @@ -1248,7 +1248,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) ! ----------------------------------------------------------------- call UpdateCohortBioPhysRates(currentCohort) - call FuseVegLayerMem(currentCohort,nextc,currentPatch,currentSite%spread) + ! c13disc_acc calculation; weighted mean by GPP if ((currentCohort%n * currentCohort%gpp_acc + nextc%n * nextc%gpp_acc) .eq. 0.0_r8) then @@ -1491,6 +1491,11 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) currentCohort%ddbhdt = (currentCohort%n*currentCohort%ddbhdt + & nextc%n*nextc%ddbhdt)/newn + ! Try to do this last, it relies on many things like actual leaf carbon + ! and crown area. We have to call the tree_lai and tree_sai routines + ! to get the updated layering. + call FuseVegLayerMem(currentCohort,nextc,currentPatch,currentSite%spread) + end if !(currentCohort%isnew) currentCohort%n = newn @@ -1855,9 +1860,10 @@ subroutine copy_cohort( currentCohort,copyc ) n%resp_tstep = o%resp_tstep n%resp_acc = o%resp_acc n%resp_acc_hold = o%resp_acc_hold - n%year_net_uptake = o%year_net_uptake - n%ts_net_uptake = o%ts_net_uptake + n%year_net_uptake = o%year_net_uptake + n%is_trimmable = o%is_trimmable + n%daily_nh4_uptake = o%daily_nh4_uptake n%daily_no3_uptake = o%daily_no3_uptake n%daily_p_uptake = o%daily_p_uptake @@ -2000,24 +2006,28 @@ subroutine FuseVegLayerMem(ccohort,ncohort,cpatch,site_spread) do iv = iv0 ,iv1 if(joint_wgt(iv) nlevleaf)then @@ -857,7 +859,9 @@ subroutine trim_canopy( currentSite ) ! that are not zero'd here, are not allowed to trim (because ! they need full year) currentCohort%year_net_uptake(:) = 0.0_r8 - + currentCohort%is_trimmable = .true. + + ! Add to trim fraction if cohort not trimmed at all if ( (.not.trimmed) .and.currentCohort%canopy_trim < 1.0_r8)then currentCohort%canopy_trim = currentCohort%canopy_trim + EDPftvarcon_inst%trim_inc(ipft) diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 79b37c4271..48d930d849 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -802,9 +802,9 @@ real(r8) function tree_sai(pft, dbh, crowndamage, canopy_trim, c_area, nplant, c end function tree_sai - subroutine GetNLevVeg(dbh, leaf_c, site_spread, ipft, canopy_trim, & - canopy_layer,canopy_layer_tlai, vcmax25top, & - nveg_act, nveg_max) + subroutine GetNLevVeg(dbh, leaf_c, c_area, ipft, nplant, & + crowndamage, canopy_trim, canopy_layer, canopy_layer_tlai, & + vcmax25top, nveg_act, nveg_max) ! Determine the maximum number of leaf (vegetation) layers ! for a cohort. This is based off of allometry, and assuming the plant @@ -812,9 +812,11 @@ subroutine GetNLevVeg(dbh, leaf_c, site_spread, ipft, canopy_trim, & real(r8), intent(inout) :: dbh real(r8), intent(in) :: leaf_c - real(r8), intent(in) :: site_spread + real(r8), intent(in) :: c_area integer, intent(in) :: ipft + real(r8), intent(in) :: nplant real(r8), intent(in) :: canopy_trim + integer, intent(in) :: crowndamage integer, intent(in) :: canopy_layer real(r8), intent(in) :: canopy_layer_tlai(:) real(r8), intent(in) :: vcmax25top @@ -822,32 +824,29 @@ subroutine GetNLevVeg(dbh, leaf_c, site_spread, ipft, canopy_trim, & integer, intent(out) :: nveg_max integer :: nv - real(r8) :: c_area real(r8) :: leaf_c_target real(r8) :: max_lai, max_sai - real(r8), parameter :: nplant_dum = 1.0_r8 ! Number of plants falls out (dummy) - call carea_allom(dbh,nplant_dum,site_spread,ipft,c_area) max_lai = tree_lai(leaf_c, ipft, c_area, & - nplant_dum, canopy_layer, canopy_layer_tlai, vcmax25top ) + nplant, canopy_layer, canopy_layer_tlai, vcmax25top ) - max_sai = tree_sai(ipft, dbh, canopy_trim, & - c_area, nplant_dum, canopy_layer, canopy_layer_tlai, max_lai, & - vcmax25top, 0) + max_sai = tree_sai(ipft, dbh, crowndamage, canopy_trim, & + c_area, nplant, canopy_layer, canopy_layer_tlai, max_lai, & + vcmax25top, 0) - nveg_act = count((currentCohort%treelai+currentCohort%treesai) .gt. dlower_vai(:)) + 1 + nveg_act = count((max_lai+max_sai) .gt. dlower_vai(:)) + 1 - call bleaf(dbh,ipft,canopy_trim,leaf_c_target) + call bleaf(dbh,ipft,crowndamage,canopy_trim,leaf_c_target) max_lai = tree_lai(leaf_c_target, ipft, c_area, & - nplant_dum, canopy_layer, canopy_layer_tlai, vcmax25top ) + nplant, canopy_layer, canopy_layer_tlai, vcmax25top ) - max_sai = tree_sai(ipft, dbh, canopy_trim, & - c_area, nplant_dum, canopy_layer, canopy_layer_tlai, max_lai, & + max_sai = tree_sai(ipft, dbh, crowndamage, canopy_trim, & + c_area, nplant, canopy_layer, canopy_layer_tlai, max_lai, & vcmax25top, 0) - nveg_max = count((currentCohort%treelai+currentCohort%treesai) .gt. dlower_vai(:)) + 1 + nveg_max = count((max_lai+max_sai) .gt. dlower_vai(:)) + 1 return end subroutine GetNLevVeg diff --git a/biogeophys/EDAccumulateFluxesMod.F90 b/biogeophys/EDAccumulateFluxesMod.F90 index 50796b5d51..3aa6bbf7e8 100644 --- a/biogeophys/EDAccumulateFluxesMod.F90 +++ b/biogeophys/EDAccumulateFluxesMod.F90 @@ -57,6 +57,7 @@ subroutine AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) integer :: c ! clm/alm column integer :: s ! ed site integer :: ifp ! index fates patch + integer :: ifc !---------------------------------------------------------------------- do s = 1, nsites @@ -67,11 +68,15 @@ subroutine AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) do while (associated(cpatch)) if(cpatch%nocomp_pft_label.ne.0)then ifp = ifp+1 - + if( bc_in(s)%filter_photo_pa(ifp) == 3 ) then + + ifc=0 ccohort => cpatch%shortest do while(associated(ccohort)) + ifc=ifc+1 + ! Accumulate fluxes from hourly to daily values. ! _tstep fluxes are KgC/indiv/timestep _acc are KgC/indiv/day @@ -96,7 +101,7 @@ subroutine AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) if( ccohort%is_trimmable ) then - do iv=1,min(ccohort%nveg_max,nlevleafmem) + do iv=1,min(ccohort%nveg_act,nlevleafmem) ccohort%year_net_uptake(iv) = ccohort%year_net_uptake(iv) + ccohort%ts_net_uptake(iv) end do end if diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 012ec2dee6..d8d3f218ef 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -261,6 +261,8 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) real(r8) :: sapw_n_undamaged ! nitrogen in sapwood of undamaged tree integer :: ivm, iv, nvm, nva ! Loop indices and bounds for cohort vegetation layers + integer :: ivt, ivb ! Top and bottom leaf layer indices that match + ! viable leaf memory bins integer :: cl, s, j, ps, ft, ifp ! indices integer :: NCL_p ! number of canopy layers in patch integer :: iage ! loop counter for leaf age classes @@ -652,13 +654,23 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) if(nva>0) then if(.not.stretch_cbalmem_vert)then - do ivm = min(nlevleafmem,nvm-nva+1),nlevleafmem - iv = nva - ivm + 1 - currentCohort%ts_net_uptake(ivm) = anet_av_z(iv,ft,cl) * umolC_to_kgC - end do + + ivb = nva + ivt = max(1,nvm-nlevleafmem+1) + if(ivb>=ivt) then + do iv = ivt,ivb + ivm = nvm - iv + 1 + ! Its possible that the leaves aren't flushed enough + ! to lign up with any of the memory bins + if(ivm>0 .and. ivm.le.nlevleafmem) then + currentCohort%ts_net_uptake(ivm) = anet_av_z(iv,ft,cl) * umolC_to_kgC + end if + end do + end if + else - do ivm = 1,nlevleafmem + do ivm = 1,min(nlevleafmem,nva) ! This is the fraction of actual leaf bins versus the ! current allometric maximum diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index b8697ed6c7..65c2accb1f 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -22,7 +22,8 @@ module EDInitMod use EDPatchDynamicsMod , only : set_patchno use EDPhysiologyMod , only : assign_cohort_sp_properties use ChecksBalancesMod , only : SiteMassStock - use FatesInterfaceTypesMod , only : hlm_day_of_year + use FatesInterfaceTypesMod , only : hlm_day_of_year, hlm_days_per_year + use FatesConstantsMod , only : trim_day_of_year, max_days_since_trim_coldstart use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use EDTypesMod , only : numWaterMem use EDTypesMod , only : num_vegtemp_mem @@ -503,6 +504,7 @@ subroutine init_patches( nsites, sites, bc_in) integer :: start_patch integer :: num_new_patches integer :: nocomp_pft + integer :: days_since_trim ! Number of days since trimming date real(r8) :: newparea real(r8) :: tota !check on area integer :: is_first_patch @@ -510,6 +512,7 @@ subroutine init_patches( nsites, sites, bc_in) type(ed_site_type), pointer :: sitep type(ed_patch_type), pointer :: newppft(:) type(ed_patch_type), pointer :: newp + type(ed_cohort_type),pointer :: newc type(ed_patch_type), pointer :: currentPatch ! List out some nominal patch values that are used for Near Bear Ground initializations @@ -637,6 +640,39 @@ subroutine init_patches( nsites, sites, bc_in) end if end do !no new patches + ! ---------------------------------------------------------- + ! If this is a cold start and the date is after the date + ! we trim on, but fairly close, then we say its close + ! enough and flag the newly created cohorts as trimmable + ! the next time we encounter trimming, and zero out + ! the incrementing carbon balance array + ! ----------------------------------------------------------- + + if(hlm_day_of_year.ge.trim_day_of_year) then + days_since_trim = hlm_day_of_year-trim_day_of_year + else + days_since_trim = (hlm_day_of_year+hlm_days_per_year)-trim_day_of_year + end if + + if(days_since_trim .le. max_days_since_trim_coldstart)then + + write(fates_log(),*) 'FATES Cold-start initiated on a day-of-year ',hlm_day_of_year + write(fates_log(),*) 'after and reasonably close to the trim date ',trim_day_of_year + write(fates_log(),*) 'Will initialize trimming arrays in new cohorts.' + + newp => sites(s)%oldest_patch + do while (associated(newp)) + newc => newp%tallest + do while(associated(newc)) + newc%is_trimmable = .true. + newc%year_net_uptake(:) = 0._r8 + newc => newc%shorter + end do + newp=>newp%younger + end do + end if + + !check if the total area adds to the same as site area tota = 0.0_r8 newp => sites(s)%oldest_patch diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index c99b89bcc3..47a8d09a68 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -75,6 +75,7 @@ module EDMainMod use FatesConstantsMod , only : nearzero use FatesConstantsMod , only : m2_per_ha use FatesConstantsMod , only : sec_per_day + use FatesConstantsMod , only : trim_day_of_year use FatesPlantHydraulicsMod , only : do_growthrecruiteffects use FatesPlantHydraulicsMod , only : UpdateSizeDepPlantHydProps use FatesPlantHydraulicsMod , only : UpdateSizeDepPlantHydStates @@ -831,10 +832,9 @@ subroutine ed_update_site( currentSite, bc_in, bc_out ) ! rooting mass, distributions, respiration rates and NPP call PrepCH4BCs(currentSite,bc_in,bc_out) - ! FIX(RF,032414). This needs to be monthly, not annual ! If this is the second to last day of the year, then perform trimming - if( hlm_day_of_year == hlm_days_per_year-1) then + if( hlm_day_of_year == trim_day_of_year) then if(hlm_use_sp.eq.ifalse)then call trim_canopy(currentSite) endif diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index 0ec836fe0f..e737758982 100644 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -60,7 +60,17 @@ module FatesConstantsMod ! plants are not coupled with below ground chemistry. In ! this situation, we send token boundary condition information. + + integer, public, parameter :: trim_day_of_year = 364 ! Trim on this date (2nd to last day of year) + + ! If this is a cold start and the date is after the date + ! we trim on, but fairly close, then we say its close + ! enough and flag the newly created cohorts as trimmable + ! the next time we encounter trimming, and zero out + ! the incrementing carbon balance array + integer, public, parameter :: max_days_since_trim_coldstart = 5 + ! This flag specifies the scaling of how we present ! nutrient competitors to the HLM's soil BGC model diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 9fbd45cd0c..4af64b19b9 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -3872,7 +3872,7 @@ subroutine update_history_hifrq(this,nc,nsites,sites,bc_in,dt_tstep) ! after rapid timescale productivity calculations (gpp and respiration). ! --------------------------------------------------------------------------------- - use EDTypesMod , only : nclmax, nlevleaf + use EDTypesMod , only : nclmax, nlevleaf, nlevleafmem ! ! Arguments class(fates_history_interface_type) :: this @@ -4167,10 +4167,10 @@ subroutine update_history_hifrq(this,nc,nsites,sites,bc_in,dt_tstep) !!! canopy leaf carbon balance ican = ccohort%canopy_layer - do ileaf=1, min(nlevleafmem,currentCohort%nveg_max) + do ileaf=1, min(nlevleafmem,ccohort%nveg_max) ! The actual leaf layer index (going from top to bottom) - z = currentCohort%nveg_max - ileaf + 1 + z = ccohort%nveg_max - ileaf + 1 cnlf_indx = z + (ican-1) * nlevleaf hio_ts_net_uptake_si_cnlf(io_si, cnlf_indx) = hio_ts_net_uptake_si_cnlf(io_si, cnlf_indx) + & From adec913281d64ee191faa4e9c70f2d5fa84ac439 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 29 Nov 2022 20:51:00 -0500 Subject: [PATCH 06/12] Reversed leaf layer loop during trimming --- biogeochem/EDPhysiologyMod.F90 | 36 +++++++++++++++++++--------------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 9911d4c9c6..af0c25683c 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -704,7 +704,9 @@ subroutine trim_canopy( currentSite ) ! Lets loop over the last few leaf layers, only the number of layers ! necessary for evaluating the trend - leafmem_loop: do zm = 1, min(nll,currentCohort%nveg_max) + cumulative_lai_cohort = 0._r8 + + leafmem_loop: do zm = min(nll,currentCohort%nveg_max),1,-1 ! The actual leaf layer index (going from top to bottom) z = currentCohort%nveg_max - zm + 1 @@ -755,7 +757,6 @@ subroutine trim_canopy( currentSite ) bfr_per_bleaf / prt_params%root_long(ipft) endif - currentCohort%leaf_cost = currentCohort%leaf_cost * & (prt_params%grperc(ipft) + 1._r8) else !evergreen costs @@ -764,7 +765,6 @@ subroutine trim_canopy( currentSite ) currentCohort%leaf_cost = 1.0_r8/(sla_levleaf* & sum(prt_params%leaf_long(ipft,:))*1000.0_r8) !convert from sla in m2g-1 to m2kg-1 - if ( int(prt_params%allom_fmode(ipft)) .eq. 1 ) then ! if using trimmed leaf for fine root biomass allometry, add the cost of the root increment ! to the leaf increment; otherwise do not. @@ -780,19 +780,23 @@ subroutine trim_canopy( currentSite ) ! if at least nll leaf layers are present in the current cohort and only for the bottom nll ! leaf layers. - ! Build the A matrix for the LHS of the linear system. A = [n sum(x); sum(x) sum(x^2)] - ! where n = nll and x = yearly_net_uptake-leafcost - nnu_clai_a(1,1) = nnu_clai_a(1,1) + 1 ! Increment for each layer used - nnu_clai_a(1,2) = nnu_clai_a(1,2) + currentCohort%year_net_uptake(zm) - currentCohort%leaf_cost - nnu_clai_a(2,1) = nnu_clai_a(1,2) - nnu_clai_a(2,2) = nnu_clai_a(2,2) + (currentCohort%year_net_uptake(zm) - currentCohort%leaf_cost)**2 - - ! Build the B matrix for the RHS of the linear system. B = [sum(y); sum(x*y)] - ! where x = yearly_net_uptake-leafcost and y = cumulative_lai_cohort - nnu_clai_b(1,1) = nnu_clai_b(1,1) + cumulative_lai_cohort - nnu_clai_b(2,1) = nnu_clai_b(2,1) + (cumulative_lai_cohort * & - (currentCohort%year_net_uptake(zm) - currentCohort%leaf_cost)) - + if( currentCohort%nveg_max >= nll ) then + + ! Build the A matrix for the LHS of the linear system. A = [n sum(x); sum(x) sum(x^2)] + ! where n = nll and x = yearly_net_uptake-leafcost + nnu_clai_a(1,1) = nnu_clai_a(1,1) + 1 ! Increment for each layer used + nnu_clai_a(1,2) = nnu_clai_a(1,2) + currentCohort%year_net_uptake(zm) - currentCohort%leaf_cost + nnu_clai_a(2,1) = nnu_clai_a(1,2) + nnu_clai_a(2,2) = nnu_clai_a(2,2) + (currentCohort%year_net_uptake(zm) - currentCohort%leaf_cost)**2 + + ! Build the B matrix for the RHS of the linear system. B = [sum(y); sum(x*y)] + ! where x = yearly_net_uptake-leafcost and y = cumulative_lai_cohort + nnu_clai_b(1,1) = nnu_clai_b(1,1) + cumulative_lai_cohort + nnu_clai_b(2,1) = nnu_clai_b(2,1) + (cumulative_lai_cohort * & + (currentCohort%year_net_uptake(zm) - currentCohort%leaf_cost)) + + end if + ! Check leaf cost against the yearly net uptake for that cohort leaf layer if (currentCohort%year_net_uptake(zm) < currentCohort%leaf_cost) then ! Make sure the cohort trim fraction is great than the pft trim limit From 19c9fbb25346a963c4b4fbb3e6670200ab0b2eba Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Sun, 18 Dec 2022 22:42:05 -0500 Subject: [PATCH 07/12] reshaping leaflayermem --- biogeochem/DamageMainMod.F90 | 15 --- biogeochem/EDCanopyStructureMod.F90 | 51 ++++----- biogeochem/EDCohortDynamicsMod.F90 | 119 +++++++++++++++---- biogeochem/EDPhysiologyMod.F90 | 32 +++--- biogeochem/FatesAllometryMod.F90 | 116 +++++++------------ biogeophys/EDAccumulateFluxesMod.F90 | 2 +- biogeophys/FatesPlantRespPhotosynthMod.F90 | 59 ++++++---- main/EDInitMod.F90 | 3 +- main/EDMainMod.F90 | 2 +- main/EDPftvarcon.F90 | 3 +- main/EDTypesMod.F90 | 126 ++++++++++++++------- main/FatesConstantsMod.F90 | 10 -- main/FatesRestartInterfaceMod.F90 | 4 - 13 files changed, 299 insertions(+), 243 deletions(-) diff --git a/biogeochem/DamageMainMod.F90 b/biogeochem/DamageMainMod.F90 index f0a05f7ee6..69cc6dbfae 100644 --- a/biogeochem/DamageMainMod.F90 +++ b/biogeochem/DamageMainMod.F90 @@ -29,13 +29,11 @@ module DamageMainMod character(len=*), parameter, private :: sourcefile = & __FILE__ - public :: GetCrownReduction public :: GetDamageFrac public :: IsItDamageTime public :: damage_time public :: GetDamageMortality - logical :: debug = .false. ! for debugging @@ -163,20 +161,7 @@ end subroutine GetDamageFrac !------------------------------------------------------- - subroutine GetCrownReduction(crowndamage, crown_reduction) - - !------------------------------------------------------------------ - ! This subroutine takes the crown damage class of a cohort (integer) - ! and returns the fraction of the crown that is lost. - !------------------------------------------------------------------- - integer(i4), intent(in) :: crowndamage ! crown damage class of the cohort - real(r8), intent(out) :: crown_reduction ! fraction of crown lost from damage - - crown_reduction = ED_val_history_damage_bin_edges(crowndamage)/100.0_r8 - - return - end subroutine GetCrownReduction !---------------------------------------------------------------------------------------- diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 0eb1166334..86244562ec 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -19,11 +19,12 @@ module EDCanopyStructureMod use EDCohortDynamicsMod , only : InitPRTBoundaryConditions use FatesAllometryMod , only : tree_lai use FatesAllometryMod , only : tree_sai - use FatesAllometryMod , only : GetNLevVeg + use EDCohortDynamicsMod , only : GetNLevVeg use EDtypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use EDTypesMod , only : nclmax use EDTypesMod , only : nlevleaf - use EDtypesMod , only : AREA + use EDtypesMod , only : area, dinc_vai, dlower_vai + use EDTypesMod , only : n_hite_bins use EDLoggingMortalityMod , only : UpdateHarvestC use FatesGlobals , only : endrun => fates_endrun use FatesInterfaceTypesMod , only : hlm_days_per_year @@ -31,7 +32,7 @@ module EDCanopyStructureMod use FatesInterfaceTypesMod , only : hlm_use_cohort_age_tracking use FatesInterfaceTypesMod , only : hlm_use_sp use FatesInterfaceTypesMod , only : numpft - use FatesInterfaceTypesMod, only : bc_in_type + use FatesInterfaceTypesMod, only : bc_in_type,bc_out_type use FatesPlantHydraulicsMod, only : UpdateH2OVeg,InitHydrCohort, RecruitWaterStorage use PRTGenericMod, only : leaf_organ use PRTGenericMod, only : leaf_organ @@ -1229,7 +1230,6 @@ subroutine canopy_spread( currentSite ) ! Calculates the spatial spread of tree canopies based on canopy closure. ! ! !USES: - use EDTypesMod , only : AREA use EDParamsMod, only : ED_val_canopy_closure_thresh ! ! !ARGUMENTS @@ -1268,7 +1268,7 @@ subroutine canopy_spread( currentSite ) ! If the canopy area is approaching closure, ! squash the tree canopies and make them taller and thinner - if( sitelevel_canopyarea/AREA .gt. ED_val_canopy_closure_thresh ) then + if( sitelevel_canopyarea/area .gt. ED_val_canopy_closure_thresh ) then currentSite%spread = currentSite%spread - inc else currentSite%spread = currentSite%spread + inc @@ -1494,7 +1494,7 @@ subroutine leaf_area_profile( currentSite ) ! !USES: - use EDtypesMod , only : area, dinc_vai, dlower_vai, hitemax, n_hite_bins + ! ! !ARGUMENTS @@ -1794,11 +1794,6 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out) ! to vegetation coverage to the host land model. ! ---------------------------------------------------------------------------------- - use EDTypesMod , only : ed_patch_type, ed_cohort_type, & - ed_site_type, AREA - use FatesInterfaceTypesMod , only : bc_out_type - - ! ! !ARGUMENTS integer, intent(in) :: nsites type(ed_site_type), intent(inout), target :: sites(nsites) @@ -1904,17 +1899,17 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out) ! The alternative is to assume it is all spatial distinct from tree canopies. ! In which case, the bare area would have to be reduced by the grass area... ! currentPatch%total_canopy_area/currentPatch%area is fraction of this patch cover by plants - ! currentPatch%area/AREA is the fraction of the soil covered by this patch. + ! currentPatch%area/area is the fraction of the soil covered by this patch. if(currentPatch%area.gt.0.0_r8)then bc_out(s)%canopy_fraction_pa(ifp) = & - min(1.0_r8,currentPatch%total_canopy_area/currentPatch%area)*(currentPatch%area/AREA) + min(1.0_r8,currentPatch%total_canopy_area/currentPatch%area)*(currentPatch%area/area) else bc_out(s)%canopy_fraction_pa(ifp) = 0.0_r8 endif bare_frac_area = (1.0_r8 - min(1.0_r8,currentPatch%total_canopy_area/currentPatch%area)) * & - (currentPatch%area/AREA) + (currentPatch%area/area) total_patch_area = total_patch_area + bc_out(s)%canopy_fraction_pa(ifp) + bare_frac_area @@ -1947,7 +1942,7 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out) else ! nocomp or SP, and currentPatch%nocomp_pft_label .eq. 0 - total_patch_area = total_patch_area + currentPatch%area/AREA + total_patch_area = total_patch_area + currentPatch%area/area end if currentPatch => currentPatch%younger @@ -2120,7 +2115,6 @@ subroutine UpdatePatchLAI(currentPatch) ! --------------------------------------------------------------------------------------------- ! Uses - use EDtypesMod, only : dlower_vai ! Arguments type(ed_patch_type),intent(inout), target :: currentPatch @@ -2147,8 +2141,7 @@ subroutine UpdatePatchLAI(currentPatch) ft = currentCohort%pft ! Update the cohort level lai and related variables - call UpdateCohortLAI(currentCohort,currentPatch%canopy_layer_tlai, & - currentPatch%total_canopy_area) + call UpdateCohortLAI(currentCohort,currentPatch) currentPatch%ncan(cl,ft) = max(currentPatch%ncan(cl,ft),currentCohort%nveg_act) @@ -2165,17 +2158,14 @@ subroutine UpdatePatchLAI(currentPatch) end subroutine UpdatePatchLAI ! =============================================================================================== - subroutine UpdateCohortLAI(currentCohort, canopy_layer_tlai, total_canopy_area) + subroutine UpdateCohortLAI(currentCohort, currentPatch) ! Update LAI and related variables for a given cohort ! Uses - use EDtypesMod, only : dlower_vai - ! Arguments - type(ed_cohort_type),intent(inout), target :: currentCohort - real(r8), intent(in) :: canopy_layer_tlai(nclmax) ! total leaf area index of each canopy layer - real(r8), intent(in) :: total_canopy_area ! either patch%total_canopy_area or patch%area + type(ed_cohort_type) :: currentCohort + type(ed_patch_type) :: currentPatch ! Local variables real(r8) :: leaf_c ! leaf carbon [kg] @@ -2187,23 +2177,20 @@ subroutine UpdateCohortLAI(currentCohort, canopy_layer_tlai, total_canopy_area) ! Note that tree_lai has an internal check on the canopy locatoin currentCohort%treelai = tree_lai(leaf_c, currentCohort%pft, currentCohort%c_area, & currentCohort%n, currentCohort%canopy_layer, & - canopy_layer_tlai,currentCohort%vcmax25top ) + currentPatch%canopy_layer_tlai, currentCohort%vcmax25top ) if (hlm_use_sp .eq. ifalse) then currentCohort%treesai = tree_sai(currentCohort%pft, currentCohort%dbh, currentCohort%crowndamage, & currentCohort%canopy_trim, & currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, & - canopy_layer_tlai, currentCohort%treelai , & - currentCohort%vcmax25top,4) + currentPatch%canopy_layer_tlai, currentCohort%treelai , & + sum(dinc_vai),currentCohort%vcmax25top,4) end if ! Number of actual vegetation layers in this cohort's crown - call GetNLevVeg(currentCohort%dbh, leaf_c, currentCohort%c_area, currentCohort%pft, currentCohort%n, & - currentCohort%crowndamage, currentCohort%canopy_trim, currentCohort%canopy_layer, & - canopy_layer_tlai, currentCohort%vcmax25top, & - currentCohort%nveg_act, currentCohort%nveg_max) + call GetNLevVeg(currentCohort,currentPatch) - + return end subroutine UpdateCohortLAI ! =============================================================================================== diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 3730e2d045..417844c0c6 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -39,11 +39,14 @@ Module EDCohortDynamicsMod use EDTypesMod , only : min_n_safemath use EDTypesMod , only : nlevleaf use EDTypesMod , only : nlevleafmem + use EDTypesMod , only : use_relative_leafmem_layers + use EDTypesMod , only : GetLeafFromMemLayer use PRTGenericMod , only : max_nleafage use EDTypesMod , only : ican_upper use EDTypesMod , only : site_fluxdiags_type - use PRTGenericMod , only : num_elements + use PRTGenericMod , only : num_elements use EDTypesMod , only : leaves_on + use EDTypesMod , only : dinc_vai,dlower_vai use EDParamsMod , only : ED_val_cohort_age_fusion_tol use FatesInterfaceTypesMod , only : hlm_use_planthydro use FatesInterfaceTypesMod , only : hlm_parteh_mode @@ -74,7 +77,6 @@ Module EDCohortDynamicsMod use FatesAllometryMod , only : tree_lai, tree_sai use FatesAllometryMod , only : set_root_fraction use PRTGenericMod, only : prt_carbon_allom_hyp - use FatesAllometryMod , only : GetNLevVeg use PRTGenericMod, only : prt_cnp_flex_allom_hyp use PRTGenericMod, only : prt_vartypes use PRTGenericMod, only : carbon12_element @@ -134,6 +136,7 @@ Module EDCohortDynamicsMod public :: DeallocateCohort public :: EvaluateAndCorrectDBH public :: DamageRecovery + public :: GetNLevVeg logical, parameter :: debug = .false. ! local debug flag @@ -299,7 +302,8 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, coage, dbh, & new_cohort%treesai = tree_sai(new_cohort%pft, new_cohort%dbh, & new_cohort%crowndamage, new_cohort%canopy_trim, & new_cohort%c_area, new_cohort%n, new_cohort%canopy_layer, & - patchptr%canopy_layer_tlai, new_cohort%treelai,new_cohort%vcmax25top,2 ) + patchptr%canopy_layer_tlai, new_cohort%treelai, sum(dinc_vai), & + new_cohort%vcmax25top,2 ) end if @@ -2013,31 +2017,34 @@ subroutine FuseVegLayerMem(ccohort,ncohort,cpatch,site_spread) real(r8) :: joint_net_uptake(nlevleaf) real(r8) :: joint_wgt(nlevleaf) - real(r8) :: leaf_c integer :: iv ! veg layer index integer :: ivm ! veg layer memory index integer :: iv0,iv1 ! lowest and highest veg layer index matching memory joint_net_uptake(:) = 0._r8 joint_wgt(:) = 0._r8 + + ccohort%canopy_trim = (ccohort%n*ccohort%canopy_trim & + + ncohort%n*ncohort%canopy_trim)/(ccohort%n+ncohort%n) + if(ccohort%is_trimmable .and. ncohort%is_trimmable) then - ccohort%canopy_trim = (ccohort%n*ccohort%canopy_trim & - + ncohort%n*ncohort%canopy_trim)/(ccohort%n+ncohort%n) + !ccohort%canopy_trim = (ccohort%n*ccohort%canopy_trim & + ! + ncohort%n*ncohort%canopy_trim)/(ccohort%n+ncohort%n) do ivm = 1, min(nlevleafmem,ccohort%nveg_max) - iv = ccohort%nveg_max - ivm + 1 + iv = GetLeafFromMemLayer(ccohort,ivm) joint_net_uptake(iv) = ccohort%n*ccohort%year_net_uptake(ivm) joint_wgt(iv) = ccohort%n end do do ivm = 1, min(nlevleafmem,ncohort%nveg_max) - iv = ncohort%nveg_max - ivm + 1 + iv = GetLeafFromMemLayer(ncohort,ivm) joint_net_uptake(iv) = joint_net_uptake(iv)+ncohort%n*ncohort%year_net_uptake(ivm) joint_wgt(iv) = joint_wgt(iv)+ncohort%n end do - iv0 = min(max(ccohort%nveg_max-nlevleafmem+1,1),max(ncohort%nveg_max-nlevleafmem+1,1)) + iv0 = max(1,min(ccohort%nveg_max,ncohort%nveg_max)-nlevleafmem+1) iv1 = max(ccohort%nveg_max,ncohort%nveg_max) if(debug)then @@ -2053,6 +2060,8 @@ subroutine FuseVegLayerMem(ccohort,ncohort,cpatch,site_spread) if(joint_wgt(iv)0) then + cohort%year_net_uptake(1+di:nlevleafmem) = cohort%year_net_uptake(1:nlevleafmem-di) + cohort%year_net_uptake(1:di) = 0._r8 + else + cohort%year_net_uptake(1:nlevleafmem+di) = cohort%year_net_uptake(1-di:nlevleafmem) + cohort%year_net_uptake(nlevleafmem+di+1:nlevleafmem) = 0._r8 + end if + end if + end if + + cohort%nveg_max = nveg_max + + return + end subroutine GetNLevVeg + end module EDCohortDynamicsMod diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 7125b5d4ff..781501084f 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -33,6 +33,7 @@ module EDPhysiologyMod use EDCohortDynamicsMod , only : InitPRTObject use EDCohortDynamicsMod , only : InitPRTBoundaryConditions use EDCohortDynamicsMod , only : copy_cohort + use EDCohortDynamicsMod , only : GetNLevVeg use FatesAllometryMod , only : tree_lai use FatesAllometryMod , only : tree_sai use FatesAllometryMod , only : leafc_from_treelai @@ -42,6 +43,7 @@ module EDPhysiologyMod use EDTypesMod , only : numlevsoil_max use EDTypesMod , only : numWaterMem use EDTypesMod , only : dl_sf, dinc_vai, dlower_vai, area_inv + use EDTypesMod , only : GetLeafFromMemLayer use EDTypesMod , only : AREA use FatesLitterMod , only : ncwd use FatesLitterMod , only : ndcmpy @@ -70,6 +72,7 @@ module EDPhysiologyMod use EDTypesMod , only : phen_dstat_moiston use EDTypesMod , only : phen_dstat_timeon use EDTypesMod , only : init_recruit_trim + use EDTypesMod , only : use_relative_leafmem_layers use shr_log_mod , only : errMsg => shr_log_errMsg use FatesGlobals , only : fates_log use FatesGlobals , only : endrun => fates_endrun @@ -93,7 +96,6 @@ module EDPhysiologyMod use FatesAllometryMod , only : carea_allom use FatesAllometryMod , only : CheckIntegratedAllometries use FatesAllometryMod, only : set_root_fraction - use FatesAllometryMod, only : GetNLevVeg use PRTGenericMod, only : prt_carbon_allom_hyp use PRTGenericMod, only : prt_cnp_flex_allom_hyp use PRTGenericMod, only : prt_vartypes @@ -115,7 +117,7 @@ module EDPhysiologyMod use PRTLossFluxesMod, only : PRTDamageLosses use PRTGenericMod, only : StorageNutrientTarget use DamageMainMod, only : damage_time - use DamageMainMod, only : GetCrownReduction + use FatesAllometryMod, only : GetCrownReduction use DamageMainMod, only : GetDamageFrac use SFParamsMod, only : SF_val_CWD_frac use FatesParameterDerivedMod, only : param_derived @@ -664,20 +666,14 @@ subroutine trim_canopy( currentSite ) currentCohort%canopy_trim, & currentCohort%c_area, currentCohort%n,currentCohort%canopy_layer,& currentPatch%canopy_layer_tlai, currentCohort%treelai, & - currentCohort%vcmax25top,0 ) - - call GetNLevVeg(currentCohort%dbh, leaf_c, currentCohort%c_area, & - currentCohort%pft, currentCohort%n, & - currentCohort%crowndamage,currentCohort%canopy_trim, & - currentCohort%canopy_layer, & - currentPatch%canopy_layer_tlai, currentCohort%vcmax25top, & - currentCohort%nveg_act,currentCohort%nveg_max) - + sum(dinc_vai),currentCohort%vcmax25top,0 ) + + call GetNLevVeg(currentCohort,currentPatch) if (currentCohort%nveg_max > nlevleaf)then - write(fates_log(),*) 'nv > nlevleaf',currentCohort%nveg_max, & + write(fates_log(),*) 'nveg_max > nlevleaf',currentCohort%nveg_max, & currentCohort%treelai,currentCohort%treesai, & - currentCohort%c_area,currentCohort%n,leaf_c + currentCohort%c_area,currentCohort%n call endrun(msg=errMsg(sourcefile, __LINE__)) endif @@ -709,10 +705,10 @@ subroutine trim_canopy( currentSite ) cumulative_lai_cohort = 0._r8 - leafmem_loop: do zm = min(nll,currentCohort%nveg_max),1,-1 + leafmem_loop: do zm = min(nlevleafmem,currentCohort%nveg_max),1,-1 ! The actual leaf layer index (going from top to bottom) - z = currentCohort%nveg_max - zm + 1 + z = GetLeafFromMemLayer(currentCohort,zm) ! Calculate the cumulative total vegetation area index (no snow occlusion, stems and leaves) @@ -799,10 +795,11 @@ subroutine trim_canopy( currentSite ) (currentCohort%year_net_uptake(zm) - currentCohort%leaf_cost)) end if - + ! Check leaf cost against the yearly net uptake for that cohort leaf layer if (currentCohort%year_net_uptake(zm) < currentCohort%leaf_cost) then ! Make sure the cohort trim fraction is great than the pft trim limit + if (currentCohort%canopy_trim > EDPftvarcon_inst%trim_limit(ipft)) then ! keep trimming until none of the canopy is in negative carbon balance. @@ -878,6 +875,9 @@ subroutine trim_canopy( currentSite ) write(fates_log(),*) 'trimming:',currentCohort%canopy_trim endif + ! Update the number of vegetation layers, max and actual + call GetNLevVeg(currentCohort,currentPatch) + currentCohort => currentCohort%shorter icohort = icohort + 1 enddo diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 1c791c0241..0907cd721d 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -96,14 +96,10 @@ module FatesAllometryMod use FatesGlobals , only : fates_log use FatesGlobals , only : endrun => fates_endrun use FatesGlobals , only : FatesWarn,N2S,A2S,I2S - use EDTypesMod , only : nlevleaf, dinc_vai - use EDTypesMod , only : nclmax, dlower_vai - use DamageMainMod , only : GetCrownReduction implicit none private - public :: GetNLevVeg ! Number of vegetation layers for the plant public :: h2d_allom ! Generic height to diameter wrapper public :: h_allom ! Generic diameter to height wrapper public :: bagw_allom ! Generic AGWB (above grnd. woody bio) wrapper @@ -126,7 +122,8 @@ module FatesAllometryMod public :: set_root_fraction ! Generic wrapper to calculate normalized ! root profiles public :: leafc_from_treelai ! Calculate target leaf carbon for a given treelai for SP mode - + public :: GetCrownReduction + logical , parameter :: verbose_logging = .false. character(len=*), parameter :: sourcefile = __FILE__ @@ -367,7 +364,6 @@ end subroutine h_allom subroutine bagw_allom(d,ipft,crowndamage, bagw,dbagwdd) - use DamageMainMod, only : GetCrownReduction use FatesParameterDerivedMod, only : param_derived real(r8),intent(in) :: d ! plant diameter [cm] @@ -421,6 +417,26 @@ subroutine bagw_allom(d,ipft,crowndamage, bagw,dbagwdd) return end subroutine bagw_allom + ! ============================================================================= + + subroutine GetCrownReduction(crowndamage, crown_reduction) + + use EDParamsMod , only : ED_val_history_damage_bin_edges + + !------------------------------------------------------------------ + ! This subroutine takes the crown damage class of a cohort (integer) + ! and returns the fraction of the crown that is lost. + !------------------------------------------------------------------- + + integer(i4), intent(in) :: crowndamage ! crown damage class of the cohort + real(r8), intent(out) :: crown_reduction ! fraction of crown lost from damage + + crown_reduction = ED_val_history_damage_bin_edges(crowndamage)/100.0_r8 + + return + end subroutine GetCrownReduction + + ! ============================================================================ ! Generic diameter to maximum leaf biomass interface ! ============================================================================ @@ -548,8 +564,6 @@ subroutine bleaf(d,ipft,crowndamage,canopy_trim,bl,dbldd) ! this routine is not name-spaced with allom_ ! ------------------------------------------------------------------------- - use DamageMainMod , only : GetCrownReduction - real(r8),intent(in) :: d ! plant diameter [cm] integer(i4),intent(in) :: ipft ! PFT index integer(i4),intent(in) :: crowndamage ! crown damage class [1: undamaged, >1: damaged] @@ -625,7 +639,7 @@ real(r8) function tree_lai( leaf_c, pft, c_area, nplant, cl, canopy_lai, vcmax25 real(r8), intent(in) :: c_area ! areal extent of canopy (m2) real(r8), intent(in) :: nplant ! number of individuals in cohort per ha integer, intent(in) :: cl ! canopy layer index - real(r8), intent(in) :: canopy_lai(nclmax) ! total leaf area index of + real(r8), intent(in) :: canopy_lai(:) ! total leaf area index of ! each canopy layer real(r8), intent(in) :: vcmax25top ! maximum carboxylation rate at canopy ! top, ref 25C @@ -742,7 +756,7 @@ end function tree_lai ! ============================================================================ real(r8) function tree_sai(pft, dbh, crowndamage, canopy_trim, c_area, nplant, cl, & - canopy_lai, treelai, vcmax25top, call_id ) + canopy_lai, treelai, max_treevai, vcmax25top, call_id ) ! ============================================================================ ! SAI of individual trees is a function of the LAI of individual trees @@ -751,16 +765,19 @@ real(r8) function tree_sai(pft, dbh, crowndamage, canopy_trim, c_area, nplant, c integer, intent(in) :: pft real(r8), intent(in) :: dbh integer, intent(in) :: crowndamage - real(r8), intent(in) :: canopy_trim ! trimming function (0-1) - real(r8), intent(in) :: c_area ! crown area (m2) - real(r8), intent(in) :: nplant ! number of plants - integer, intent(in) :: cl ! canopy layer index - real(r8), intent(in) :: canopy_lai(nclmax) ! total leaf area index of - ! each canopy layer - real(r8), intent(in) :: treelai ! tree LAI for checking purposes only - real(r8), intent(in) :: vcmax25top ! maximum carboxylation rate at top of crown - integer,intent(in) :: call_id ! flag specifying where this is called - ! from + real(r8), intent(in) :: canopy_trim ! trimming function (0-1) + real(r8), intent(in) :: c_area ! crown area (m2) + real(r8), intent(in) :: nplant ! number of plants + integer, intent(in) :: cl ! canopy layer index + real(r8), intent(in) :: canopy_lai(:) ! total leaf area index of + ! each canopy layer + real(r8), intent(in) :: treelai ! tree LAI for checking purposes only + real(r8), intent(in) :: max_treevai ! maximum tree LAI+SAI that is allowed based + ! on the constraints of the + ! incrementation arrays (ie dinc_vai) + real(r8), intent(in) :: vcmax25top ! maximum carboxylation rate at top of crown + integer,intent(in) :: call_id ! flag specifying where this is called + ! from real(r8) :: h real(r8) :: target_lai real(r8) :: target_bleaf @@ -772,7 +789,7 @@ real(r8) function tree_sai(pft, dbh, crowndamage, canopy_trim, c_area, nplant, c tree_sai = prt_params%allom_sai_scaler(pft) * target_lai - if( (treelai + tree_sai) > (sum(dinc_vai)) )then + if( (treelai + tree_sai) > max_treevai )then call h_allom(dbh,pft,h) @@ -783,8 +800,6 @@ real(r8) function tree_sai(pft, dbh, crowndamage, canopy_trim, c_area, nplant, c write(fates_log(),*) 'target_bleaf: ', target_bleaf write(fates_log(),*) 'area: ', c_area write(fates_log(),*) 'target_lai: ',target_lai - write(fates_log(),*) 'dinc_vai:',dinc_vai - write(fates_log(),*) 'nlevleaf,sum(dinc_vai):',nlevleaf,sum(dinc_vai) write(fates_log(),*) 'pft: ',pft write(fates_log(),*) 'call id: ',call_id write(fates_log(),*) 'n: ',nplant @@ -802,56 +817,6 @@ real(r8) function tree_sai(pft, dbh, crowndamage, canopy_trim, c_area, nplant, c return end function tree_sai - - subroutine GetNLevVeg(dbh, leaf_c, c_area, ipft, nplant, & - crowndamage, canopy_trim, canopy_layer, canopy_layer_tlai, & - vcmax25top, nveg_act, nveg_max) - - ! Determine the maximum number of leaf (vegetation) layers - ! for a cohort. This is based off of allometry, and assuming the plant - ! has the maximum leaf for its size and trimming. - - real(r8), intent(inout) :: dbh - real(r8), intent(in) :: leaf_c - real(r8), intent(in) :: c_area - integer, intent(in) :: ipft - real(r8), intent(in) :: nplant - real(r8), intent(in) :: canopy_trim - integer, intent(in) :: crowndamage - integer, intent(in) :: canopy_layer - real(r8), intent(in) :: canopy_layer_tlai(:) - real(r8), intent(in) :: vcmax25top - integer, intent(out) :: nveg_act - integer, intent(out) :: nveg_max - - integer :: nv - real(r8) :: leaf_c_target - real(r8) :: max_lai, max_sai - - - max_lai = tree_lai(leaf_c, ipft, c_area, & - nplant, canopy_layer, canopy_layer_tlai, vcmax25top ) - - max_sai = tree_sai(ipft, dbh, crowndamage, canopy_trim, & - c_area, nplant, canopy_layer, canopy_layer_tlai, max_lai, & - vcmax25top, 0) - - nveg_act = count((max_lai+max_sai) .gt. dlower_vai(:)) + 1 - - call bleaf(dbh,ipft,crowndamage,canopy_trim,leaf_c_target) - - max_lai = tree_lai(leaf_c_target, ipft, c_area, & - nplant, canopy_layer, canopy_layer_tlai, vcmax25top ) - - max_sai = tree_sai(ipft, dbh, crowndamage, canopy_trim, & - c_area, nplant, canopy_layer, canopy_layer_tlai, max_lai, & - vcmax25top, 0) - - nveg_max = count((max_lai+max_sai) .gt. dlower_vai(:)) + 1 - - return - end subroutine GetNLevVeg - ! ===================================================================================== real(r8) function leafc_from_treelai( treelai, pft, c_area, nplant, cl, vcmax25top) @@ -937,17 +902,12 @@ end function leafc_from_treelai ! ===================================================================================== - - - - ! ============================================================================ ! Generic sapwood biomass interface ! ============================================================================ subroutine bsap_allom(d,ipft,crowndamage,canopy_trim,sapw_area,bsap,dbsapdd) - use DamageMainMod , only : GetCrownReduction use FatesParameterDerivedMod, only : param_derived real(r8),intent(in) :: d ! plant diameter [cm] diff --git a/biogeophys/EDAccumulateFluxesMod.F90 b/biogeophys/EDAccumulateFluxesMod.F90 index 9aa50ab988..5cad3141e7 100644 --- a/biogeophys/EDAccumulateFluxesMod.F90 +++ b/biogeophys/EDAccumulateFluxesMod.F90 @@ -103,7 +103,7 @@ subroutine AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) if( ccohort%is_trimmable ) then - do iv=1,min(ccohort%nveg_act,nlevleafmem) + do iv=1,nlevleafmem ccohort%year_net_uptake(iv) = ccohort%year_net_uptake(iv) + ccohort%ts_net_uptake(iv) end do end if diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 1752e4e375..1763e315ee 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -39,6 +39,8 @@ module FATESPlantRespPhotosynthMod use EDTypesMod, only : nlevleaf use EDTypesMod, only : nlevleafmem use EDTypesMod, only : nclmax + use EDTypesMod, only : GetLeafFromMemLayer + use EDTypesMod, only : GetMemFromLeafLayer use PRTGenericMod, only : max_nleafage use EDTypesMod, only : do_fates_salinity use EDParamsMod, only : q10_mr @@ -136,8 +138,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) use FatesAllometryMod, only : storage_fraction_of_target use FatesAllometryMod, only : set_root_fraction use FatesAllometryMod, only : decay_coeff_kn - - use DamageMainMod, only : GetCrownReduction + use FatesAllometryMod, only : GetCrownReduction use FatesInterfaceTypesMod, only : hlm_use_tree_damage @@ -613,13 +614,13 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) end do leaf_layer_loop ! Zero cohort flux accumulators. - currentCohort%npp_tstep = 0.0_r8 - currentCohort%resp_tstep = 0.0_r8 - currentCohort%gpp_tstep = 0.0_r8 - currentCohort%rdark = 0.0_r8 - currentCohort%resp_m = 0.0_r8 - currentCohort%ts_net_uptake = 0.0_r8 - currentCohort%c13disc_clm = 0.0_r8 + currentCohort%npp_tstep = 0.0_r8 + currentCohort%resp_tstep = 0.0_r8 + currentCohort%gpp_tstep = 0.0_r8 + currentCohort%rdark = 0.0_r8 + currentCohort%resp_m = 0.0_r8 + currentCohort%ts_net_uptake(:) = 0.0_r8 + currentCohort%c13disc_clm = 0.0_r8 ! --------------------------------------------------------------- ! Part VII: Transfer leaf flux rates (like maintenance respiration, @@ -663,18 +664,32 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) if(.not.stretch_cbalmem_vert)then - ivb = nva - ivt = max(1,nvm-nlevleafmem+1) - if(ivb>=ivt) then - do iv = ivt,ivb - ivm = nvm - iv + 1 - ! Its possible that the leaves aren't flushed enough - ! to lign up with any of the memory bins - if(ivm>0 .and. ivm.le.nlevleafmem) then - currentCohort%ts_net_uptake(ivm) = anet_av_z(iv,ft,cl) * umolC_to_kgC - end if - end do - end if + ! These are the top (ivt) and bottom (ivb) + ! leaf-layers that match up with the top and bottom + ! leaf memory layers. iv loops over these leaf + ! layers, which starts higher (vertically) and decends + ! If the leaves are fully flushed, it should end on + ! the leaf memory index of 1. + + !ivb = nvm + !ivt = max(1,nvm-nlevleafmem+1) + !if(ivb>=ivt) then + ! do iv = ivt,ivb + ! ivm = currentCohort%GetMemFromLeafLayer(iv) + ! ! Its possible that the leaves aren't flushed enough + ! ! to lign up with any of the memory bins + ! if(ivm>0 .and. ivm.le.nlevleafmem) then + ! currentCohort%ts_net_uptake(ivm) = anet_av_z(iv,ft,cl) * umolC_to_kgC + ! end if + ! end do + !end if + + do ivm = 1,nlevleafmem + iv = GetLeafFromMemLayer(currentCohort,ivm) + if(iv>0 .and. iv<=nva) then + currentCohort%ts_net_uptake(ivm) = anet_av_z(iv,ft,cl) * umolC_to_kgC + end if + end do else @@ -687,7 +702,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! We are matching the memory bins with the actual ! bins by stretching out the actual bins to cover the same ! depth as the bin space defined by allometry - iv = min(nva,int(real(nvm - ivm + 1,r8) * flush_frac)) + iv = min(nva,int(real(GetLeafFromMemLayer(currentCohort,ivm),r8) * flush_frac)) currentCohort%ts_net_uptake(ivm) = flush_frac*anet_av_z(iv,ft,cl) * umolC_to_kgC diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 4c65ffe374..37e7e8971a 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -54,7 +54,7 @@ module EDInitMod use FatesInterfaceTypesMod , only : nlevdamage use FatesInterfaceTypesMod , only : hlm_use_nocomp use FatesInterfaceTypesMod , only : nlevage - + use EDCohortDynamicsMod , only : GetNLevVeg use FatesAllometryMod , only : h2d_allom use FatesAllometryMod , only : bagw_allom use FatesAllometryMod , only : bbgw_allom @@ -672,6 +672,7 @@ subroutine init_patches( nsites, sites, bc_in) do while(associated(newc)) newc%is_trimmable = .true. newc%year_net_uptake(:) = 0._r8 + call GetNLevVeg(newc,newp) newc => newc%shorter end do newp=>newp%younger diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index 7eb8985f87..8416a37002 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -86,7 +86,7 @@ module EDMainMod use FatesPlantHydraulicsMod , only : InitPlantHydStates use FatesPlantHydraulicsMod , only : UpdateSizeDepRhizHydProps use FatesPlantHydraulicsMod , only : AccumulateMortalityWaterStorage - use FatesAllometryMod , only : h_allom,tree_sai,tree_lai + use FatesAllometryMod , only : h_allom use EDLoggingMortalityMod , only : IsItLoggingTime use DamageMainMod , only : IsItDamageTime use EDPatchDynamicsMod , only : get_frac_site_primary diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index b3ec36ecc7..18ec019349 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -6,8 +6,7 @@ module EDPftvarcon ! read and initialize vegetation (PFT) constants. ! ! !USES: - use EDTypesMod , only : maxSWb, ivis, inir - use EDTypesMod , only : n_uptake_mode, p_uptake_mode + use EDTypesMod, only : maxSWb, ivis, inir use FatesConstantsMod, only : r8 => fates_r8 use FatesConstantsMod, only : nearzero use FatesConstantsMod, only : itrue, ifalse diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index d58a7c83f6..bba121a13a 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -4,6 +4,7 @@ module EDTypesMod use FatesConstantsMod, only : ifalse use FatesConstantsMod, only : itrue use FatesGlobals, only : fates_log + use FatesGlobals, only : endrun => fates_endrun use FatesHydraulicsMemMod, only : ed_cohort_hydr_type use FatesHydraulicsMemMod, only : ed_site_hydr_type use PRTGenericMod, only : prt_vartypes @@ -22,11 +23,17 @@ module EDTypesMod use FatesRunningMeanMod, only : rmean_type use FatesInterfaceTypesMod,only : bc_in_type use FatesInterfaceTypesMod,only : bc_out_type + use shr_log_mod ,only : errMsg => shr_log_errMsg implicit none private ! By default everything is private save + character(len=*), parameter, private :: sourcefile = & + __FILE__ + + logical, parameter :: debug = .true. ! Local debug flag + integer, parameter, public :: nclmax = 2 ! Maximum number of canopy layers integer, parameter, public :: ican_upper = 1 ! Nominal index for the upper canopy integer, parameter, public :: ican_ustory = 2 ! Nominal index for diagnostics that refer @@ -41,6 +48,8 @@ module EDTypesMod real(r8), parameter, public :: init_recruit_trim = 0.8_r8 ! This is the initial trimming value that ! new recruits start with + logical, parameter, public :: use_relative_leafmem_layers = .false. + ! ------------------------------------------------------------------------------------- ! Radiation parameters ! These should be part of the radiation module, but since we only have one option @@ -54,36 +63,11 @@ module EDTypesMod integer, parameter, public :: idiffuse = 2 ! This is the array index for diffuse radiation ! parameters that govern the VAI (LAI+SAI) bins used in radiative transfer code - integer, parameter, public :: nlevleafmem = 5 ! number of layers for veg memory used for trimming + integer, parameter, public :: nlevleafmem = 7 ! number of layers for veg memory used for trimming integer, parameter, public :: nlevleaf = 30 ! number of leaf+stem layers in canopy layer real(r8), public :: dinc_vai(nlevleaf) = fates_unset_r8 ! VAI bin widths array real(r8), public :: dlower_vai(nlevleaf) = fates_unset_r8 ! lower edges of VAI bins - ! TODO: we use this cp_maxSWb only because we have a static array q(size=2) of - ! land-ice abledo for vis and nir. This should be a parameter, which would - ! get us on track to start using multi-spectral or hyper-spectral (RGK 02-2017) - - integer, parameter, public :: maxSWb = 2 ! maximum number of broad-bands in the - ! shortwave spectrum cp_numSWb <= cp_maxSWb - ! this is just for scratch-array purposes - ! if cp_numSWb is larger than this value - ! simply bump this number up as needed - - integer, parameter, public :: ivis = 1 ! This is the array index for short-wave - ! radiation in the visible spectrum, as expected - ! in boundary condition files and parameter - ! files. This will be compared with - ! the HLM's expectation in FatesInterfaceMod - integer, parameter, public :: inir = 2 ! This is the array index for short-wave - ! radiation in the near-infrared spectrum, as expected - ! in boundary condition files and parameter - ! files. This will be compared with - ! the HLM's expectation in FatesInterfaceMod - - integer, parameter, public :: ipar = ivis ! The photosynthetically active band - ! can be approximated to be equal to the visible band - - integer, parameter, public :: leaves_on = 2 ! Flag specifying that a deciduous plant has leaves ! and should be allocating to them as well integer, parameter, public :: leaves_off = 1 ! Flag specifying that a deciduous plant has dropped @@ -119,7 +103,29 @@ module EDTypesMod integer, parameter, public :: numlevsoil_max = 30 ! This is scratch space used for static arrays ! The actual number of soil layers should not exceed this + ! TODO: we use this cp_maxSWb only because we have a static array q(size=2) of + ! land-ice abledo for vis and nir. This should be a parameter, which would + ! get us on track to start using multi-spectral or hyper-spectral (RGK 02-2017) + + integer, parameter, public :: maxSWb = 2 ! maximum number of broad-bands in the + ! shortwave spectrum cp_numSWb <= cp_maxSWb + ! this is just for scratch-array purposes + ! if cp_numSWb is larger than this value + ! simply bump this number up as needed + + integer, parameter, public :: ivis = 1 ! This is the array index for short-wave + ! radiation in the visible spectrum, as expected + ! in boundary condition files and parameter + ! files. This will be compared with + ! the HLM's expectation in FatesInterfaceMod + integer, parameter, public :: inir = 2 ! This is the array index for short-wave + ! radiation in the near-infrared spectrum, as expected + ! in boundary condition files and parameter + ! files. This will be compared with + ! the HLM's expectation in FatesInterfaceMod + integer, parameter, public :: ipar = ivis ! The photosynthetically active band + ! can be approximated to be equal to the visible band ! BIOLOGY/BIOGEOCHEMISTRY integer , parameter, public :: num_vegtemp_mem = 10 ! Window of time over which we track temp for cold sensecence (days) @@ -140,15 +146,13 @@ module EDTypesMod integer, parameter, public :: phen_dstat_moiston = 2 ! Leaves on due to moisture avail (drought phenology) integer, parameter, public :: phen_dstat_timeon = 3 ! Leaves on due to time exceedance (drought phenology) + integer, parameter, public :: NFSC = NCWD+2 ! number fuel size classes (4 cwd size classes, leaf litter, and grass) + integer, parameter, public :: tw_sf = 1 ! array index of twig pool for spitfire + integer, parameter, public :: lb_sf = 3 ! array index of large branch pool for spitfire + integer, parameter, public :: tr_sf = 4 ! array index of dead trunk pool for spitfire + integer, parameter, public :: dl_sf = 5 ! array index of dead leaf pool for spitfire (dead grass and dead leaves) + integer, parameter, public :: lg_sf = 6 ! array index of live grass pool for spitfire - ! SPITFIRE - - integer, parameter, public :: NFSC = NCWD+2 ! number fuel size classes (4 cwd size classes, leaf litter, and grass) - integer, parameter, public :: tw_sf = 1 ! array index of twig pool for spitfire - integer, parameter, public :: lb_sf = 3 ! array index of large branch pool for spitfire - integer, parameter, public :: tr_sf = 4 ! array index of dead trunk pool for spitfire - integer, parameter, public :: dl_sf = 5 ! array index of dead leaf pool for spitfire (dead grass and dead leaves) - integer, parameter, public :: lg_sf = 6 ! array index of live grass pool for spitfire ! PATCH FUSION real(r8), parameter, public :: force_patchfuse_min_biomass = 0.005_r8 ! min biomass (kg / m2 patch area) below which to force-fuse patches @@ -235,7 +239,10 @@ module EDTypesMod real(r8) :: excl_weight ! How much of this cohort is demoted each year, as a proportion of all cohorts:- real(r8) :: prom_weight ! How much of this cohort is promoted each year, as a proportion of all cohorts:- integer :: nveg_act ! Number of vegetation layers (actual) at any point in time - integer :: nveg_max ! Maximum number of vegetation layers based on allometry + integer :: nveg_max ! Maximum number of vegetation layers based on allometry + !integer :: nveg_mem ! The vegetation layer that matches memory layer 1, we set this layer + ! lower down in the vegetation layers, in-case the plant grows beyond + ! its current nveg_max in between trimming events integer :: status_coh ! growth status of plant (2 = leaves on , 1 = leaves off) real(r8) :: c_area ! areal extent of canopy (m2) real(r8) :: treelai ! lai of an individual within cohort leaf area (m2) / crown area (m2) @@ -414,9 +421,11 @@ module EDTypesMod !class(rmean_type), pointer :: tveg_lpa ! exponential moving average of leaf temperature at the ! leaf photosynthetic acclimation time-scale [K] - end type ed_cohort_type + + + !************************************ !** Patch type structure ** !************************************ @@ -906,12 +915,17 @@ module EDTypesMod public :: dump_patch public :: dump_cohort public :: dump_cohort_hydr - public :: CanUpperUnder - contains + public :: CanUpperUnder + public :: GetMemFromLeafLayer + public :: GetLeafFromMemLayer + +contains - ! ===================================================================================== - function CanUpperUnder(ccohort) result(can_position) + + ! ===================================================================================== + + function CanUpperUnder(ccohort) result(can_position) ! This simple function is used to determine if a ! cohort's crown position is in the upper portion (ie the canopy) @@ -931,7 +945,12 @@ function CanUpperUnder(ccohort) result(can_position) end if end function CanUpperUnder - + + + + + + ! ===================================================================================== subroutine ZeroFluxDiags(this) @@ -1211,4 +1230,29 @@ subroutine dump_cohort_hydr(ccohort) return end subroutine dump_cohort_hydr + ! ===================================================================================== + + function GetMemFromLeafLayer(cohort,ileaf) result(imem) + + type(ed_cohort_type) :: cohort + integer :: ileaf + integer :: imem + + imem = cohort%nveg_max - ileaf + 1 + + end function GetMemFromLeafLayer + + ! ===================================================================================== + + function GetLeafFromMemLayer(cohort,imem) result(ileaf) + + type(ed_cohort_type) :: cohort + integer :: ileaf + integer :: imem + + ileaf = cohort%nveg_max - imem + 1 + + end function GetLeafFromMemLayer + + end module EDTypesMod diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index fb132a99d6..ee3bd6cfb3 100644 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -53,17 +53,7 @@ module FatesConstantsMod ! plants are not coupled with below ground chemistry. In ! this situation, we send token boundary condition information. - - integer, public, parameter :: trim_day_of_year = 364 ! Trim on this date (2nd to last day of year) - - ! If this is a cold start and the date is after the date - ! we trim on, but fairly close, then we say its close - ! enough and flag the newly created cohorts as trimmable - ! the next time we encounter trimming, and zero out - ! the incrementing carbon balance array - integer, public, parameter :: max_days_since_trim_coldstart = 5 - ! This flag specifies the scaling of how we present ! nutrient competitors to the HLM's soil BGC model diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 64a671bac4..908019b7be 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -1848,7 +1848,6 @@ subroutine set_restart_vectors(this,nc,nsites,sites) integer :: io_idx_si_pft ! each site-pft index integer :: io_idx_si_vtmem ! indices for veg-temp memory at site integer :: io_idx_pa_ncl ! each canopy layer within each patch - integer :: io_idx_co_nll ! each leaf layer within each cohort ! Some counters (for checking mostly) integer :: totalcohorts ! total cohort count on this thread (diagnostic) @@ -2005,7 +2004,6 @@ subroutine set_restart_vectors(this,nc,nsites,sites) io_idx_si_pfcl = io_idx_co_1st io_idx_si_vtmem = io_idx_co_1st io_idx_pa_ncl = io_idx_co_1st - io_idx_co_nll = io_idx_co_1st ! Hydraulics counters lyr = hydraulic layer, shell = rhizosphere shell io_idx_si_lyr_shell = io_idx_co_1st @@ -2754,7 +2752,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) integer :: io_idx_si_cdpf ! damage x size x pft within site integer :: io_idx_pa_ncl ! each canopy layer within each patch - integer :: io_idx_co_nll ! each leaf layer within each cohort ! Some counters (for checking mostly) integer :: totalcohorts ! total cohort count on this thread (diagnostic) @@ -2892,7 +2889,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) io_idx_si_pfcl = io_idx_co_1st io_idx_si_vtmem = io_idx_co_1st io_idx_pa_ncl = io_idx_co_1st - io_idx_co_nll = io_idx_co_1st ! Hydraulics counters lyr = hydraulic layer, shell = rhizosphere shell io_idx_si_lyr_shell = io_idx_co_1st From 378b7d8beb882c6b65af034fdd4b3abae6d89a15 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Sun, 18 Dec 2022 23:30:28 -0500 Subject: [PATCH 08/12] small fix --- biogeochem/EDPhysiologyMod.F90 | 2 +- main/FatesConstantsMod.F90 | 10 ++++++++++ 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 781501084f..775ca02553 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -779,7 +779,7 @@ subroutine trim_canopy( currentSite ) ! if at least nll leaf layers are present in the current cohort and only for the bottom nll ! leaf layers. - if( currentCohort%nveg_max >= nll ) then + if( currentCohort%nveg_max >= nll .and. zm <= nll) then ! Build the A matrix for the LHS of the linear system. A = [n sum(x); sum(x) sum(x^2)] ! where n = nll and x = yearly_net_uptake-leafcost diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index ee3bd6cfb3..fb132a99d6 100644 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -53,7 +53,17 @@ module FatesConstantsMod ! plants are not coupled with below ground chemistry. In ! this situation, we send token boundary condition information. + + integer, public, parameter :: trim_day_of_year = 364 ! Trim on this date (2nd to last day of year) + + ! If this is a cold start and the date is after the date + ! we trim on, but fairly close, then we say its close + ! enough and flag the newly created cohorts as trimmable + ! the next time we encounter trimming, and zero out + ! the incrementing carbon balance array + integer, public, parameter :: max_days_since_trim_coldstart = 5 + ! This flag specifies the scaling of how we present ! nutrient competitors to the HLM's soil BGC model From f2d4ed1bbda920ea9ce346613f90adff58f49236 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 20 Dec 2022 16:22:14 -0500 Subject: [PATCH 09/12] Updates to the veg memory reductions --- biogeochem/EDCohortDynamicsMod.F90 | 20 +++++++++++++++----- main/EDTypesMod.F90 | 3 --- 2 files changed, 15 insertions(+), 8 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 417844c0c6..0bb157daef 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -2015,11 +2015,18 @@ subroutine FuseVegLayerMem(ccohort,ncohort,cpatch,site_spread) type(ed_patch_type) :: cpatch ! current patch (both cohorts are on this patch) real(r8),intent(in) :: site_spread - real(r8) :: joint_net_uptake(nlevleaf) - real(r8) :: joint_wgt(nlevleaf) - integer :: iv ! veg layer index - integer :: ivm ! veg layer memory index - integer :: iv0,iv1 ! lowest and highest veg layer index matching memory + real(r8) :: joint_net_uptake(nlevleaf) ! temp space for joint carbon balance data + real(r8) :: joint_wgt(nlevleaf) ! weighting factor (plant num) for joint data + integer :: iv ! veg layer index + integer :: ivm ! veg layer memory index + integer :: iv0,iv1 ! lowest and highest veg layer index matching memory + real(r8) :: nplant ! We save the old number of plants in the + ! receiver cohort. We do this because we want to pass the joint + ! number to a routine, yet we want to revert to the unfused + ! number again in case this routine is not the last thing + ! called during cohort fusion, because the number + ! density is fused last. + joint_net_uptake(:) = 0._r8 joint_wgt(:) = 0._r8 @@ -2101,7 +2108,10 @@ subroutine FuseVegLayerMem(ccohort,ncohort,cpatch,site_spread) ! Update the number of vegetation layers, cohort%nveg_act and cohort%nveg_max + nplant = ccohort%n + ccohort%n = ccohort%n + ncohort%n call GetNLevVeg(ccohort,cpatch) + ccohort%n = nplant end subroutine FuseVegLayerMem diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index bba121a13a..9ae53b1114 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -240,9 +240,6 @@ module EDTypesMod real(r8) :: prom_weight ! How much of this cohort is promoted each year, as a proportion of all cohorts:- integer :: nveg_act ! Number of vegetation layers (actual) at any point in time integer :: nveg_max ! Maximum number of vegetation layers based on allometry - !integer :: nveg_mem ! The vegetation layer that matches memory layer 1, we set this layer - ! lower down in the vegetation layers, in-case the plant grows beyond - ! its current nveg_max in between trimming events integer :: status_coh ! growth status of plant (2 = leaves on , 1 = leaves off) real(r8) :: c_area ! areal extent of canopy (m2) real(r8) :: treelai ! lai of an individual within cohort leaf area (m2) / crown area (m2) From 26cadf4bd32c07da005c316b10dac1dbcaf53c6a Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 21 Dec 2022 15:47:19 -0500 Subject: [PATCH 10/12] final tweaks to reduced memory year_net_uptake arrays that gets qualitatively similar results as main --- biogeochem/EDCohortDynamicsMod.F90 | 69 +++++++++++++++++-------- biogeochem/EDPhysiologyMod.F90 | 77 ++++++++++++++++++---------- biogeophys/EDAccumulateFluxesMod.F90 | 13 ++++- main/EDTypesMod.F90 | 10 +++- main/FatesHistoryInterfaceMod.F90 | 26 +++++----- 5 files changed, 133 insertions(+), 62 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 0bb157daef..78e83e5446 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -28,6 +28,7 @@ Module EDCohortDynamicsMod use FatesParameterDerivedMod, only : param_derived use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use EDTypesMod , only : nclmax + use EDTypesMod , only : match_old_trim_method use PRTGenericMod , only : element_list use PRTGenericMod , only : StorageNutrientTarget use FatesLitterMod , only : ncwd @@ -39,7 +40,6 @@ Module EDCohortDynamicsMod use EDTypesMod , only : min_n_safemath use EDTypesMod , only : nlevleaf use EDTypesMod , only : nlevleafmem - use EDTypesMod , only : use_relative_leafmem_layers use EDTypesMod , only : GetLeafFromMemLayer use PRTGenericMod , only : max_nleafage use EDTypesMod , only : ican_upper @@ -2031,15 +2031,18 @@ subroutine FuseVegLayerMem(ccohort,ncohort,cpatch,site_spread) joint_net_uptake(:) = 0._r8 joint_wgt(:) = 0._r8 - ccohort%canopy_trim = (ccohort%n*ccohort%canopy_trim & - + ncohort%n*ncohort%canopy_trim)/(ccohort%n+ncohort%n) - + if(match_old_trim_method) then + ccohort%canopy_trim = (ccohort%n*ccohort%canopy_trim & + + ncohort%n*ncohort%canopy_trim)/(ccohort%n+ncohort%n) + end if if(ccohort%is_trimmable .and. ncohort%is_trimmable) then - !ccohort%canopy_trim = (ccohort%n*ccohort%canopy_trim & - ! + ncohort%n*ncohort%canopy_trim)/(ccohort%n+ncohort%n) - + if(.not.match_old_trim_method) then + ccohort%canopy_trim = (ccohort%n*ccohort%canopy_trim & + + ncohort%n*ncohort%canopy_trim)/(ccohort%n+ncohort%n) + end if + do ivm = 1, min(nlevleafmem,ccohort%nveg_max) iv = GetLeafFromMemLayer(ccohort,ivm) joint_net_uptake(iv) = ccohort%n*ccohort%year_net_uptake(ivm) @@ -2099,7 +2102,7 @@ subroutine FuseVegLayerMem(ccohort,ncohort,cpatch,site_spread) ! thus, nothing needs to be done !!! if(ncohort%is_trimmable) then - !ccohort%canopy_trim = ncohort%canopy_trim + if (.not.match_old_trim_method) ccohort%canopy_trim = ncohort%canopy_trim ccohort%year_net_uptake(:) = ncohort%year_net_uptake(:) ccohort%is_trimmable = .true. end if @@ -2573,6 +2576,16 @@ subroutine GetNLevVeg(cohort,patch) ! net carbon uptake array, if necessary integer :: di ! differential memory index used for ! shifting yearly uptake array to a new offset + + ! If the max layer has changed, then + ! we need to shift the year_net_uptake + ! array. These holds the previous + ! values in the vertically lowest (1) + ! and highest (nlevleafmem) positions + ! which may be used to initialize new positions + real(r8) :: year_net_uptake_low + real(r8) :: year_net_uptake_high + lai = tree_lai(cohort%prt%GetState(leaf_organ, carbon12_element), & cohort%pft, cohort%c_area, & @@ -2598,20 +2611,36 @@ subroutine GetNLevVeg(cohort,patch) nveg_max = count((lai+sai) .gt. dlower_vai(:)) + 1 - ! If the plant recently grew past the matchpoint on the veg layer memory - ! array, then we have to shift everything up one and initialize the new - ! memory layer. - - if( (.not.use_relative_leafmem_layers) .and. & - (cohort%nveg_max.ne.fates_unset_int) ) then - if(nveg_max .ne. cohort%nveg_max) then + ! If the plant recently changed its nveg_max + ! then the matchpoint on the veg layer memory + ! array is different and the array needs to shift + + if (cohort%nveg_max.ne.fates_unset_int) then + if (nveg_max .ne. cohort%nveg_max) then di = nveg_max - cohort%nveg_max - if(di>0) then - cohort%year_net_uptake(1+di:nlevleafmem) = cohort%year_net_uptake(1:nlevleafmem-di) - cohort%year_net_uptake(1:di) = 0._r8 + + year_net_uptake_low = cohort%year_net_uptake(1) + year_net_uptake_high = cohort%year_net_uptake(nlevleafmem) + + if (match_old_trim_method) then + if(di>0) then + ! We need a new layer + cohort%year_net_uptake(1:di) = 0._r8 + cohort%year_net_uptake(1+di:nlevleafmem) = cohort%year_net_uptake(1:nlevleafmem-di) + else + ! We need to remove a layer + cohort%year_net_uptake(nlevleafmem+di+1:nlevleafmem) = year_net_uptake_high + cohort%year_net_uptake(1:nlevleafmem+di) = cohort%year_net_uptake(1-di:nlevleafmem) + end if else - cohort%year_net_uptake(1:nlevleafmem+di) = cohort%year_net_uptake(1-di:nlevleafmem) - cohort%year_net_uptake(nlevleafmem+di+1:nlevleafmem) = 0._r8 + if(di>0) then + ! We need a new layer + cohort%year_net_uptake(1+di:nlevleafmem) = cohort%year_net_uptake(1:nlevleafmem-di) + cohort%year_net_uptake(1:di) = year_net_uptake_low + else + cohort%year_net_uptake(1:nlevleafmem+di) = cohort%year_net_uptake(1-di:nlevleafmem) + cohort%year_net_uptake(nlevleafmem+di+1:nlevleafmem) = year_net_uptake_high + end if end if end if end if diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 775ca02553..805fe257d9 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -43,8 +43,9 @@ module EDPhysiologyMod use EDTypesMod , only : numlevsoil_max use EDTypesMod , only : numWaterMem use EDTypesMod , only : dl_sf, dinc_vai, dlower_vai, area_inv - use EDTypesMod , only : GetLeafFromMemLayer + use EDTypesMod , only : GetLeafFromMemLayer,GetMemFromLeafLayer use EDTypesMod , only : AREA + use EDTypesMod , only : match_old_trim_method use FatesLitterMod , only : ncwd use FatesLitterMod , only : ndcmpy use FatesLitterMod , only : ilabile @@ -72,7 +73,6 @@ module EDPhysiologyMod use EDTypesMod , only : phen_dstat_moiston use EDTypesMod , only : phen_dstat_timeon use EDTypesMod , only : init_recruit_trim - use EDTypesMod , only : use_relative_leafmem_layers use shr_log_mod , only : errMsg => shr_log_errMsg use FatesGlobals , only : fates_log use FatesGlobals , only : endrun => fates_endrun @@ -600,26 +600,30 @@ subroutine trim_canopy( currentSite ) ! b is the y-intercept i.e. the cumulative lai that has zero net-net uptake ! m is the slope of the linear fit - integer :: nll = 3 ! Number of leaf layers to fit a regression - ! to for calculating the optimum lai - character(1) :: trans = 'N' ! Input matrix is not transposed + integer :: nll = 3 ! Number of leaf layers to fit a regression + ! to for calculating the optimum lai + character(1) :: trans = 'N' ! Input matrix is not transposed - integer, parameter :: m = 2, n = 2 ! Number of rows and columns, respectively, in matrix A - integer, parameter :: nrhs = 1 ! Number of columns in matrix B and X - integer, parameter :: workmax = 100 ! Maximum iterations to minimize work + integer, parameter :: m = 2, n = 2 ! Number of rows and columns, respectively, in matrix A + integer, parameter :: nrhs = 1 ! Number of columns in matrix B and X + integer, parameter :: workmax = 100 ! Maximum iterations to minimize work - integer :: lda = m, ldb = n ! Leading dimension of A and B, respectively - integer :: lwork ! Dimension of work array - integer :: info ! Procedure diagnostic ouput - - real(r8) :: nnu_clai_a(m,n) ! LHS of linear least squares fit, A matrix - real(r8) :: nnu_clai_b(m,nrhs) ! RHS of linear least squares fit, B matrix - real(r8) :: work(workmax) ! work array - - real(r8) :: initial_trim ! Initial trim - real(r8) :: optimum_trim ! Optimum trim value - - real(r8) :: target_c_area + integer :: lda = m, ldb = n ! Leading dimension of A and B, respectively + integer :: lwork ! Dimension of work array + integer :: info ! Procedure diagnostic ouput + integer :: z0,z1 ! Veg layer loop bounds + + real(r8) :: nnu_clai_a(m,n) ! LHS of linear least squares fit, A matrix + real(r8) :: nnu_clai_b(m,nrhs) ! RHS of linear least squares fit, B matrix + real(r8) :: work(workmax) ! work array + + real(r8) :: initial_trim ! Initial trim + real(r8) :: optimum_trim ! Optimum trim value + logical :: build_matrix ! Logical used to decid if we build out the + ! trim optimization matrix, which + ! depends on which version of trimming + ! (new vs old) we use + !---------------------------------------------------------------------- @@ -705,11 +709,19 @@ subroutine trim_canopy( currentSite ) cumulative_lai_cohort = 0._r8 - leafmem_loop: do zm = min(nlevleafmem,currentCohort%nveg_max),1,-1 - - ! The actual leaf layer index (going from top to bottom) - z = GetLeafFromMemLayer(currentCohort,zm) + if(match_old_trim_method) then + z0 = 1 + z1 = currentCohort%nveg_act + else + z0 = GetLeafFromMemLayer(currentCohort,min(nlevleafmem,currentCohort%nveg_max)) + z1 = GetLeafFromMemLayer(currentCohort,1) + end if + leafmem_loop: do z = z0,z1 + + ! Get the memory layer index + zm = GetMemFromLeafLayer(currentCohort,z) + ! Calculate the cumulative total vegetation area index (no snow occlusion, stems and leaves) leaf_inc = dinc_vai(z) * & @@ -775,11 +787,24 @@ subroutine trim_canopy( currentSite ) (prt_params%grperc(ipft) + 1._r8) endif + if(match_old_trim_method) then + if (currentCohort%nveg_act > nll .and. currentCohort%nveg_act - z < nll) then + build_matrix = .true. + else + build_matrix = .false. + end if + else + if( currentCohort%nveg_max >= nll .and. zm <= nll) then + build_matrix = .true. + else + build_matrix = .false. + end if + end if + ! Construct the arrays for a least square fit of the net_net_uptake versus the cumulative lai ! if at least nll leaf layers are present in the current cohort and only for the bottom nll ! leaf layers. - - if( currentCohort%nveg_max >= nll .and. zm <= nll) then + if (build_matrix) then ! Build the A matrix for the LHS of the linear system. A = [n sum(x); sum(x) sum(x^2)] ! where n = nll and x = yearly_net_uptake-leafcost diff --git a/biogeophys/EDAccumulateFluxesMod.F90 b/biogeophys/EDAccumulateFluxesMod.F90 index 5cad3141e7..db5866c787 100644 --- a/biogeophys/EDAccumulateFluxesMod.F90 +++ b/biogeophys/EDAccumulateFluxesMod.F90 @@ -14,7 +14,7 @@ module EDAccumulateFluxesMod use shr_log_mod , only : errMsg => shr_log_errMsg use FatesConstantsMod , only : r8 => fates_r8 use EDTypesMod, only : nlevleafmem - + use EDTypesMod, only : match_old_trim_method implicit none private @@ -101,7 +101,16 @@ subroutine AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) (ccohort%gpp_acc + ccohort%gpp_tstep) endif - + ! Force v1 version of trimming memory + if(match_old_trim_method) then + if( .not.ccohort%is_trimmable ) then + ccohort%is_trimmable = .true. + do iv=1,nlevleafmem + ccohort%year_net_uptake(iv) = 0._r8 + end do + end if + end if + if( ccohort%is_trimmable ) then do iv=1,nlevleafmem ccohort%year_net_uptake(iv) = ccohort%year_net_uptake(iv) + ccohort%ts_net_uptake(iv) diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 9ae53b1114..d4d05e3a58 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -48,7 +48,13 @@ module EDTypesMod real(r8), parameter, public :: init_recruit_trim = 0.8_r8 ! This is the initial trimming value that ! new recruits start with - logical, parameter, public :: use_relative_leafmem_layers = .false. + + logical, parameter, public :: match_old_trim_method = .true. ! There are some subtle + ! updates to trimming logic. For + ! instance, does a cohort need to + ! have a year of year_net_uptake + ! data to be allowed to trim (new) + ! or not (old) ! ------------------------------------------------------------------------------------- ! Radiation parameters @@ -63,7 +69,7 @@ module EDTypesMod integer, parameter, public :: idiffuse = 2 ! This is the array index for diffuse radiation ! parameters that govern the VAI (LAI+SAI) bins used in radiative transfer code - integer, parameter, public :: nlevleafmem = 7 ! number of layers for veg memory used for trimming + integer, parameter, public :: nlevleafmem = 5 ! number of layers for veg memory used for trimming integer, parameter, public :: nlevleaf = 30 ! number of leaf+stem layers in canopy layer real(r8), public :: dinc_vai(nlevleaf) = fates_unset_r8 ! VAI bin widths array real(r8), public :: dlower_vai(nlevleaf) = fates_unset_r8 ! lower edges of VAI bins diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index d0ceafaf2c..dec7e978ae 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -14,6 +14,7 @@ module FatesHistoryInterfaceMod use FatesGlobals , only : endrun => fates_endrun use EDTypesMod , only : nclmax use EDTypesMod , only : ican_upper + use EDTypesMod , only : GetMemFromLeafLayer use PRTGenericMod , only : element_pos use PRTGenericMod , only : num_elements use PRTGenericMod , only : prt_cnp_flex_allom_hyp @@ -4442,19 +4443,20 @@ subroutine update_history_hifrq(this,nc,nsites,sites,bc_in,dt_tstep) endif end associate endif - - !!! canopy leaf carbon balance - ican = ccohort%canopy_layer - do ileaf=1, min(nlevleafmem,ccohort%nveg_max) - - ! The actual leaf layer index (going from top to bottom) - z = ccohort%nveg_max - ileaf + 1 - cnlf_indx = z + (ican-1) * nlevleaf + + ! canopy leaf carbon balance + if(nlevleafmem == nlevleaf) then + ican = ccohort%canopy_layer + do z = 1,ccohort%nveg_act + ! The actual leaf layer index (going from top to bottom) + ileaf = GetMemFromLeafLayer(ccohort,z) + cnlf_indx = z + (ican-1) * nlevleaf + + hio_ts_net_uptake_si_cnlf(io_si, cnlf_indx) = hio_ts_net_uptake_si_cnlf(io_si, cnlf_indx) + & + ccohort%ts_net_uptake(ileaf) * per_dt_tstep * ccohort%c_area * area_inv + end do + end if - hio_ts_net_uptake_si_cnlf(io_si, cnlf_indx) = hio_ts_net_uptake_si_cnlf(io_si, cnlf_indx) + & - ccohort%ts_net_uptake(ileaf) * per_dt_tstep * ccohort%c_area * area_inv - end do - ccohort => ccohort%taller enddo ! cohort loop From 225adc5c259acf0750c7bd27f11de727597a70b4 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 21 Dec 2022 16:24:54 -0500 Subject: [PATCH 11/12] removed old comments --- biogeophys/FatesPlantRespPhotosynthMod.F90 | 17 ----------------- 1 file changed, 17 deletions(-) diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 1763e315ee..d6572b439d 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -664,25 +664,8 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) if(.not.stretch_cbalmem_vert)then - ! These are the top (ivt) and bottom (ivb) - ! leaf-layers that match up with the top and bottom - ! leaf memory layers. iv loops over these leaf - ! layers, which starts higher (vertically) and decends ! If the leaves are fully flushed, it should end on ! the leaf memory index of 1. - - !ivb = nvm - !ivt = max(1,nvm-nlevleafmem+1) - !if(ivb>=ivt) then - ! do iv = ivt,ivb - ! ivm = currentCohort%GetMemFromLeafLayer(iv) - ! ! Its possible that the leaves aren't flushed enough - ! ! to lign up with any of the memory bins - ! if(ivm>0 .and. ivm.le.nlevleafmem) then - ! currentCohort%ts_net_uptake(ivm) = anet_av_z(iv,ft,cl) * umolC_to_kgC - ! end if - ! end do - !end if do ivm = 1,nlevleafmem iv = GetLeafFromMemLayer(currentCohort,ivm) From 6755cf7bed340b3ba4fbdaba2e2c358fefc86569 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Sun, 1 Jan 2023 20:26:42 -0700 Subject: [PATCH 12/12] Removed clause that disallowed zero weight when fusing leaf layer carbon balance --- biogeochem/EDCohortDynamicsMod.F90 | 17 +++-------------- 1 file changed, 3 insertions(+), 14 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 78e83e5446..6908410652 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -2067,20 +2067,9 @@ subroutine FuseVegLayerMem(ccohort,ncohort,cpatch,site_spread) end if do iv = iv0 ,iv1 - if(joint_wgt(iv)nearzero) then joint_net_uptake(iv) = joint_net_uptake(iv)/joint_wgt(iv) end if end do