From d364f51258098f7729ead32efe4caebbf0b83b4c Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 10 Jul 2018 16:18:01 -0700 Subject: [PATCH 1/5] Shortened some lines. Changed a trimmed flag to use a logical --- biogeochem/EDPhysiologyMod.F90 | 41 ++++++++++++------------ biogeochem/FatesAllometryMod.F90 | 53 ++++++++++++++++++++------------ 2 files changed, 55 insertions(+), 39 deletions(-) diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 26e8af3fef..6c3050bb61 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -173,29 +173,30 @@ 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) :: 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 !---------------------------------------------------------------------- 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, & @@ -312,7 +313,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 +321,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,12 +2473,14 @@ 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 diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index e5c33a570f..3c8c034f37 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -136,7 +136,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 @@ -545,27 +550,34 @@ end subroutine storage_fraction_of_target real(r8) function tree_lai( bl, status_coh, pft, c_area, n, cl, canopy_layer_tai ) - ! ============================================================================ - ! 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_tai(nclmax) ! total 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) :: 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 !---------------------------------------------------------------------- if( bl < 0._r8 .or. pft == 0 ) then @@ -646,7 +658,8 @@ real(r8) function tree_lai( bl, status_coh, pft, c_area, n, cl, canopy_layer_tai ! 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) 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, laican write(fates_log(),*) 'Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) endif From 3eab20d900e90eafce6b2d390b72d854444ea42c Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 11 Jul 2018 16:09:58 -0700 Subject: [PATCH 2/5] converted the cumulative tvai scheme to a function call. changed some variable names to be more accurate. made all calculations of the decay kn to use the function --- biogeochem/EDCanopyStructureMod.F90 | 16 ++-- biogeochem/EDCohortDynamicsMod.F90 | 3 +- biogeochem/EDPatchDynamicsMod.F90 | 4 +- biogeochem/EDPhysiologyMod.F90 | 71 ++++------------ biogeochem/FatesAllometryMod.F90 | 96 +++++++++++++++++++++- biogeophys/FatesPlantRespPhotosynthMod.F90 | 54 ++++-------- main/EDTypesMod.F90 | 5 +- 7 files changed, 141 insertions(+), 108 deletions(-) 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 6c3050bb61..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 @@ -182,14 +183,15 @@ subroutine trim_canopy( currentSite ) 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) :: 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 - real(r8) :: vai_to_lai ! ratio of vegetation area index (ie. lai+sai) : lai for individual tree !---------------------------------------------------------------------- + write(fates_log(),*) 'trimming' + currentPatch => currentSite%youngest_patch do while(associated(currentPatch)) @@ -200,7 +202,7 @@ subroutine trim_canopy( currentSite ) 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 @@ -219,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 @@ -2476,26 +2461,6 @@ 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 3c8c034f37..6c159a4716 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 @@ -599,8 +601,9 @@ real(r8) function tree_lai( bl, status_coh, pft, c_area, n, cl, canopy_layer_tai ! 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 @@ -2064,6 +2067,97 @@ 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 + + cum_tvai = cum_tvai + tvai0 + 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 From c0dbd7d7a85b2191a5766eeee900525347fc25c7 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 11 Jul 2018 16:21:27 -0700 Subject: [PATCH 3/5] Updated some variable names and descriptions in tree_lai --- biogeochem/FatesAllometryMod.F90 | 49 ++++++++++++++++---------------- 1 file changed, 25 insertions(+), 24 deletions(-) diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 6c159a4716..1293425612 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -550,7 +550,7 @@ 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 @@ -558,19 +558,20 @@ real(r8) function tree_lai( bl, status_coh, pft, c_area, n, cl, canopy_layer_tai ! ---------------------------------------------------------------------------------- ! !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) + 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 ha - 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: 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) :: 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 @@ -591,11 +592,11 @@ 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: @@ -606,10 +607,10 @@ real(r8) function tree_lai( bl, status_coh, pft, c_area, n, cl, canopy_layer_tai 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 @@ -621,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 @@ -651,18 +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 + leafc_per_unitarea, leafc_slamax, clim, pft, can_tvai write(fates_log(),*) 'Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) endif From 516670a48901251bfcc968e20f7eed0cd6267c1a Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 11 Jul 2018 16:26:25 -0700 Subject: [PATCH 4/5] Updated cumulative TVAI function to estimate the mid-point of layer of interest --- biogeochem/FatesAllometryMod.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 1293425612..8ced4209dc 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -2155,7 +2155,8 @@ function CumulativeLayerTVAI(icanlayer, & call endrun(msg=errMsg(sourcefile, __LINE__)) end if - cum_tvai = cum_tvai + tvai0 + tvai + ! Apply the cumulative TVAI at the current layer mid-point (ie 0.5*) + cum_tvai = cum_tvai + tvai0 + 0.5*tvai return end function CumulativeLayerTVAI From ac9085dd4f08dd9b13121d814ea065ac65a8149f Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 11 Jul 2018 16:27:26 -0700 Subject: [PATCH 5/5] Updated hard-constant to be an r8 --- biogeochem/FatesAllometryMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 8ced4209dc..6f7f499f23 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -2156,7 +2156,7 @@ function CumulativeLayerTVAI(icanlayer, & end if ! Apply the cumulative TVAI at the current layer mid-point (ie 0.5*) - cum_tvai = cum_tvai + tvai0 + 0.5*tvai + cum_tvai = cum_tvai + tvai0 + 0.5_r8*tvai return end function CumulativeLayerTVAI