diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 3de3ceb94e..e113874046 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1018,7 +1018,7 @@ subroutine canopy_summarization( nsites, sites, bc_in ) currentCohort%pft,currentCohort%c_area) currentCohort%treelai = tree_lai(currentCohort%bl, currentCohort%status_coh, & currentCohort%pft, currentCohort%c_area, currentCohort%n, & - currentCohort%canopy_layer, currentPatch%canopy_layer_tai ) + currentCohort%canopy_layer, currentPatch%canopy_layer_tvai ) canopy_leaf_area = canopy_leaf_area + currentCohort%treelai *currentCohort%c_area @@ -1088,7 +1088,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) ! ! The following patch level diagnostics are updated here: ! - ! currentPatch%canopy_layer_tai(cl) ! TAI of each canopy layer + ! currentPatch%canopy_layer_tvai(cl) ! total vegetated area index of layer ! currentPatch%ncan(cl,ft) ! number of vegetation layers needed ! ! in this patch's pft/canopy-layer ! currentPatch%nrad(cl,ft) ! same as ncan, but does not include @@ -1158,7 +1158,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) ! calculate tree lai and sai. ! -------------------------------------------------------------------------------- - currentPatch%canopy_layer_tai(:) = 0._r8 + currentPatch%canopy_layer_tvai(:) = 0._r8 currentPatch%ncan(:,:) = 0 currentPatch%nrad(:,:) = 0 patch_lai = 0._r8 @@ -1184,25 +1184,25 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) cl = currentCohort%canopy_layer currentCohort%treelai = tree_lai(currentCohort%bl, currentCohort%status_coh, currentCohort%pft, & - currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, currentPatch%canopy_layer_tai ) + currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, currentPatch%canopy_layer_tvai ) currentCohort%treesai = tree_sai(currentCohort%pft, currentCohort%treelai) 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) + 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_tai(cl) = currentPatch%canopy_layer_tai(cl) + & +! currentPatch%canopy_layer_tvai(cl) = currentPatch%canopy_layer_tvai(cl) + & ! currentCohort%lai + currentCohort%sai do cl = 1,nclmax-1 if(currentCohort%canopy_layer == cl)then - currentPatch%canopy_layer_tai(cl) = currentPatch%canopy_layer_tai(cl) + & + currentPatch%canopy_layer_tvai(cl) = currentPatch%canopy_layer_tvai(cl) + & currentCohort%lai + currentCohort%sai endif enddo @@ -1384,7 +1384,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) if(iv==currentCohort%NV) then remainder = (currentCohort%treelai + currentCohort%treesai) - & - (dinc_ed*dble(currentCohort%NV-1.0_r8)) + (dinc_ed*dble(currentCohort%nv-1.0_r8)) if(remainder > dinc_ed )then write(fates_log(), *)'ED: issue with remainder', & currentCohort%treelai,currentCohort%treesai,dinc_ed, & diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 41cb064f6f..9dd9f46f66 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -21,7 +21,6 @@ module EDCohortDynamicsMod use EDTypesMod , only : min_npm2, min_nppatch use EDTypesMod , only : min_n_safemath use EDTypesMod , only : nlevleaf - use EDTypesMod , only : dinc_ed use FatesInterfaceMod , only : hlm_use_planthydro use FatesPlantHydraulicsMod, only : FuseCohortHydraulics use FatesPlantHydraulicsMod, only : CopyCohortHydraulics @@ -153,7 +152,7 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine call carea_allom(new_cohort%dbh,new_cohort%n,spread,new_cohort%pft,new_cohort%c_area) new_cohort%treelai = tree_lai(new_cohort%bl, new_cohort%status_coh, new_cohort%pft, & - new_cohort%c_area, new_cohort%n, new_cohort%canopy_layer, patchptr%canopy_layer_tai ) + new_cohort%c_area, new_cohort%n, new_cohort%canopy_layer, patchptr%canopy_layer_tvai ) new_cohort%lai = new_cohort%treelai * new_cohort%c_area/patchptr%area new_cohort%treesai = tree_sai(new_cohort%pft, new_cohort%treelai) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index a120835d4c..341c8c9de8 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -1267,7 +1267,7 @@ subroutine zero_patch(cp_p) currentPatch%age = nan currentPatch%age_class = 1 currentPatch%area = nan - currentPatch%canopy_layer_tai(:) = nan + currentPatch%canopy_layer_tvai(:) = nan currentPatch%total_canopy_area = nan currentPatch%tlai_profile(:,:,:) = nan @@ -1347,7 +1347,7 @@ subroutine zero_patch(cp_p) currentPatch%burnt_frac_litter(:) = 0.0_r8 currentPatch%btran_ft(:) = 0.0_r8 - currentPatch%canopy_layer_tai(:) = 0.0_r8 + currentPatch%canopy_layer_tvai(:) = 0.0_r8 currentPatch%seeds_in(:) = 0.0_r8 currentPatch%seed_decay(:) = 0.0_r8 diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 26e8af3fef..3860c872fe 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -20,6 +20,8 @@ module EDPhysiologyMod use EDCohortDynamicsMod , only : create_cohort, sort_cohorts use FatesAllometryMod , only : tree_lai use FatesAllometryMod , only : tree_sai + use FatesAllometryMod , only : decay_coeff_kn + use FatesAllometryMod , only : CumulativeLayerTVAI use EDTypesMod , only : numWaterMem use EDTypesMod , only : dl_sf, dinc_ed @@ -73,7 +75,6 @@ module EDPhysiologyMod private :: seed_decay private :: seed_germination public :: flux_into_litter_pools - public :: decay_coeff_kn logical, parameter :: DEBUG = .false. ! local debug flag @@ -173,33 +174,35 @@ subroutine trim_canopy( currentSite ) type (ed_cohort_type) , pointer :: currentCohort type (ed_patch_type) , pointer :: currentPatch - integer :: z ! leaf layer - integer :: ipft ! pft index - integer :: trimmed ! was this layer trimmed in this year? If not expand the canopy. - real(r8) :: tar_bl ! target leaf biomass (leaves flushed, trimmed) - real(r8) :: tar_bfr ! target fine-root biomass (leaves flushed, trimmed) - real(r8) :: bfr_per_bleaf ! ratio of fine root per leaf biomass - real(r8) :: sla_levleaf ! sla at leaf level z + integer :: z ! leaf layer + integer :: ipft ! pft index + logical :: trimmed ! was this layer trimmed in this year? If not expand the canopy. + real(r8) :: tar_bl ! target leaf biomass (leaves flushed, trimmed) + real(r8) :: tar_bfr ! target fine-root biomass (leaves flushed, trimmed) + real(r8) :: bfr_per_bleaf ! ratio of fine root per leaf biomass + real(r8) :: sla_levleaf ! sla at leaf level z real(r8) :: nscaler_levleaf ! nscaler value at leaf level z - integer :: cl ! canopy layer index - real(r8) :: laican ! canopy sum of lai_z - real(r8) :: vai ! leaf and stem area in this layer - real(r8) :: kn ! nitrogen decay coefficient - real(r8) :: sla_max ! Observational constraint on how large sla (m2/gC) can become - real(r8) :: vai_to_lai ! ratio of vegetation area index (ie. lai+sai) : lai for individual tree + integer :: cl ! canopy layer index + real(r8) :: cum_tvai ! cumulative total vegetation area index at arbitrary leaf layers + real(r8) :: vai ! leaf and stem area in this layer + real(r8) :: kn ! nitrogen decay coefficient + real(r8) :: sla_max ! Observational constraint on how large sla (m2/gC) can become !---------------------------------------------------------------------- + write(fates_log(),*) 'trimming' + currentPatch => currentSite%youngest_patch do while(associated(currentPatch)) + currentCohort => currentPatch%tallest do while (associated(currentCohort)) - trimmed = 0 + trimmed = .false. ipft = currentCohort%pft call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread,currentCohort%pft,currentCohort%c_area) currentCohort%treelai = tree_lai(currentCohort%bl, currentCohort%status_coh, currentCohort%pft, & - currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, currentPatch%canopy_layer_tai ) + currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, currentPatch%canopy_layer_tvai ) currentCohort%treesai = tree_sai(currentCohort%pft, currentCohort%treelai) currentCohort%nv = ceiling((currentCohort%treelai+currentCohort%treesai)/dinc_ed) if (currentCohort%nv > nlevleaf)then @@ -218,45 +221,28 @@ subroutine trim_canopy( currentSite ) ! Identify current canopy layer (cl) cl = currentCohort%canopy_layer - ! Calculate the lai+sai of overlying canopy layers (laican) - if (cl==1) then !are we in the top canopy layer or a shaded layer? - laican = 0._r8 - else - laican = sum(currentPatch%canopy_layer_tai(1:cl-1)) - end if - ! Observational constraint for maximum sla value (m2/gC): - ! m2/gC = m2/gBiomass *kgC/kgBiomass + ! m2/gC = m2/gBiomass * kgBiomass/kgC sla_max = sla_max_drymass * EDPftvarcon_inst%c2b(ipft) - ! Ratio of vegetation area index (ie. lai+sai) to lai for individual tree: - vai_to_lai = 1.0_r8 + (EDPftvarcon_inst%allom_sai_scaler(ipft)/ & - EDPftvarcon_inst%slatop(ipft)) !Leaf cost vs netuptake for each leaf layer. - do z = 1,nlevleaf - - ! Vegetation area index as a function of tree_lai+tree_sai - if(currentCohort%treelai >= (z-1)*dinc_ed)then - ! if tree_lai in this layer is greater than 0 - if (z == 1) then - ! If in first layer, add value equal to vai halfway through the layer - ! Here, vai = layer z's tree_lai + tree_sai = vai_to_lai *dinc_ed - laican = laican + 0.5_r8 * vai_to_lai * dinc_ed - else - ! Since starting from halfway through first layer when z=1, - ! we can add vai for an entire layer - ! to remain at halfway value for subsequent layers. - laican = laican + vai_to_lai * dinc_ed - end if - else ! If tree_lai does not extend into this layer, do not add to laican - end if + do z = 1, currentCohort%nv + + ! Calculate the cumulative total vegetation area index (no snow occlusion, stems and leaves) + cum_tvai = CumulativeLayerTVAI(cl, & + z, & + ipft, & + currentCohort%treelai+currentCohort%treesai, & + currentPatch%canopy_layer_tvai(:)) + + if (currentCohort%year_net_uptake(z) /= 999._r8)then !there was activity this year in this leaf layer. ! Scale for leaf nitrogen profile kn = decay_coeff_kn(ipft) ! Nscaler value at leaf level z - nscaler_levleaf = exp(-kn * laican) + nscaler_levleaf = exp(-kn * cum_tvai) ! Sla value at leaf level z after nitrogen profile scaling (m2/gC) sla_levleaf = EDPftvarcon_inst%slatop(ipft)/nscaler_levleaf @@ -312,7 +298,7 @@ subroutine trim_canopy( currentSite ) currentCohort%laimemory = currentCohort%laimemory * & (1.0_r8 - EDPftvarcon_inst%trim_inc(ipft)) endif - trimmed = 1 + trimmed = .true. endif endif endif @@ -320,7 +306,7 @@ subroutine trim_canopy( currentSite ) enddo !z currentCohort%year_net_uptake(:) = 999.0_r8 - if (trimmed == 0.and.currentCohort%canopy_trim < 1.0_r8)then + if ( (.not.trimmed) .and.currentCohort%canopy_trim < 1.0_r8)then currentCohort%canopy_trim = currentCohort%canopy_trim + EDPftvarcon_inst%trim_inc(ipft) endif @@ -2472,27 +2458,9 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) end subroutine flux_into_litter_pools - - real(r8) function decay_coeff_kn(pft) - - ! ============================================================================ - ! Decay coefficient (kn) is a function of vcmax25top for each pft. - ! ============================================================================ - - !ARGUMENTS - integer, intent(in) :: pft - - !LOCAL VARIABLES - ! ----------------------------------------------------------------------------------- + ! =================================================================================== - ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 used - ! kn = 0.11. Here, we derive kn from vcmax25 as in Lloyd et al - ! (2010) Biogeosciences, 7, 1833-1859 - - decay_coeff_kn = exp(0.00963_r8 * EDPftvarcon_inst%vcmax25top(pft) - 2.43_r8) - return - end function decay_coeff_kn end module EDPhysiologyMod diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index e5c33a570f..6f7f499f23 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -113,6 +113,8 @@ module FatesAllometryMod public :: bdead_allom ! Generic bdead wrapper public :: carea_allom ! Generic crown area wrapper public :: bstore_allom ! Generic maximum storage carbon wrapper + public :: CumulativeLayerTVAI + public :: decay_coeff_kn public :: StructureResetOfDH ! Method to set DBH to sync with structure biomass public :: CheckIntegratedAllometries public :: set_root_fraction ! Generic wrapper to calculate normalized @@ -136,7 +138,12 @@ module FatesAllometryMod integer, parameter, public :: i_hydro_rootprof_context = 1 integer, parameter, public :: i_biomass_rootprof_context = 2 - real, parameter, public :: sla_max_drymass = 0.0477_r8 ! Max sla (m2/g dry mass) obs for all plants (Kattge et al. 2011) + + + ! Max sla (m2/g dry mass) obs for all plants (Kattge et al. 2011) + + real, parameter, public :: sla_max_drymass = 0.0477_r8 + ! If testing b4b with older versions, do not remove sapwood ! Our old methods with saldarriaga did not remove sapwood from the @@ -543,29 +550,37 @@ end subroutine storage_fraction_of_target ! ===================================================================================== - real(r8) function tree_lai( bl, status_coh, pft, c_area, n, cl, canopy_layer_tai ) + real(r8) function tree_lai( bl, status_coh, pft, c_area, n, cl, canopy_layer_tvai ) - ! ============================================================================ - ! LAI of individual trees is a function of the total leaf area and the total canopy area. - ! ============================================================================ - ! !ARGUMENTS - real(r8), intent(in) :: bl ! plant leaf biomass [kg] - integer, intent(in) :: status_coh ! growth status of plant (2 = leaves on , 1 = leaves off) + ! ----------------------------------------------------------------------------------- + ! LAI of individual trees is a function of the total leaf area and the total + ! canopy area. + ! ---------------------------------------------------------------------------------- + + ! !ARGUMENTS + real(r8), intent(in) :: bl ! plant leaf biomass [kg] + integer, intent(in) :: status_coh ! growth status of plant + ! (2 = leaves on , 1 = leaves off) integer, intent(in) :: pft - real(r8), intent(in) :: c_area ! areal extent of canopy (m2) - real(r8), intent(in) :: n ! number of individuals in cohort per 'area'(10000m2 default) - integer, intent(in) :: cl ! canopy layer index - real(r8), intent(in) :: canopy_layer_tai(nclmax) ! total area index of each canopy layer + real(r8), intent(in) :: c_area ! areal extent of canopy (m2) + real(r8), intent(in) :: n ! number of individuals in cohort per ha + integer, intent(in) :: cl ! canopy layer index + real(r8), intent(in) :: canopy_layer_tvai(nclmax) ! total vegetated area index of + ! each canopy layer - ! !LOCAL VARIABLES: + ! !LOCAL VARIABLES: real(r8) :: leafc_per_unitarea ! KgC of leaf per m2 area of ground. - real(r8) :: slat ! the sla of the top leaf layer. m2/kgC - real(r8) :: laican ! lai + sai of canopy layer overlying this tree - real(r8) :: vai_to_lai ! ratio of vegetation area index (ie. sai+lai) to lai for individual tree - real(r8) :: kn ! coefficient for exponential decay of 1/sla and vcmax with canopy depth - real(r8) :: sla_max ! Observational constraint on how large sla (m2/gC) can become - real(r8) :: leafc_slamax ! Leafc_per_unitarea at which sla_max is reached - real(r8) :: clim ! Upper limit for leafc_per_unitarea in exponential tree_lai function + real(r8) :: slat ! the sla of the top leaf layer. m2/kgC + real(r8) :: can_tvai ! total lai + sai (vai) of canopy layer overlying this tree + real(r8) :: vai_to_lai ! ratio of vegetation area index (ie. sai+lai) + ! to lai for individual tree + real(r8) :: kn ! coefficient for exponential decay of 1/sla and + ! vcmax with canopy depth + real(r8) :: sla_max ! Observational constraint on how large sla + ! (m2/gC) can become + real(r8) :: leafc_slamax ! Leafc_per_unitarea at which sla_max is reached + real(r8) :: clim ! Upper limit for leafc_per_unitarea in exponential + ! tree_lai function !---------------------------------------------------------------------- if( bl < 0._r8 .or. pft == 0 ) then @@ -577,24 +592,25 @@ real(r8) function tree_lai( bl, status_coh, pft, c_area, n, cl, canopy_layer_tai if(leafc_per_unitarea > 0.0_r8)then - ! Calculate laican = LAI + SAI of overlying canopy layer + ! Calculate can_tvai = LAI + SAI of overlying canopy layer if (cl==1) then ! if in we are in the canopy (top) layer) - laican = 0._r8 + can_tvai = 0._r8 else - laican = sum(canopy_layer_tai(1:cl-1)) + can_tvai = sum(canopy_layer_tvai(1:cl-1)) end if ! Ratio of vegetation area index (ie. lai+sai) to lai for individual tree: vai_to_lai = 1.0_r8 + (EDPftvarcon_inst%allom_sai_scaler(pft)/ & EDPftvarcon_inst%slatop(pft)) + ! Coefficient for exponential decay of 1/sla with canopy depth: - kn = exp(0.00963_r8 * EDPftvarcon_inst%vcmax25top(pft) - 2.43_r8) + kn = decay_coeff_kn(pft) ! Observational constraint for maximum sla value (m2/kgC): - ! m2/kgC = g/kg* m2/gBiomass *kgC/kgBiomass + ! m2/kgC = g/kg* m2/gBiomass * kgBiomass/kgC sla_max = g_per_kg * sla_max_drymass * EDPftvarcon_inst%c2b(pft) ! Leafc_per_unitarea at which sla_max is reached due to exponential sla profile in canopy: - leafc_slamax = (slat - sla_max * exp(-1.0_r8 * kn * laican)) / & + leafc_slamax = (slat - sla_max * exp(-1.0_r8 * kn * can_tvai)) / & (-1.0_r8 * kn * vai_to_lai * slat * sla_max) if(leafc_slamax < 0.0_r8)then leafc_slamax = 0.0_r8 @@ -606,23 +622,23 @@ real(r8) function tree_lai( bl, status_coh, pft, c_area, n, cl, canopy_layer_tai ! sla with depth in the canopy will not exceed sla_max. ! In this case, we can use an exponential profile for sla throughout the entire canopy. ! The exponential profile for sla is given by: - ! sla(at a given canopy depth) = slat / exp(-kn (laican + vai_to_lai * tree_lai) + ! sla(at a given canopy depth) = slat / exp(-kn (can_tvai + vai_to_lai * tree_lai) ! where vai_to_lai * tree_lai = tree_lai + tree_sai ! We can solve for tree_lai using the above function for the sla profile and first setting - ! leafc_per_unitarea = integral of e^(-kn(vai_to_lai * x + laican)) / slatop + ! leafc_per_unitarea = integral of e^(-kn(vai_to_lai * x + can_tvai)) / slatop ! over x = 0 to tree_lai ! Then, rearranging the equation to solve for tree_lai. if (leafc_per_unitarea <= leafc_slamax)then - tree_lai = (log(exp(-1.0_r8 * kn * laican) - & + tree_lai = (log(exp(-1.0_r8 * kn * can_tvai) - & kn * vai_to_lai * slat * leafc_per_unitarea) + & - (kn * laican)) / (-1.0_r8 * kn * vai_to_lai) + (kn * can_tvai)) / (-1.0_r8 * kn * vai_to_lai) ! If leafc_per_unitarea becomes too large, tree_lai becomes an imaginary number ! (because the tree_lai equation requires us to take the natural log of something >0) ! Thus, we include the following error message in case leafc_per_unitarea becomes too large. - clim = (exp(-1.0_r8 * kn * laican)) / (kn * vai_to_lai * slat) + clim = (exp(-1.0_r8 * kn * can_tvai)) / (kn * vai_to_lai * slat) if(leafc_per_unitarea >= clim)then - write(fates_log(),*) 'too much leafc_per_unitarea' , leafc_per_unitarea, clim, pft, laican + write(fates_log(),*) 'too much leafc_per_unitarea' , leafc_per_unitarea, clim, pft, can_tvai write(fates_log(),*) 'Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) endif @@ -636,17 +652,18 @@ real(r8) function tree_lai( bl, status_coh, pft, c_area, n, cl, canopy_layer_tai ! Add exponential and linear portions of tree_lai ! Exponential term for leafc = leafc_slamax; ! Linear term (static sla = sla_max) for portion of leafc > leafc_slamax - tree_lai = ((log(exp(-1.0_r8 * kn * laican) - & + tree_lai = ((log(exp(-1.0_r8 * kn * can_tvai) - & kn * vai_to_lai * slat * leafc_slamax) + & - (kn * laican)) / (-1.0_r8 * kn * vai_to_lai)) + & + (kn * can_tvai)) / (-1.0_r8 * kn * vai_to_lai)) + & (leafc_per_unitarea - leafc_slamax) * sla_max ! if leafc_slamax becomes too large, tree_lai_exp becomes an imaginary number ! (because the tree_lai equation requires us to take the natural log of something >0) ! Thus, we include the following error message in case leafc_slamax becomes too large. - clim = (exp(-1.0_r8 * kn * laican)) / (kn * vai_to_lai * slat) + clim = (exp(-1.0_r8 * kn * can_tvai)) / (kn * vai_to_lai * slat) if(leafc_slamax >= clim)then - write(fates_log(),*) 'too much leafc_slamax' , leafc_per_unitarea, leafc_slamax, clim, pft, laican + write(fates_log(),*) 'too much leafc_slamax' , & + leafc_per_unitarea, leafc_slamax, clim, pft, can_tvai write(fates_log(),*) 'Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) endif @@ -2051,6 +2068,98 @@ subroutine jackson_beta_root_profile(root_fraction, ft, zi) return end subroutine jackson_beta_root_profile + ! ===================================================================================== + + + real(r8) function decay_coeff_kn(pft) + + ! --------------------------------------------------------------------------------- + ! This function estimates the decay coefficient used to estimate vertical + ! attenuation of properties in the canopy. + ! + ! Decay coefficient (kn) is a function of vcmax25top for each pft. + ! + ! Currently, this decay is applied to vcmax attenuation, and SLA (optionally) + ! + ! --------------------------------------------------------------------------------- + + !ARGUMENTS + integer, intent(in) :: pft + + !LOCAL VARIABLES + ! ----------------------------------------------------------------------------------- + + ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 used + ! kn = 0.11. Here, we derive kn from vcmax25 as in Lloyd et al + ! (2010) Biogeosciences, 7, 1833-1859 + + decay_coeff_kn = exp(0.00963_r8 * EDPftvarcon_inst%vcmax25top(pft) - 2.43_r8) + + return + end function decay_coeff_kn + + ! ===================================================================================== + + function CumulativeLayerTVAI(icanlayer, & + ileaflayer, & + ipft, & + current_tvai, & + canopy_layer_tvai) result(cum_tvai) + + ! ----------------------------------------------------------------------------------- + ! This function calculates the cumulative (top-down) vegetation area index + ! for a given leaf layer in a canopy layer. + ! + ! A top canopy layer (icanlayer==1) has now leaf+stem area above it. + ! + ! Important note: This subroutine operates on total area indices, and not effective. + ! This is done to promote consistency with calculations of tree_tai + ! which is an optional method of calculating the SLA decay in + ! the profile. + ! ----------------------------------------------------------------------------------- + + integer, intent(in) :: icanlayer ! Layer index for the current canopy + integer, intent(in) :: ileaflayer ! Layer index for the current leaf layer + integer, intent(in) :: ipft ! PFT index + real(r8), intent(in) :: current_tvai ! This is the total vegetation area index + ! for the current canopy layer, for the + ! entity of interest. Note + ! depending on where this routine + ! is called, this entity may be for the current + ! cohort, or for all pfts in this layer + real(r8), intent(in) :: canopy_layer_tvai(:) ! The total vegetation index of + ! each canopy layer + + real(r8) :: cum_tvai ! Resulting cumulative vegetation + ! area index + + real(r8) :: tvai0 ! tvai of leaf-layers up to current + real(r8) :: tvai ! tvai of current layer + + + ! Calculate the tvai of canopy layers above the current + + if (icanlayer==1) then + cum_tvai = 0._r8 + else + cum_tvai = sum(canopy_layer_tvai(1:icanlayer-1)) + end if + + tvai0 = dble(ileaflayer-1)*dinc_ed + + tvai = min(dinc_ed, current_tvai-tvai0) + + if (tvai<0.0_r8) then + write(fates_log(),*) 'A leaf layer greater than max for this cohort was specified' + write(fates_log(),*) tvai,tvai0,current_tvai,cum_tvai,icanlayer,ileaflayer + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + ! Apply the cumulative TVAI at the current layer mid-point (ie 0.5*) + cum_tvai = cum_tvai + tvai0 + 0.5_r8*tvai + + return + end function CumulativeLayerTVAI ! ===================================================================================== diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 20f1021b82..2d6efd1eac 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -29,7 +29,6 @@ module FATESPlantRespPhotosynthMod use EDTypesMod, only : maxpft use EDTypesMod, only : nlevleaf use EDTypesMod, only : nclmax - use EDTypesMod, only : dinc_ed ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg @@ -85,6 +84,8 @@ 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 : i_hydro_rootprof_context + use FatesAllometryMod, only : CumulativeLayerTVAI + use FatesAllometryMod, only : decay_coeff_kn ! ARGUMENTS: ! ----------------------------------------------------------------------------------- @@ -149,7 +150,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) real(r8) :: ceair ! vapor pressure of air, constrained (Pa) real(r8) :: nscaler ! leaf nitrogen scaling coefficient real(r8) :: leaf_frac ! ratio of to leaf biomass to total alive biomass - real(r8) :: laican ! canopy sum of lai_z + real(r8) :: cum_tvai ! cumulative total vegetation area index of descrete leaf layers real(r8) :: tcsoi ! Temperature response function for root respiration. real(r8) :: tcwood ! Temperature response function for wood @@ -295,13 +296,9 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 used ! kn = 0.11. Here, derive kn from vcmax25 as in Lloyd et al ! (2010) Biogeosciences, 7, 1833-1859 - ! Remove daylength factor from vcmax25 so that kn is based on maximum vcmax25 - if (bc_in(s)%dayl_factor_pa(ifp) == 0._r8) then - kn(ft) = 0._r8 - else - kn(ft) = exp(0.00963_r8 * EDPftvarcon_inst%vcmax25top(ft) - 2.43_r8) - end if + kn(ft) = decay_coeff_kn(ft) + ! This is probably unnecessary and already calculated ! ALSO, THIS ROOTING PROFILE IS USED TO CALCULATE RESPIRATION @@ -337,6 +334,8 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! ------------------------------------------------------------------------ rate_mask_z(:,1:numpft,:) = .false. + write(fates_log(),*) 'photo' + if(currentPatch%countcohorts > 0.0)then ! Ignore empty patches currentCohort => currentPatch%tallest @@ -355,14 +354,6 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! are there any leaves of this pft in this layer? if(currentPatch%canopy_mask(cl,ft) == 1)then - if(cl==1)then !are we in the top canopy layer or a shaded layer? - laican = 0._r8 - else - - laican = sum(currentPatch%canopy_layer_tai(1:cl-1)) - - end if - ! Loop over leaf-layers do iv = 1,currentCohort%nv @@ -392,33 +383,17 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) bbb = max (bbbopt(nint(c3psn(ft)))*currentPatch%btran_ft(ft), 1._r8) btran_eff = currentPatch%btran_ft(ft) end if + - ! Ratio of vegetation area index (ie. lai + sai) to lai: - vai_to_lai = 1.0_r8 + (EDPftvarcon_inst%allom_sai_scaler(ft)/ & - EDPftvarcon_inst%slatop(ft)) + cum_tvai = CumulativeLayerTVAI(cl, & + iv, & + ft, & + currentCohort%treelai+currentCohort%treesai, & + currentPatch%canopy_layer_tvai(:)) - ! Add this leaf layer's vegetation area index to laican - ! For consistency between the vcmax profile (based on elai+esai) - ! and the sla profile (based on tree_lai+tree_sai), - ! we add the vai value halfway through the layer. - if(currentPatch%elai_profile(cl,ft,iv) > 0.0_r8)then - ! if elai in this layer is greater than 0 - if (iv == 1) then - ! If in first layer, - ! add value equal to exposed vai halfway through the layer - ! here a full layer's exposed vai = ratio of lai+sai : lai * dinc_ed - laican = laican + 0.5_r8 * vai_to_lai * dinc_ed - else - ! Since starting from halfway through first layer when iv=1, - ! We can add vai for an entire layer - ! to remain at halway value for subsequent layers - laican = laican + vai_to_lai * dinc_ed - end if - else ! if elai does not extend into this layer, do not add to laican - end if ! Scale for leaf nitrogen profile - nscaler = exp(-kn(ft) * laican) + nscaler = exp(-kn(ft) * cum_tvai) ! Part VII: Calculate dark respiration (leaf maintenance) for this layer call LeafLayerMaintenanceRespiration( param_derived%lmr25top(ft),& ! in @@ -1122,7 +1097,6 @@ subroutine ScaleLeafLayerFluxToCohort(nv, & ! in currentCohort%nv ! ------------------------------------------------------------------------------------ use FatesConstantsMod, only : umolC_to_kgC - use EDTypesMod, only : dinc_ed ! Arguments integer, intent(in) :: nv ! number of active leaf layers diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 4aa90757a0..6ad2765cd7 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -81,7 +81,8 @@ module EDTypesMod ! BIOLOGY/BIOGEOCHEMISTRY integer , parameter :: external_recruitment = 0 ! external recruitment flag 1=yes integer , parameter :: SENES = 10 ! Window of time over which we track temp for cold sensecence (days) - real(r8), parameter :: DINC_ED = 1.0_r8 ! size of LAI bins. + real(r8), parameter :: dinc_ed = 1.0_r8 ! size of VAI bins (LAI+SAI) [CHANGE THIS NAME WITH NEXT INTERFACE + ! UPDATE] integer , parameter :: N_DIST_TYPES = 3 ! Disturbance Modes 1) tree-fall, 2) fire, 3) logging integer , parameter :: dtype_ifall = 1 ! index for naturally occuring tree-fall generated event integer , parameter :: dtype_ifire = 2 ! index for fire generated disturbance event @@ -312,7 +313,7 @@ module EDTypesMod ! LEAF ORGANIZATION real(r8) :: pft_agb_profile(maxpft,n_dbh_bins) ! binned above ground biomass, for patch fusion: KgC/m2 - real(r8) :: canopy_layer_tai(nclmax) ! total area index of each canopy layer + real(r8) :: canopy_layer_tvai(nclmax) ! total vegetation area index of each canopy layer ! used to determine attenuation of parameters during ! photosynthesis m2 veg / m2 of canopy area (patch without bare ground) real(r8) :: total_canopy_area ! area that is covered by vegetation : m2