Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

change cohort fusion to conserve crown area #447

Merged
merged 6 commits into from
Apr 5, 2019
Merged

change cohort fusion to conserve crown area #447

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
e48065b
changed cohort fusion to conserve crown area
ckoven Nov 17, 2018
7cdca13
Some minor syntax updates to inverse crown area allometry. Minor opti…
rgknox Dec 13, 2018
3c01f12
Merge pull request #2 from rgknox/cohortfusion_crownareaconserv_offma…
ckoven Dec 13, 2018
a722053
merged master and resolved conflicts
ckoven Feb 16, 2019
7d49bef
reincorporated old conservation-of-dbh logic as an option in a case s…
ckoven Feb 16, 2019
36ea2c0
merged master and resolved conflict
ckoven Feb 19, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
100 changes: 82 additions & 18 deletions biogeochem/EDCohortDynamicsMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ module EDCohortDynamicsMod
use FatesConstantsMod , only : r8 => fates_r8
use FatesConstantsMod , only : fates_unset_int
use FatesConstantsMod , only : itrue,ifalse
use FatesConstantsMod , only : fates_unset_r8
use FatesConstantsMod , only : nearzero
use FatesInterfaceMod , only : hlm_days_per_year
use FatesInterfaceMod , only : nleafage
Expand Down Expand Up @@ -98,6 +99,12 @@ module EDCohortDynamicsMod
character(len=*), parameter, private :: sourcefile = &
__FILE__


integer, parameter, private :: conserve_crownarea_and_number_not_dbh = 1
integer, parameter, private :: conserve_dbh_and_number_not_crownarea = 2

integer, parameter, private :: cohort_fusion_conservation_method = conserve_crownarea_and_number_not_dbh

! 10/30/09: Created by Rosie Fisher
!-------------------------------------------------------------------------------------!

Expand Down Expand Up @@ -829,6 +836,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in)
real(r8) :: leaf_c_next ! Leaf carbon * plant density of current (for weighting)
real(r8) :: leaf_c_curr ! Leaf carbon * plant density of next (for weighting)
real(r8) :: dynamic_fusion_tolerance
real(r8) :: dbh
real(r8) :: leaf_c ! leaf carbon [kg]

integer :: largersc, smallersc, sc_i ! indices for tracking the growth flux caused by fusion
Expand Down Expand Up @@ -895,6 +903,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in)
if( currentCohort%isnew.eqv.nextc%isnew ) then

newn = currentCohort%n + nextc%n

fusion_took_place = 1

if ( fuse_debug .and. currentCohort%isnew ) then
Expand All @@ -918,31 +927,86 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in)


! Fuse all mass pools
call currentCohort%prt%WeightedFusePRTVartypes(nextc%prt, currentCohort%n/newn )
call currentCohort%prt%WeightedFusePRTVartypes(nextc%prt, &
currentCohort%n/newn )

currentCohort%laimemory = (currentCohort%n*currentCohort%laimemory &
+ nextc%n*nextc%laimemory)/newn

currentCohort%dbh = (currentCohort%n*currentCohort%dbh &
+ nextc%n*nextc%dbh)/newn

call h_allom(currentCohort%dbh,currentCohort%pft,currentCohort%hite)

currentCohort%canopy_trim = (currentCohort%n*currentCohort%canopy_trim &
+ nextc%n*nextc%canopy_trim)/newn

! -----------------------------------------------------------------
! If fusion pushed structural biomass to be larger than
! the allometric target value derived by diameter, we
! then increase diameter and height until the allometric
! target matches actual bdead. (if it is the other way around
! we then just let the carbon pools grow to fill-out allometry)
! -----------------------------------------------------------------

