Skip to content

Commit

Permalink
Merge pull request #3 from ckoven/kovenock_slaprofile
Browse files Browse the repository at this point in the history 
make slamax a PFT profile
  • Loading branch information
@kovenock
kovenock authored Jul 17, 2018
2 parents 08ab6c6 + 1fa12fa commit 770ac36
Show file tree
Hide file tree
Showing 5 changed files with 105 additions and 112 deletions.
37 changes: 13 additions & 24 deletions biogeochem/EDPhysiologyMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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. &
Expand Down
157 changes: 70 additions & 87 deletions biogeochem/FatesAllometryMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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(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
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)
Expand Down
12 changes: 11 additions & 1 deletion main/EDPftvarcon.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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(:)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
6 changes: 6 additions & 0 deletions parameter_files/fates_params_14pfts.cdl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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" ;
Expand Down Expand Up @@ -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 ;

Expand Down
5 changes: 5 additions & 0 deletions parameter_files/fates_params_default.cdl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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" ;
Expand Down Expand Up @@ -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 ;
Expand Down

0 comments on commit 770ac36

Please sign in to comment.