From e48065bcd435f7c91b446b10cbdcad750f83df34 Mon Sep 17 00:00:00 2001 From: Charles Koven Date: Fri, 16 Nov 2018 16:28:33 -0800 Subject: [PATCH] changed cohort fusion to conserve crown area --- biogeochem/EDCohortDynamicsMod.F90 | 37 +++++++++++-------- biogeochem/FatesAllometryMod.F90 | 59 ++++++++++++++++++++++-------- main/FatesConstantsMod.F90 | 2 +- 3 files changed, 66 insertions(+), 32 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 71a1416cfc..6bef2c35ad 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -11,6 +11,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 FatesInterfaceMod , only : hlm_days_per_year use EDPftvarcon , only : EDPftvarcon_inst use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type @@ -729,6 +730,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) real(r8) :: newn real(r8) :: diff real(r8) :: dynamic_fusion_tolerance + real(r8) :: dbh integer :: largersc, smallersc, sc_i ! indices for tracking the growth flux caused by fusion real(r8) :: larger_n, smaller_n @@ -794,6 +796,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 @@ -822,26 +825,28 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) 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) - ! ----------------------------------------------------------------- + ! 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 + + currentCohort%c_area = currentCohort%c_area + nextc%c_area + + call carea_allom(dbh,newn,currentSite%spread,currentCohort%pft,& + currentCohort%c_area,inverse=.true.) - 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 + if (dbh .eq. fates_unset_r8) 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) call sizetype_class_index(currentCohort%dbh,currentCohort%pft, & currentCohort%size_class,currentCohort%size_by_pft_class) diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 81881a190a..74f08b6072 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -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 @@ -448,15 +449,17 @@ end subroutine blmax_allom ! Generic crown area allometry wrapper ! ============================================================================ - subroutine carea_allom(d,nplant,site_spread,ipft,c_area) + subroutine carea_allom(d,nplant,site_spread,ipft,c_area,inverse) - real(r8),intent(in) :: d ! plant diameter [cm] + 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(out) :: c_area ! crown area per cohort (m2) + 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 associate( dbh_maxh => EDPftvarcon_inst%allom_dbh_maxheight(ipft), & allom_lmode => EDPftvarcon_inst%allom_lmode(ipft), & @@ -465,24 +468,45 @@ 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 + endif + + if (do_inverse) then + c_area = c_area / nplant + 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) + call carea_2pwr(d_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(d,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) + call carea_2pwr(d_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 (d_eff .lt. dbh_maxh) then + d = d_eff + else + d = fates_unset_r8 + endif + endif + c_area = c_area * nplant end associate return @@ -1848,20 +1872,21 @@ end subroutine h2d_martcano ! ============================================================================= - subroutine carea_2pwr(d,spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area) + subroutine carea_2pwr(d,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) :: d ! diameter [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 @@ -1885,7 +1910,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 * d ** crown_area_to_dbh_exponent + else + d = (c_area / spreadterm) ** (1./crown_area_to_dbh_exponent) + endif end subroutine carea_2pwr diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index f585c53915..eadc1df658 100644 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -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