if( EDPftvarcon_inst%woody(currentCohort%pft) == itrue ) then
call StructureResetOfDH( currentCohort%prt%GetState(struct_organ,all_carbon_elements), currentCohort%pft, &
currentCohort%canopy_trim, currentCohort%dbh, currentCohort%hite )
end if
select case(cohort_fusion_conservation_method)
!
! -----------------------------------------------------------------
! Because cohort fusion is an unavoidable but non-physical process,
! and because of the various nonlinear allometric relationships,
! it isn't possible to simultaneously conserve all of the allometric
! relationships during cohort fusion. We will always conserve carbon,
! but there are choices to made about what else to conserve or not.
! In particular, there is a choice to be made of conservation amongst
! the number density, stem diameter, and crown area. Below,
! some different conservation relationships can be chosen during fusion.
! -----------------------------------------------------------------
!
case(conserve_crownarea_and_number_not_dbh)
!
! -----------------------------------------------------------------
! conserve total crown area during the fusion step, and then calculate
! dbh of the fused cohort as that which conserves both crown area and
! the dbh to crown area allometry. dbh will be updated in the next
! growth step in the (likely) event that dbh to structural iomass
! allometry is exceeded. if using a capped crown area allometry and
! above the cap, then calculate as the weighted average of fusing
! cohorts' dbh
! -----------------------------------------------------------------
!
currentCohort%c_area = currentCohort%c_area + nextc%c_area
!
call carea_allom(dbh,newn,currentSite%spread,currentCohort%pft,&
currentCohort%c_area,inverse=.true.)
!
if (abs(dbh-fates_unset_r8)<nearzero) then
currentCohort%dbh = (currentCohort%n*currentCohort%dbh &
+ nextc%n*nextc%dbh)/newn
else
currentCohort%dbh = dbh
endif
!
call h_allom(currentCohort%dbh,currentCohort%pft,currentCohort%hite)
!
case(conserve_dbh_and_number_not_crownarea)
!
! -----------------------------------------------------------------
! Here we conserve the mean stem diameter of the trees in the cohorts
! rather than the crown area of the cohort
! -----------------------------------------------------------------
!
currentCohort%dbh = (currentCohort%n*currentCohort%dbh &
+ nextc%n*nextc%dbh)/newn
!
call h_allom(currentCohort%dbh,currentCohort%pft,currentCohort%hite)
!
! -----------------------------------------------------------------
! If fusion pushed structural biomass to be larger than
! the allometric target value derived by diameter, we
! then increase diameter and height until the allometric
! target matches actual bdead. (if it is the other way around
! we then just let the carbon pools grow to fill out allometry)
! -----------------------------------------------------------------
!
if( EDPftvarcon_inst%woody(currentCohort%pft) == itrue ) then
call StructureResetOfDH( currentCohort%prt%GetState(struct_organ,all_carbon_elements), currentCohort%pft, &
currentCohort%canopy_trim, currentCohort%dbh, currentCohort%hite )
end if
!
call carea_allom(currentCohort%dbh,newn,currentSite%spread,currentCohort%pft,&
currentCohort%c_area,inverse=.false.)
!
case(default)
write(fates_log(),*) 'FATES: Invalid choice for cohort_fusion_conservation_method'
call endrun(msg=errMsg(sourcefile, __LINE__))
end select

call sizetype_class_index(currentCohort%dbh,currentCohort%pft, &
currentCohort%size_class,currentCohort%size_by_pft_class)
Expand Down
82 changes: 60 additions & 22 deletions biogeochem/FatesAllometryMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,7 @@ module FatesAllometryMod
use FatesConstantsMod, only : cm2_per_m2
use FatesConstantsMod, only : kg_per_Megag
use FatesConstantsMod, only : calloc_abs_error
use FatesConstantsMod, only : fates_unset_r8
use shr_log_mod , only : errMsg => shr_log_errMsg
use FatesGlobals , only : fates_log
use FatesGlobals , only : endrun => fates_endrun
Expand Down Expand Up @@ -449,15 +450,22 @@ end subroutine blmax_allom
! Generic crown area allometry wrapper
! ============================================================================

subroutine carea_allom(d,nplant,site_spread,ipft,c_area)
subroutine carea_allom(dbh,nplant,site_spread,ipft,c_area,inverse)

real(r8),intent(in) :: d ! plant diameter [cm]
real(r8),intent(in) :: site_spread ! site level spread factor (crowdedness)
real(r8),intent(in) :: nplant ! number of plants [1/ha]
integer(i4),intent(in) :: ipft ! PFT index
real(r8),intent(out) :: c_area ! crown area per cohort (m2)

real(r8) :: d_eff ! Effective diameter (cm)
real(r8),intent(inout) :: dbh ! plant diameter at breast (reference) height [cm]
real(r8),intent(in) :: site_spread ! site level spread factor (crowdedness)
real(r8),intent(in) :: nplant ! number of plants [1/ha]
integer(i4),intent(in) :: ipft ! PFT index
real(r8),intent(inout) :: c_area ! crown area per cohort (m2)
logical,optional,intent(in) :: inverse ! if true, calculate dbh from crown area
! instead of crown area from dbh

real(r8) :: dbh_eff ! Effective diameter (cm)
logical :: do_inverse ! local copy of the inverse argument
! defaults to false
logical :: capped_allom ! if we are using an allometry that caps
! crown area at height, we need to make
! special considerations

