Skip to content

Commit

Permalink
Some minor syntax updates to inverse crown area allometry. Minor opti…
Browse files Browse the repository at this point in the history 
…mization to inverse crown area allometry. Implemented safer logic against a real.
  • Loading branch information
@rgknox
rgknox committed Dec 13, 2018
1 parent e48065b commit 7cdca13
Show file tree
Hide file tree
Showing 2 changed files with 47 additions and 34 deletions.
18 changes: 11 additions & 7 deletions biogeochem/EDCohortDynamicsMod.F90
View file
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
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

0 comments on commit 7cdca13

Please sign in to comment.