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

integrate minor updates to crown-conservation #2

Merged
merged 1 commit into from
Dec 13, 2018
Merged

integrate minor updates to crown-conservation #2

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
18 changes: 11 additions & 7 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 : 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 EDPftvarcon , only : EDPftvarcon_inst
use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type
Expand Down Expand Up @@ -820,26 +821,29 @@ 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%canopy_trim = (currentCohort%n*currentCohort%canopy_trim &
+ nextc%n*nextc%canopy_trim)/newn

! 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
! biomass allometry is exceeded. if using a capped crown area allometry and above the cap,
! then calculate as the weighted average of fusing cohorts' 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 (dbh .eq. fates_unset_r8) then
if (abs(dbh-fates_unset_r8)<nearzero) then
currentCohort%dbh = (currentCohort%n*currentCohort%dbh &
+ nextc%n*nextc%dbh)/newn
else
Expand Down
63 changes: 36 additions & 27 deletions biogeochem/FatesAllometryMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -449,17 +449,22 @@ end subroutine blmax_allom
! Generic crown area allometry wrapper
! ============================================================================

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

real(r8),intent(inout) :: 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(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) :: d_eff ! Effective diameter (cm)
logical :: do_inverse, capped_allom
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 @@ -472,23 +477,22 @@ subroutine carea_allom(d,nplant,site_spread,ipft,c_area,inverse)
do_inverse = .false.
else
do_inverse = inverse
endif

if (do_inverse) then
c_area = c_area / nplant
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,do_inverse)
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,do_inverse)
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,do_inverse)
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: ', &
Expand All @@ -499,10 +503,15 @@ subroutine carea_allom(d,nplant,site_spread,ipft,c_area,inverse)


if (capped_allom .and. do_inverse) then
if (d_eff .lt. dbh_maxh) then
d = d_eff
if (dbh_eff .lt. dbh_maxh) then
dbh = dbh_eff
else
d = fates_unset_r8
! 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

Expand Down Expand Up @@ -1872,20 +1881,20 @@ end subroutine h2d_martcano
! =============================================================================


subroutine carea_2pwr(d,spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,inverse)
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(inout) :: 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(inout) :: 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
Expand All @@ -1911,9 +1920,9 @@ subroutine carea_2pwr(d,spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,inver
spreadterm = spread * d2ca_max + (1._r8 - spread) * d2ca_min

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

end subroutine carea_2pwr
Expand Down