associate( dbh_maxh => EDPftvarcon_inst%allom_dbh_maxheight(ipft), &
allom_lmode => EDPftvarcon_inst%allom_lmode(ipft), &
Expand All @@ -466,24 +474,49 @@ subroutine carea_allom(d,nplant,site_spread,ipft,c_area)
d2ca_min => EDPftvarcon_inst%allom_d2ca_coefficient_min(ipft), &
d2ca_max => EDPftvarcon_inst%allom_d2ca_coefficient_max(ipft))

if( .not. present(inverse) ) then
do_inverse = .false.
else
do_inverse = inverse
if (do_inverse) then
c_area = c_area / nplant
endif
endif

select case(int(allom_lmode))
case(1)
d_eff = min(d,dbh_maxh)
call carea_2pwr(d_eff,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area)
dbh_eff = min(dbh,dbh_maxh)
call carea_2pwr(dbh_eff,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,do_inverse)
capped_allom = .true.
case(2) ! "2par_pwr")
call carea_2pwr(d,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area)
call carea_2pwr(dbh,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,do_inverse)
capped_allom = .false.
case(3)
d_eff = min(d,dbh_maxh)
call carea_2pwr(d_eff,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area)
dbh_eff = min(dbh,dbh_maxh)
call carea_2pwr(dbh_eff,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,do_inverse)
capped_allom = .true.
case DEFAULT
write(fates_log(),*) 'An undefined leaf allometry was specified: ', &
write(fates_log(),*) 'An undefined leaf allometry was specified: ', &
allom_lmode
write(fates_log(),*) 'Aborting'
call endrun(msg=errMsg(sourcefile, __LINE__))
end select
write(fates_log(),*) 'Aborting'
call endrun(msg=errMsg(sourcefile, __LINE__))
end select


c_area = c_area * nplant
if (capped_allom .and. do_inverse) then
if (dbh_eff .lt. dbh_maxh) then
dbh = dbh_eff
else
! In this situation, we are signaling to the
! calling routine that we we cannot calculate
! dbh from crown area, because we have already
! hit the area cap, and the two are not proportional
! anymore. hopefully, the calling routine has an alternative
dbh = fates_unset_r8
endif
endif

c_area = c_area * nplant

end associate
return
Expand Down Expand Up @@ -1880,20 +1913,21 @@ end subroutine CrownDepth
! =============================================================================


subroutine carea_2pwr(d,spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area)
subroutine carea_2pwr(dbh,spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,inverse)

! ============================================================================
! Calculate area of ground covered by entire cohort. (m2)
! Function of DBH (cm) canopy spread (m/cm) and number of individuals.
! ============================================================================

real(r8),intent(in) :: d ! diameter [cm]
real(r8),intent(inout) :: dbh ! diameter at breast (refernce) height [cm]
real(r8),intent(in) :: spread ! site level relative spread score [0-1]
real(r8),intent(in) :: d2bl_p2 ! parameter 2 in the diameter->bleaf allometry (exponent)
real(r8),intent(in) :: d2bl_ediff ! area difference factor in the diameter-bleaf allometry (exponent)
real(r8),intent(in) :: d2ca_min ! minimum diameter to crown area scaling factor
real(r8),intent(in) :: d2ca_max ! maximum diameter to crown area scaling factor
real(r8),intent(out) :: c_area ! crown area for one plant [m2]
real(r8),intent(inout) :: c_area ! crown area for one plant [m2]
logical,intent(in) :: inverse ! if true, calculate dbh from crown area rather than its reverse

real(r8) :: crown_area_to_dbh_exponent
real(r8) :: spreadterm ! Effective 2bh to crown area scaling factor
Expand All @@ -1917,7 +1951,11 @@ subroutine carea_2pwr(d,spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area)

spreadterm = spread * d2ca_max + (1._r8 - spread) * d2ca_min

c_area = spreadterm * d ** crown_area_to_dbh_exponent
if ( .not. inverse) then
c_area = spreadterm * dbh ** crown_area_to_dbh_exponent
else
dbh = (c_area / spreadterm) ** (1./crown_area_to_dbh_exponent)
endif

end subroutine carea_2pwr

Expand Down
2 changes: 1 addition & 1 deletion main/FatesConstantsMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ module FatesConstantsMod

! Unset and various other 'special' values
integer, parameter :: fates_unset_int = -9999

real(fates_r8), parameter :: fates_unset_r8 = -9999._fates_r8

! Integer equivalent of true (in case some compilers dont auto convert)
integer, parameter :: itrue = 1
Expand Down