From cc4e819380e50110c4361f5caf890d70f12ed8df Mon Sep 17 00:00:00 2001 From: Charlie Koven Date: Mon, 16 Jul 2018 11:46:23 -0700 Subject: [PATCH 1/2] revert use_slaprofile switch logic, and replace with PFT-level control via the slamax trait --- biogeochem/EDPhysiologyMod.F90 | 37 ++---- biogeochem/FatesAllometryMod.F90 | 157 ++++++++++------------- main/EDPftvarcon.F90 | 12 +- parameter_files/fates_params_14pfts.cdl | 6 + parameter_files/fates_params_default.cdl | 5 + 5 files changed, 105 insertions(+), 112 deletions(-) diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 67c7cb327d..8174e34c18 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -164,8 +164,6 @@ subroutine trim_canopy( currentSite ) ! Canopy trimming / leaf optimisation. Removes leaves in negative annual carbon balance. ! ! !USES: - use FatesAllometryMod , only : sla_max_drymass, use_slaprofile - use FatesConstantsMod , only : g_per_kg ! !ARGUMENTS type (ed_site_type),intent(inout), target :: currentSite @@ -221,9 +219,8 @@ subroutine trim_canopy( currentSite ) ! Identify current canopy layer (cl) cl = currentCohort%canopy_layer - ! Observational constraint for maximum sla value (m2/gC): - ! m2/gC = m2/gBiomass * kgBiomass/kgC - sla_max = sla_max_drymass * EDPftvarcon_inst%c2b(ipft) + ! PFT-level maximum SLA value, even if under a thick canopy (same units as slatop) + sla_max = EDPftvarcon_inst%slamax(ipft) !Leaf cost vs netuptake for each leaf layer. do z = 1, currentCohort%nv @@ -239,27 +236,19 @@ subroutine trim_canopy( currentSite ) if (currentCohort%year_net_uptake(z) /= 999._r8)then !there was activity this year in this leaf layer. - if(use_slaprofile == 0)then - - ! Use slatop for sla value at all leaf levels - sla_levleaf = EDPftvarcon_inst%slatop(ipft) - - elseif(use_slaprofile == 1)then - - ! Calculate sla_levleaf following the sla profile with overlying leaf area - ! Scale for leaf nitrogen profile - kn = decay_coeff_kn(ipft) - ! Nscaler value at leaf level z - 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 - - if(sla_levleaf > sla_max)then - sla_levleaf = sla_max - end if + ! Calculate sla_levleaf following the sla profile with overlying leaf area + ! Scale for leaf nitrogen profile + kn = decay_coeff_kn(ipft) + ! Nscaler value at leaf level z + 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 + + if(sla_levleaf > sla_max)then + sla_levleaf = sla_max end if - + !Leaf Cost kgC/m2/year-1 !decidous costs. if (EDPftvarcon_inst%season_decid(ipft) == 1.or. & diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 23c96e8152..2d19e7987a 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -140,16 +140,6 @@ module FatesAllometryMod integer, parameter, public :: i_biomass_rootprof_context = 2 - ! Use sla profile (ie. increase sla with overlying tvai)? - ! Value of 1 = use sla profile, 0 = do NOT use sla profile - - integer, parameter, public :: use_slaprofile = 1 - - ! 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 ! bdead pool. But newer allometries are providing total agb @@ -596,90 +586,83 @@ real(r8) function tree_lai( bl, status_coh, pft, c_area, n, cl, canopy_layer_tva leafc_per_unitarea = bl/(c_area/n) !KgC/m2 if(leafc_per_unitarea > 0.0_r8)then - - if(use_slaprofile == 0) then - - tree_lai = leafc_per_unitarea * slat !kg/m2 * m2/kg = unitless LAI - - else if(use_slaprofile == 1) then - - ! Calculate can_tvai = LAI + SAI of overlying canopy layer - if (cl==1) then ! if in we are in the canopy (top) layer) - can_tvai = 0._r8 - else - 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 = decay_coeff_kn(pft) - - ! Observational constraint for maximum sla value (m2/kgC): - ! 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 * can_tvai)) / & - (-1.0_r8 * kn * vai_to_lai * slat * sla_max) - if(leafc_slamax < 0.0_r8)then - leafc_slamax = 0.0_r8 + ! Calculate can_tvai = LAI + SAI of overlying canopy layer + if (cl==1) then ! if in we are in the canopy (top) layer) + can_tvai = 0._r8 + else + 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 = decay_coeff_kn(pft) + + ! take PFT-level maximum SLA value, even if under a thick canopy (which has units of m2/gC), + ! and put into units of m2/kgC + sla_max = g_per_kg *EDPftvarcon_inst%slamax(ipft) + ! 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 * can_tvai)) / & + (-1.0_r8 * kn * vai_to_lai * slat * sla_max) + if(leafc_slamax < 0.0_r8)then + leafc_slamax = 0.0_r8 + endif + + ! Calculate tree_lai (m2 leaf area /m2 ground) = unitless LAI + !---------------------------------------------------------------------- + ! If leafc_per_unitarea is less than leafc_slamax, + ! 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 (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 + 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 * can_tvai) - & + kn * vai_to_lai * slat * leafc_per_unitarea) + & + (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 * 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, can_tvai + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) endif - - ! Calculate tree_lai (m2 leaf area /m2 ground) = unitless LAI - !---------------------------------------------------------------------- - ! If leafc_per_unitarea is less than leafc_slamax, - ! 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 (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 + 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 * can_tvai) - & - kn * vai_to_lai * slat * leafc_per_unitarea) + & - (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 * 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, can_tvai - write(fates_log(),*) 'Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) - endif - + ! When leafc_per_unitarea is greater than leafc_slamax, ! tree_lai could become so great that the sla profile surpasses sla_max at depth. ! In this case, we use the exponential profile to calculate tree_lai until ! we reach the maximum allowed value for sla (sla_max). ! Then, calculate the remaining tree_lai using a linear function of sla_max and the remaining leafc. - else if(leafc_per_unitarea > leafc_slamax)then - ! 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 * can_tvai) - & - kn * vai_to_lai * slat * leafc_slamax) + & - (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 * 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, can_tvai - write(fates_log(),*) 'Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) - endif - end if ! (leafc_per_unitarea > leafc_slamax) - end if ! (use_slaprofile) + else if(leafc_per_unitarea > leafc_slamax)then + ! 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 * can_tvai) - & + kn * vai_to_lai * slat * leafc_slamax) + & + (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 * 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, can_tvai + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + endif + end if ! (leafc_per_unitarea > leafc_slamax) else tree_lai = 0.0_r8 endif ! (leafc_per_unitarea > 0.0_r8) diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index b72cf67c30..a5c70ca2a4 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -51,6 +51,7 @@ module EDPftvarcon real(r8), allocatable :: stress_decid(:) real(r8), allocatable :: season_decid(:) real(r8), allocatable :: evergreen(:) + real(r8), allocatable :: slamax(:) real(r8), allocatable :: slatop(:) real(r8), allocatable :: leaf_long(:) real(r8), allocatable :: roota_par(:) @@ -358,6 +359,10 @@ subroutine Register_PFT(this, fates_params) call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_leaf_slamax' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_leaf_slatop' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -787,6 +792,10 @@ subroutine Receive_PFT(this, fates_params) call fates_params%RetreiveParameterAllocate(name=name, & data=this%evergreen) + name = 'fates_leaf_slamax' + call fates_params%RetreiveParameterAllocate(name=name, & + data=this%slamax) + name = 'fates_leaf_slatop' call fates_params%RetreiveParameterAllocate(name=name, & data=this%slatop) @@ -1496,7 +1505,8 @@ subroutine FatesReportPFTParams(is_master) write(fates_log(),fmt0) 'stress_decid = ',EDPftvarcon_inst%stress_decid write(fates_log(),fmt0) 'season_decid = ',EDPftvarcon_inst%season_decid write(fates_log(),fmt0) 'evergreen = ',EDPftvarcon_inst%evergreen - write(fates_log(),fmt0) 'slatop = ',EDPftvarcon_inst%slatop + write(fates_log(),fmt0) 'slamax = ',EDPftvarcon_inst%slamax + write(fates_log(),fmt0) 'slatop = ',EDPftvarcon_inst%slatop write(fates_log(),fmt0) 'leaf_long = ',EDPftvarcon_inst%leaf_long write(fates_log(),fmt0) 'roota_par = ',EDPftvarcon_inst%roota_par write(fates_log(),fmt0) 'rootb_par = ',EDPftvarcon_inst%rootb_par diff --git a/parameter_files/fates_params_14pfts.cdl b/parameter_files/fates_params_14pfts.cdl index a94b3e747b..d8404a9883 100644 --- a/parameter_files/fates_params_14pfts.cdl +++ b/parameter_files/fates_params_14pfts.cdl @@ -308,6 +308,9 @@ variables: float fates_leaf_long(fates_pft) ; fates_leaf_long:units = "yr" ; fates_leaf_long:long_name = "Leaf longevity (ie turnover timescale)" ; + float fates_leaf_slamax(fates_pft) ; + fates_leaf_slamax:units = "m^2/gC" ; + fates_leaf_slamax:long_name = "Maximum Specific Leaf Area (SLA), even if under a dense canopy" ; float fates_leaf_slatop(fates_pft) ; fates_leaf_slatop:units = "m^2/gC" ; fates_leaf_slatop:long_name = "Specific Leaf Area (SLA) at top of canopy, projected area basis" ; @@ -870,6 +873,9 @@ data: fates_leaf_long = 1.5, 4, 6, 1, 1.5, 1, 1, 1, 1.5, 1, 1, 1, 1, 1 ; + fates_leaf_slamax = 0.0954, 0.0954, 0.0954, 0.0954, 0.0954, 0.0954, 0.0954, 0.0954, + 0.012, 0.03, 0.03, 0.03, 0.03, 0.03 ; + fates_leaf_slatop = 0.012, 0.01, 0.008, 0.024, 0.012, 0.03, 0.03, 0.03, 0.012, 0.03, 0.03, 0.03, 0.03, 0.03 ; diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index 6bb903afd3..1a619b7ce8 100644 --- a/parameter_files/fates_params_default.cdl +++ b/parameter_files/fates_params_default.cdl @@ -248,6 +248,9 @@ variables: float fates_leaf_long(fates_pft) ; fates_leaf_long:units = "yr" ; fates_leaf_long:long_name = "Leaf longevity (ie turnover timescale)" ; + float fates_leaf_slamax(fates_pft) ; + fates_leaf_slamax:units = "m^2/gC" ; + fates_leaf_slamax:long_name = "Maximum Specific Leaf Area (SLA), even if under a dense canopy" ; float fates_leaf_slatop(fates_pft) ; fates_leaf_slatop:units = "m^2/gC" ; fates_leaf_slatop:long_name = "Specific Leaf Area (SLA) at top of canopy, projected area basis" ; @@ -758,6 +761,8 @@ data: fates_leaf_long = 1.5, 1.5 ; + fates_leaf_slamax = 0.0954, 0.0954 ; + fates_leaf_slatop = 0.012, 0.012 ; fates_leaf_stor_priority = 0.8, 0.8 ; From 1fa12faa2f2724f96b4555abfb67d06b2b5ac682 Mon Sep 17 00:00:00 2001 From: ckoven Date: Mon, 16 Jul 2018 17:59:30 -0600 Subject: [PATCH 2/2] bugfix on prior --- biogeochem/FatesAllometryMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 2d19e7987a..028b4fa5c4 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -603,7 +603,7 @@ real(r8) function tree_lai( bl, status_coh, pft, c_area, n, cl, canopy_layer_tva ! take PFT-level maximum SLA value, even if under a thick canopy (which has units of m2/gC), ! and put into units of m2/kgC - sla_max = g_per_kg *EDPftvarcon_inst%slamax(ipft) + sla_max = g_per_kg *EDPftvarcon_inst%slamax(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 * can_tvai)) / & (-1.0_r8 * kn * vai_to_lai * slat * sla_max)